数值分析课程设计

1、摘要实验一拉格朗日插及数值求解1.1实验目的了解Lagranger差值的基本原理和方法通过实例掌握用MATLAB求插值的方法根据实际计算理论，利用Lagranger插值多项式计算1.2实验原理设已知X。，XiX2,.,Xn及y.f（Xi）（i=0,i,.,n）,Ln（x）为不超过n次多项式且满足Ln(Xi)=yi(i=0,1,.n)易知Ln（X）v|0（X）y。.ln（X）yn其中，li（X）均为n次多项式，再由Xj（j*i）为n次多项式1i（X）的n个根知nli(X)=C|IX-Xjj=0.最后，由nli(Xj)=C|1(Xi-Xj)=1=jfj号1nI1(Xi-Xj)jfje,i=0,1,

2、.,n.nlOi总之，Ln(X)=i=0nX-Xj二.j=0X-Xjli(X)=j#式为n阶Lagrange插值公式，其中，li（X）（i=0,1,.n）称为n阶Lagrange插值的基函数(X-X0).(X-Xi4)(X-Xi1).(X-Xn)li(X)=,i(X-M).(Xi-Xy)(Xi-Xi1).(Xi-Xn)=0,12.,n1.3实验内容functiony=lagranger(X0,y0,X);%UNTITLEDSummaryofthisfunctiongoeshere%Detailedexplanationgoesheren=length(x0);m=length(x);fori=

3、1:mz=x(i);s=0.0;fork=1:nli=1.0;forj=1:nifj=kli=li*(z-x0(j)/(x0(k)-x0(j);endends=li*y0(k)+s;endy(i)=s;end1.4实验案例及结果分析(1)输入：x0=4,5,6;y0=10,5.25,1;x=5;y=lagranger(x0,y0,x)y-5. 2500输入：X0=1,4,8;y0=6,3.2,4;x=4;y=lagranger(x0,y0,x)y=实验二LU分解法解线性方程组2.1 实验目的1 .了解LU分解法解线性方程组的基本原理；2 .熟悉计算方法的技巧和过程，能用LU分解法解实际问题;3

4、 .用matlab实现LU分解。2.2 实验原理1 .若一个线性方程组系数矩阵为n阶方阵A且各阶顺序主子式均不为0则A的LU分解存在且唯一。将高斯消去法改写为紧凑形式，可以直接从矩阵A的元素得到计算L,U元素的递推公式，而不需任何中间步骤，这就是所谓直接三角分解法。一旦实现了矩阵A的LU分解，那么求解Ax=b的问题就等价于求解两个三角形方程组：Ly=b，求y;Ux=y，求x。2 .在满足1的条件下课推导得出以下公式j1j4uij-ajiljkukik1i4(2)yi刊-xlikykk1nxi=(yi-'Uikxk)/Uiiki1(1)lij=(aij-工likuij)/uijk13 .

5、公式（1）用于求解矩阵L、U,公式（2）用于会带求解V、X。从公式中可以看出：L对角线上元素为1,U第一行与A第一行相同。4 .LU分解的具体过程和顺序如下：(1)第一步分解：51=41(2)第二步分解：l21=a21/U11U12=a12U22-a22-l21U12(3)第三步分解：I31=a31/U11l32=(a32-I31U12)/U22U13=a13U23=a23-I21U13U33=a33-I31U13-I32U23(0)第n步分解：依次计算：In1、In2Inn,UinUnn1"U11U12.Um11211U22U2n12.3 实验内容编写一个M文件fUnctionL,

6、U,x=LU_x(A,b)n,m=size(A);ifn-=merror('TherowsandcoIUmnsofmatrixAmUStbeeqUaI!');retUrn;endfori=1:iiforj=1:iiAA(i,j)=A(i,j);endendif(det(AA)=0)error('ThematrixcannotbedividedbyLU!')return;endendAn,n=size(A);L=zeros(n,n);U=zeros(n,n);fori=1:nL(i,i)=1;endfork=1:nforj=k:nU(k,j)=A(k,j)-sum

7、(L(k,1:k-1).*U(1:k-1,j)');endfori=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)')/U(k,k);endendy(1)=d(1);fori=2:nforj=1:i-1d(i)=d(i)-L(i,j)*y(j);endy(i)=d(i)；endx(n)=y(n)/U(n,n);fori=(n-1):-1:1forj=n:-1:i+1y(i)=y(i)-U(i,j)*x(j);endx(i)=y(i)/U(i,i);end2.4 实验案例及结果分析一1070-32.0999996515101反12x2

8、2一x4x3一815.9000015-11在MATLABt令窗体输入:>>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)得到结果如下：10,0000-7.000001.0000-3.00002.10006.OOCO2.00005.DOOO-1.00005.0000700002.00001.00000之00001.0e+006*0.0000000一口，0000口.0000000.0000-2.50000.000000.0000-2.4C000.

9、OOCO0.0000一L0e+007*0,0000-0.000000.00000-a.ocoo0.OOCO0.CIOOO001.50000.5750a00.aaoo-0,0098-1.01181.01061.0157实验三龙贝格求积公式求数值积分3.1实验目的1 .熟练掌握龙贝格求积的基本思路和步骤;2 .培养编程与上机调试能力；3 .利用龙贝格求积方法求解积分。3.2 实验原理(1),b-ah=置n=1,精度要求3n(2)计算TT%).fMh2n=二hn(3)置2,并计算(0)1(0)n12?=Tn()h2nXf(a(2k-1')232n9km(2Tk49k(mT-计5mm二一(6)

10、若m=14专(7)；否贝U,置2,k=k+1,转(5)。T(k)_T(kJ1)(k)(7)若11,则停止计算(输出T1),否则转(3)3.3 实验内容functionR=romberg(f,a,b,e)%参数介绍：%f-被积函数f(x)%a-x的左区间.%b-x的右区间.%e-误差限.%结果：%R-返回Romberg表.n=1;%区间二分次数while1%在此仅代表多次二分，在后面判断循环终止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点梯形公式）.fori=1:n%按照公式

11、计算Romberg表的第一列.h=(b-a)/2Ai;s=0;fork=1:2A(i-1)s=s+feval(f,a+(2*k-1)*h);endR(i+1,0+1)=R(i-1+1,0+1)/2+h*s;endforj=1:n%计算Romberg表的其他列.fac=1/(4Aj-1);form=j:nR(m+1,j+1)=R(m+1,j-1+1)+fac*(R(m+1,j-1+1)-R(m-1+1,j-1+1);endendifabs(R(n,n)-R(n+1,n+1)<e%当精度满足设定要求时退出break;endn=n+1;%未达到指定精度继续二分endt=R(i,j)3.4 实验

12、案例及结果分析15.dx用Romberg求积法计算积分x+3,取精度要求=10-50在Matlab命令窗口输入:fun=inline('5./(3+x)','x');romberg(fun,-1,1,1e-6)输出结果如下：3.4657ans=3.75000000D354173.472200003.48513.46533.46590003.47063.46583.46573.4657003.46703,46573.46573.46573.465703.46603,46573.46"3.46573.46573.4657由输出结果可知最终积分为3.4657

13、。实验四Runge-Kutta方法求常微分方程数值解4.1 实验目的1 .熟悉Runge-Kutta常微分方程初值问题的基本原理2 .了解Runge-Kutta常微分方程初值问题的计算流程3 .能编程实现Runge-Kutta常微分方程初值问题4.2 实验原理在欧才立法中(Xn,yn,h)=f(Xn,yn)的基础上增加了计算一个右函数f的值，可望p=2。若要使得到的公式阶数p更大，中就必须包含更多的f值。实际上从与方程等价的积分形式即X4n1y(Xn1-yn)=xf(X,y(x)dxxn若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必须要增加求积节点，为此可将上述公式的右端项表

14、示为Xn1rXf(x,y(x)dx之h£Cif(Xn+九由，y(Xn+儿h)XnT一般来说，点数r越多，精度越高，上式右端相当于增量函数中(X，y，h),为得到便于计算的显示方法，可类似改进的欧拉法，将公式表示为yn1=ynh(Xn,yn,h)r(Xn,yn,h)=.cKii1iKi=f(Xnih,ynh'ijKj),i=2,.,r这里G,A,均为常数，上式称为显示龙格一库塔(Ruuge-Kutta)法，简称R-K方法。当r=1时，就是欧拉法，当r=2时，改进的欧拉法就是其中的一种。yn1=ynh*(K1七)/2Ki=f(Xn,yn)K2=f(Xnh,ynh*K1)依次类推，

15、如果在区间(为凶2内多预估几个点上的斜率值KK2、.Km,并用他们的加权平均数作为平均斜率K"的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格-库塔公式，也就是在工程中应用广泛的经典龙格-库塔算法：hyni=yn(Ki2K22K3K4)6ynyy-1Kh-2+h-2+2KJI2Xnhhyn2K2(xnh)ynhK3hK3=ynK22K4=ynhK34.3实验内容四阶龙格一库塔法的计算公式为：Ki=g(Xi,yi)hhK2=g(Xi二，yi-Ki)22hhK3=g(Xi,yiK2)22K4=g(Xih,yhK3)hyii=yiKi2K22K3K4)四

16、阶龙格一库塔公式的Matlab程序代码：functiony=DELGKT4_lungkuta(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(i)=y0;x=a:h:b;var=findsym(f);fori=2:N+iK1=Funval(f,varvec,x(i-1)y(i-1);K2=Funval(f,varvec,x(i-1)+h/2y(i-1)+K1*h/2);K3=Funval(f,varvec,x(i-1)+h/2y(i-1)+K2*h/2);K4=Funval(f,varvec,x(i-1)+hy(i-1)+h*K3

17、);y(i)=y(i-1)+h*(K1+2*K2+2*K3+K4)/6;end函数运行时需要调用下列函数：functionfv=Funval(f,varvec,varval)var=findsym(f);iflength(var)<4ifvar(1)=varvec(1)fv=subs(f,varvec(1),varval(1);elsefv=subs(f,varvec(2),varval(2);endelsefv=subs(f,varvec,varval);end4.4实验案例及结果分析四阶龙格一库塔求解一阶常微分方程应用实例，用四阶龙格一库塔法求下面常微分方程的数值解。J"y

18、3y=10MxM1在MATLABt令窗口输入:symsxy;z=4*x+y+3yy=DELGKT3_kuta(z,0.1,0,1,1,xy)输出结果如下：yy=L00000000000000L44133333333333L97114688B888892,598745836703703,334413940530384.18951S473276155.176309587856725137.603779512990729.0761436584402510.74541809935288实验五用插值方法求解带式录音机的播放时间问题5.1 实验目的进一步巩固、加强插值模型的建模、

19、求解能力。初步研究插值法的稳定性和数值解法。学习掌握用MATLAB软件求解插值的相关命令。5.2 实验原理1 .用interp1可以绘制出曲线；2 .用8次多项式做插值，求解8个系数，从而可以显示给出一个函数5.3 实验内容>>猿项式插值>>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,

20、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(原始数据','最近点插值,线性插值','样条插值','立方插值')输出结果如下：&

21、gt;>核项式拟合>>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:Polynomialisnotunique;degree>=numberofdatapoints.> Inpolyfitat720.0000-0.00000.0000-0.00480.8068000> >c1=linspace(0,10,1000);%强由范围> >t1=polyval(p,c1);%5+算在c1数据点的多项式值> >plot(c,t,'*',c,t,c1,t1,':b&