数值分析课程设计

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

2、X-Xj二.j=0X-Xjli(X)=j#式为n阶Lagrange插值公式,其中,liXi=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=1:mz=x(i);s=0.0;fork=1:nl

3、i=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 .用matlab实现LU分解.2.2 实验原理1

4、 .假设一个线性方程组系数矩阵为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 .公式1用于求解矩阵L、U,公式2用于会带求解V、X.从公

5、式中可以看出：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,UinUnn1U11U12.Um11211U22U2n12.3 实验内容编写一个M文件fUnctionL,U,x=LU_x(A,b)n,m=size(A);ifn-=merror(T

6、herowsandcoIUmnsofmatrixAmUStbeeqUaI!);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(L(k,1:k-1).*U(1:k-1,j);endfori=k+1:nL(i,k)=(A(i,k)-sum(L(

7、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反12x22一x4x3一815.9000015-11在MATLABt令窗体输入:A=10,-7,0,1;-3,2.099999,6,2;5,-1,5,

8、-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.OOCO0.0000一L0e+007*0,0000-0.000000.00000-a.ocoo0.OOCO0.CIOOO001.50000.5750a00.aaoo-0,0098-1.01

9、181.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)假设m=14专(7)；否贝U,置2,k=k+1,转(5).T(k)_T(kJ1)(k)(7)假设11,那么停止计算(输出T1),否那么转(3)3.3 实验内容functionR=romberg(f,

10、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%根据公式计算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)/

11、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 实验案例及结果分析15.dx用Romberg求积法计算积分x+3,取精度要求=10-50在Matlab命令窗口输入:fun=inline(5./(3+x),x);romberg(fun,-1,1,1e-6)

12、输出结果如下：3.4657ans=3.75000000D354173.472200003.48513.46533.46590003.47063.46583.46573.4657003.46703,46573.46573.46573.465703.46603,46573.463.46573.46573.4657由输出结果可知最终积分为3.4657.实验四Runge-Kutta方法求常微分方程数值解4.1 实验目的1 .熟悉Runge-Kutta常微分方程初值问题的根本原理2 .了解Runge-Kutta常微分方程初值问题的计算流程3 .能编程实现Runge-Kutta常微分方程初值问题4.2 实

13、验原理在欧才立法中(Xn,yn,h)=f(Xn,yn)的根底上增加了计算一个右函数f的值,可望p=2.假设要使得到的公式阶数p更大,中就必须包含更多的f值.实际上从与方程等价的积分形式即X4n1y(Xn1-yn)=xf(X,y(x)dxxn假设要使公式阶数提升,就必须使右端积分的数值求积公式精度提升,它必须要增加求积节点,为此可将上述公式的右端项表示为Xn1rXf(x,y(x)dx之hCif(Xn+九由,y(Xn+儿h)XnT一般来说,点数r越多,精度越高,上式右端相当于增量函数中(X,y,h),为得到便于计算的显示方法,可类似改良的欧拉法,将公式表示为yn1=ynh(Xn,yn,h)r(Xn

14、,yn,h)=.cKii1iKi=f(Xnih,ynhijKj),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)依次类推,如果在区间(为凶2内多预估几个点上的斜率值KK2、.Km,并用他们的加权平均数作为平均斜率K的近似值,显然能构造出具有很高精度的高阶计算公式.经数学推导、求解,可以得出四阶龙格-库塔公式,也就是在工程中应用广泛的经典龙格-库塔算法：hyni=yn(Ki2K22K3K4)6

15、ynyy-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)四阶龙格一库塔公式的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+i

16、K1=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);y(i)=y(i-1)+h*(K1+2*K2+2*K3+K4)/6;end函数运行时需要调用以下函数：functionfv=Funval(f,varvec,varval)var=findsym(f);iflength(var)猿项式插值t=205,430,677,945,1233

17、,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=

18、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) xlabel(c); ylabel(t=t(c); title(八次曲线拟合)图2八次多项式拟合图通过实验2果可