数值分析课程设计

摘要实验一拉格朗日插及数值求解1.1实验目的了解Lagranger差值的基本原理和方法通过实例掌握用MATLAB求插值的方法根据实际计算理论，利用Lagranger插值多项式计算1.2实验原理设已知0x，1x，2x，...,nx及iy=f(ix)(i=0,1,.....,n),)(xLn为不超过n次多项式且满足iinyxL)(（i=0,1,...n）.易知nnnyxlyxlxL)(...)()(00其中，)(xli均为n次多项式,再由jx（ji）为n次多项式)(xli的n个根知niijjixxcxl0)(.最后，由1)()(0nijjjijixxcxlnijjjixxc0)(1,i=0,1,...,n.总之，)(xLn=iniiyxl0)(，)(xli=.0nijjjijxxxx式为n阶Lagrange插值公式，其中，)(xli（i=0,1,...n）称为n阶Lagrange插值的基函数。nixxxxxxxxxxxxxxxxxlniiiiiiniii...,2,1,0,))...()()...(())...()()...(()(1101101.3实验内容functiony=lagranger(x0,y0,x);%UNTITLEDSummaryofthisfunctiongoeshere%Detailedexplanationgoesheren=length(x0);m=length(x);fori=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)（2）输入：X0=[1,4,8];y0=[6,3.2,4];x=4;y=lagranger(x0,y0,x)实验二LU分解法解线性方程组2.1实验目的1.了解LU分解法解线性方程组的基本原理；2.熟悉计算方法的技巧和过程，能用LU分解法解实际问题；3.用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的条件下课推导得出以下公式（1）ijjkijikijijuulal/)(1111jkkijkjiijulau（2）11ikkikiiylbyiinikkikiiuxuyx/)(13.公式（1）用于求解矩阵L、U，公式（2）用于会带求解y、x。从公式中可以看出：L对角线上元素为1，U第一行与A第一行相同。4．LU分解的具体过程和顺序如下：（1）第一步分解：1111au（2）第二步分解：112121/ual1212au12212222ulau（3）第三步分解：113131/ual2212313232/)(uulal1313au13212323ulau233213313333ululau（n）第n步分解：依次计算：1nl、2nl......1nnl，inu......nnunnnnnnuuuuuulll....1.1.1.12221121121212.3实验内容编写一个M文件function[L,U,x]=Lu_x(A,b)[n,m]=size(A);ifn~=merror('TherowsandcolumnsofmatrixAmustbeequal!');return;endforii=1:nfori=1:iiforj=1:iiAA(i,j)=A(i,j);endendif(det(AA)==0)error('ThematrixcannotbedividedbyLU!')return;endendA[n,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(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实验案例及结果分析15900001.582012151526099999.23107104321xxxx在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)得到结果如下：实验三龙贝格求积公式求数值积分3.1实验目的1.熟练掌握龙贝格求积的基本思路和步骤；2.培养编程与上机调试能力；3.利用龙贝格求积方法求解积分。3.2实验原理（1）置n=1，精度要求，nabh（2）计算)()(21)0(1bfafhT（3）置nnhh212，并计算nknnnnhkafhTT122)0()0(2))12((21（4）置m=n，n=2n，k=1。（5）计算144)1()1(2)(kkmkmkkmTTT。（6）若m=1，转（7）；否则，置2mm，k=k+1，转（5）。（7）若)1(1)(1kkTT，则停止计算（输出)(1kT），否则转（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%按照公式计算Romberg表的第一列.h=(b-a)/2^i;s=0;fork=1:2^(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/(4^j-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实验案例及结果分析用Romberg求积法计算积分dxx1135，取精度要求=10-5。在Matlab命令窗口输入：fun=inline('5./(3+x)','x');romberg(fun,-1,1,1e-6)输出结果如下：由输出结果可知最终积分为3.4657。实验四Runge-Kutta方法求常微分方程数值解4.1实验目的1.熟悉Runge-Kutta常微分方程初值问题的基本原理2.了解Runge-Kutta常微分方程初值问题的计算流程3.能编程实现Runge-Kutta常微分方程初值问题4.2实验原理在欧拉法),(),,(nnnnyxfhyx的基础上增加了计算一个右函数f的值，可望p=2。若要使得到的公式阶数p更大，就必须包含更多的f值。实际上从与方程等价的积分形式即dxxyxfyxynnxxnn1))(,()(1若要使公式阶数提高，就必须使右端积分的数值求积公式精度提高，它必须要增加求积节点，为此可将上述公式的右端项表示为))(,())(,(11hxyhxfchdxxyxfininxxriinn一般来说，点数r越多，精度越高，上式右端相当于增量函数),,(hyx，为得到便于计算的显示方法，可类似改进的欧拉法，将公式表示为),,(1hyxhyynnnniriinnKchyx1),,(riKhyhxfKijjijnini,...,2),,(1这里ijiic,,均为常数，上式称为显示龙格—库塔（Ruuge-Kutta）法，简称R-K方法。当r=1时，就是欧拉法，当r=2时，改进的欧拉法就是其中的一种。2/)(*211KKhyynn),(1nnyxfK)*,(12KhyxfKnhn依次类推，如果在区间),(1iixx内多预估几个点上的斜率值,...21mKKK、、并用他们的加权平均数作为平均斜率K的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法：334223112143211)(22222222)22(6hKyhxhKyKKhyhxKhyKKhyhxKhyKyxyKKKKKhyynnnnnnnnnnnnnn4.3实验内容四阶龙格—库塔法的计算公式为：)22(6),()2,2()2,2(),(432113423121KKKKhyyhKyhxgKKhyhxgKKhyhxgKyxgKiiiiiiiiii四阶龙格—库塔公式的Matlab程序代码：functiony=DELGKT4_lungkuta(f,h,a,b,y0,varvec)formatlong;N=(b-a)/h;y=zeros(N+1,1);y(1)=y0;x=a:h:b;var=findsym(f);fori=2:N+1K1=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)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实验案例及结果分析四阶龙格—库塔求解一阶常微分方程应用实例，用四阶龙格—库塔法求下面常微分方程的数值解。1)0(34yyxdxdy10x在MATLAB命令窗口输入：symsxy;z=4*x+y+3yy=DELGKT3_kuta(z,0.1,