数值分析报告课程设计.doc
word实验一拉格朗日插与数值求解实验目的了解Lagranger差值的根本原理和方法通过实例掌握用MATLAB求插值的方法根据实际计算理论，利用Lagranger插值多项式计算实验原理设，.,与=f()(i=0,1,.,n),为不超过n次多项式且满足i=0,1,.n.易知其中，均为n次多项式,再由ji为n次多项式的n个根知.最后，由,i=0,1,.,n.总之，=，=式为n阶Lagrange插值公式，其中，i=0,1,.n称为n阶Lagrange插值的基函数。实验内容function y = lagranger(x0,y0,x);%UNTITLED Summary of this function goes here% Detailed explanation goes heren=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n li=1.0; for j=1:n if j=k li=li*(z-x0(j)/(x0(k)-x0(j); end end s=li*y0(k)+s; end y(i)=s;end实验案例与结果分析1输入：x0=4,5,6;y0=10,5.25,1;x=5;y=lagranger(x0,y0,x)2输入：X0=1,4,8;y0=6,3.2,4;x=4;y=lagranger(x0,y0,x)实验二LU分解法解线性方程组实验目的1.了解LU分解法解线性方程组的根本原理；2.熟悉计算方法的技巧和过程，能用LU分解法解实际问题；3.用matlab实现LU分解。实验原理1.假如一个线性方程组系数矩阵为n阶方阵A且各阶顺序主子式均不为0如此A的LU分解存在且唯一。将高斯消去法改写为紧凑形式，可以直接从矩阵A的元素得到计算L，U元素的递推公式，而不需任何中间步骤，这就是所谓直接三角分解法。一旦实现了矩阵A的LU分解，那么求解Ax=b的问题就等价于求解两个三角形方程组：Ly=b，求y；Ux=y，求x。2.在满足1的条件下课推导得出以下公式123.公式1用于求解矩阵L、U，公式2用于会带求解y、x。从公式中可以看出：L对角线上元素为1，U第一行与A第一行一样。4LU分解的具体过程和顺序如下：1第一步分解：2第二步分解：3第三步分解：n第n步分解：依次计算：、.，.实验内容编写一个M文件function L,U,x=Lu_x(A,b)n,m=size(A);if n=merror('The rows and columns of matrix A must be equal!');return;endfor ii=1:nfor i=1:iifor j=1:iiAA(i,j)=A(i,j);endendif (det(AA)=0)error('The matrix can not be divided by LU!')return;endendA n,n=size(A); L=zeros(n,n); U=zeros(n,n); for i=1:n L(i,i)=1; end for k=1:n for j=k:n U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)'); end for i=k+1:n L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)')/U(k,k); end end y(1)=d(1);for i=2:n for j=1:i-1d(i)=d(i)-L(i,j)*y(j);endy(i)=d(i);endx(n)=y(n)/U(n,n);for i=(n-1):-1:1for j=n:-1:i+1y(i)=y(i)-U(i,j)*x(j);endx(i)=y(i)/U(i,i);end实验案例与结果分析在MATLAB命令窗体输入：>>A=10,-7,0,1;-3,2.099999,6,2;5,-1,5,-1;2,1,0,2;>>b=8,5.900001,5,1;>> L,U,x=Lu_x(A,b)得到结果如下：实验三龙贝格求积公式求数值积分实验目的1. 熟练掌握龙贝格求积的根本思路和步骤；2. 培养编程与上机调试能力；3.利用龙贝格求积方法求解积分。实验原理1置n=1，精度要求，2计算3置，并计算4置m=n，n=2n，k=1。5计算。6假如m=1，转7；否如此，置，k=k+1，转5。7假如，如此停止计算输出，否如此转3。实验内容function R = romberg(f, a, b, e)% 参数介绍:% f - 被积函数f(x)% a - x的左区间.% b - x的右区间.% e - 误差限.% 结果:% R - 返回Romberg 表. n = 1; %区间二分次数while 1 %在此仅代表屡次二分，在后面判断循环终止 R = zeros(n + 1, n + 1);%生成(n+1)*(n+1)的0矩阵 R(0+1, 0+1) = (b - a) / 2 * (feval(f, a) + feval(f, b); % 初始值(2点梯形公式). for i = 1 : n % 按照公式计算Romberg 表的第一列. h = (b - a) / 2i; s = 0; for k = 1 : 2(i-1) s = s + feval(f, a + (2*k - 1)*h); end R(i+1, 0+1) = R(i-1+1, 0+1)/2 + h*s; end for j = 1 : n % 计算Romberg 表的其他列. fac = 1 / (4j - 1); for m = j : n R(m+1, j+1) = R(m+1, j-1+1) + fac*(R(m+1, j-1+1) - R(m-1+1, j-1+1); end end if abs(R(n,n) - R(n+1,n+1) < e %当精度满足设定要求时退出 break; end n = n + 1; %未达到指定精度继续二分endt=R(i,j)实验案例与结果分析用Romberg求积法计算积分，取精度要求=10-5。在Matlab命令窗口输入：fun=inline('5./(3+x)','x') ;romberg(fun,-1,1,1e-6)输出结果如下：由输出结果可知最终积分为。实验四 Runge-Kutta方法求常微分方程数值解实验目的1.熟悉Runge-Kutta常微分方程初值问题的根本原理2.了解Runge-Kutta常微分方程初值问题的计算流程3.能编程实现Runge-Kutta常微分方程初值问题实验原理在欧拉法的根底上增加了计算一个右函数f的值，可望p=2。假如要使得到的公式阶数p更大，就必须包含更多的f值。实际上从与方程等价的积分形式即假如要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必须要增加求积节点，为此可将上述公式的右端项表示为一般来说，点数r越多，精度越高，上式右端相当于增量函数，为得到便于计算的显示方法，可类似改良的欧拉法，将公式表示为这里均为常数，上式称为显示龙格库塔Ruuge-Kutta法，简称R-K方法。当r=1时，就是欧拉法，当r=2时，改良的欧拉法就是其中的一种。依次类推，如果在区间内多预估几个点上的斜率值并用他们的加权平均数作为平均斜率的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格库塔公式，也就是在工程中应用广泛的经典龙格库塔算法：实验内容四阶龙格库塔法的计算公式为：四阶龙格库塔公式的Matlab程序代码：function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec)format long;N = (b-a)/h;y = zeros(N+1,1);y(1) = y0;x = a:h:b;var = findsym(f);for i=2:N+1 K1 = Funval(f,varvec,x(i-1) y(i-1); K2 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K1*h/2); K3 = Funval(f,varvec,x(i-1)+h/2 y(i-1)+K2*h/2); K4 = Funval(f,varvec,x(i-1)+h y(i-1)+h*K3); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6;end函数运行时需要调用如下函数： function fv=Funval(f, varvec, varval) var= findsym(f); if length(var)<4 if var(1)=varvec(1) fv=subs(f,varvec(1),varval(1); else fv=subs(f,varvec(2),varval(2); end else fv=subs(f,varvec,varval); end实验案例与结果分析四阶龙格库塔求解一阶常微分方程应用实例，用四阶龙格库塔法求下面常微分方程的数值解。在MATLAB命令窗口输入：syms x y;z=4*x+y+3yy=DELGKT3_kuta(z,0.1,0,1,1,x y)输出结果如下：实验五用插值方法求解带式录音机的播放时间问题实验目的进一步巩固、加强插值模型的建模、求解能力。初步研究插值法的稳定性和数值解法。学习掌握用MATLAB软件求解插值的相关命令。实验原理1.用interp1可以绘制出曲线；2.用8次多项式做插值，求解8个系数，从而可以显示给出一个函数。实验内容>> %多项式插值>> t=205,430,677,945,1233,1542,1872,2224;>> c=100,200,300,400,500,600,700,800;>> c1=0:10:1000;>> t1=interp1(c,t,c1,'nearest');>> t2=interp1(c,t,c1,'linear');>> t3=interp1(c,t,c1,'spline');>> t4=interp1(c,t,c1,'cubic');>> plot(c,t,'*',c1,t1,':b',c1,t2,'-r',c1,t3,'-g',c1,t4,'.-r')>> legend('原始数据','最近点插值','线性插值','样条插值','立方插值')输出结果如下：>> %多项式拟合>> t=205,430,677,945,1233,1542,1872,2224;>> c=100,200,300,400,500,600,700,800;>> n=8; %拟合度>> p=polyfit(c,t,n)Warning: Polynomial is not unique; degree >= number of data points.> In polyfit at 72p = 1.0e-006 * 0.0000 -0.0000 0.0000 -0.0048 0.8068 0 0 0 0>> c1=linspace(0,10,1000);%x轴X围>> t1=polyval(p,c1);%计算在c1数据点的多项式值>> plot(c,t,'*',c,t,c1,t1,':b')>> xlabel('c');>> ylabel('t=t(c)');>> title('八次曲线拟合')图2 八次多项式拟合图通过实验结果可知:(1)播放时间随计数器读数增大而增大；(2)与计数器读数具有近似线性关系；总结通过本次实验，我可以利用matlab的一些根本功能对问题进展求解，初步利用Matlab进展实验操作，通过此次实验也使我了解到MATLAB数学软件功能的强大，学会了用MATLAB进展编程，解决数值分析中的实际问题，同时在编程过程中，更加深入的了解了LU分解法，拉格朗日插值，龙贝格求积，龙格库塔的原理，对算法根本原理认识深刻，进一步加深相关算法的效用和使用X围的全面了解。没有像书中那样进展枯燥的定理论述，而是重点介绍了方法根本思想以与在MATLAB平台上的使用。课堂气氛更加活跃，使我们在不知不觉中掌握了大局部知识。还学会了设计算法不单单是设计出来，一个好的算法就要保证算法具有收敛性和较高的收敛速度，保证算法具有数值的稳定性，小心处理病态问题，注意影响收敛速度的细节。此次实验，使我进一步学习了数值分析解法，了解插值的本质，利用插值基函数求解函数。利用龙贝格求解积分的数值解，就使得原来复杂的求积过程变得简单化。总之通过本次实验，我学到了很多，加深了数值分析的理解。参考文献【1】苏金明阮沈勇，实用指南下册，：电子工业，2002【2】关治数值方法清华大学，2006【3】白峰杉数值计算引论高等教育，2004【4】李庆扬，王能超，易大义数值分析第5版清华大学，201017 / 17
