数值分析实验报告

1、图数值分析实验报告姓名：张献鹏学号：173511038专业：冶金工程班级：重冶二班1 拉格朗日插值11.1 问题背景11.2 数学模型11.3 计算方法11.4 数值分析22复化辛普森求积公式21.1 问题背景21.2 数学模型31.3 计算方法31.4 数值分析53 矩阵的LU分解53.1 问题背景53.2 数学模型63.2.1 理论基础63.2.2 实例63.3 计算方法63.4 小组元的误差84 二分法求方程的根94.1 问题背景94.2 数学模型94.3 计算方法94.4 二分法的收敛性105雅可比迭代求解方程组101.1 问题背景101.2 数学模型111.2.1 理论基础11111

2、.2.2 实例.1.3 计算方法111.4 收敛性分析136 Romber球积法136.1 问题背景136.2 数学模型：136.2.1 理论基础136.2.2 实例146.3 计算方法146.4 误差分析157 幕法157.1 问题背景157.2 数学模型167.2.1 理论基础167.2.2 实例167.3 计算方法167.4 误差分析178 改进欧拉法178.1 问题背景178.2 数学模型188.2.1 理论基础188.2.2 实例188.3 数学模型188.4 误差分析191拉格朗日插值1.1 问题背景1f(x)2对于函数1x,5X5求拉格朗日插值。n10,把f(x)和插值多项式的曲

3、线画在同一张图上进行比较，观察数值积分中的Lagrange插值。1.2 数学模型?= E?= 0取等距差值节点？?=-5+10?/n,?=0,1,.,n,构造n次lagrange插值多项?+1(?)1+?(?-?2?)?+1(?)当n=10时，十次插值多项式Lo(x)以及函数f(x)的图像可以由Matlab画出。1.3 计算方法f.m:functionf=f(x)f=1./(1+x.A2);endLagrange.mfunctiony=Lagrange(x0,y0,x);n=length(x0);m=length(x);fori=1:mz=x(i);s=0.0;fork=1:np=1.0;fo

4、rj=1:nifj=kp=p*(z-x0(j)/(x0(k)-x0(j);endends=p*y0(k)+s;endy(i)=s；End拉格朗日插值的曲线：x=-5:1:5;y=1./(1+x.A2);x0=-5:0.001:5;y0=Lagrange(x,y,x0);y1=1./(1+x0.A2);plot(x0,y0,'b')holdonplot(x0,y1,'r')运行这个文件可以得到如下拉格朗日图形曲线:1.4 数值分析Lio(x)的断误差Rio(x)=f(x)-Lio(x)在区间卜5,5的两端非常大。例如，Lio(4.8)=1.8O438，而f(4.8

5、)=0.04160。这种现象称之为龙格现象。不管n取多大，Runge现象依然存在。因此，对函数作插值多项式时，必须小心处理，不能认为差值节点取得越多，差值余项就越小。此外，当节点增多时，舍入误差的影响不能低估。为了克服高次插值的不足，应采用分段低次插值。2复化辛普森求积公式2.1 问题背景用复化Simpson公式计算定积分？?2=/?的近似值，要求误差限己=1/2X10-7,利用其余项对算法做出步长的事前估计；并将计算结果与精确值进行比较。2.2 数学模型将积分区间a,b分为n等分，h=(b-a)/n,Xk=a+kh,k=0,1,-no在每个子区问Xk,Xk+i上用Simpson公式可得：“？

6、-1?-i?f?(x)dx=E?(?)?=太汇?(?)+4?(?+i)+?(?+i)?=0?=02其中Xk+1/2=Xk+1/2h。?-1?Sn?=yE?(?)+4?(?+1)+?(?+1)=-g-?(?)6?=02?-1?-1+4汇?(?+1)+2汇?(?)+?(?)?=02?=1此式即为复化Simpson公式。设f(X)C4a,b,由Simpson公式的误差有5?=?-?=Z?=0-9q(y)?(4)(?),?,?+1o则复化Simpson公式的余项为：一?-?一4一(4)?=-880?()(?),?£?,?复化Simpson公式四阶收敛。2.3 计算方法程序1(求f(X)的n阶

7、导数)：symsxf=X*eXp(X)%定义函数f(x)n=input('输入所求导数阶数:')f2=diff(f,x,n)%tf(x)的n阶导数程序2:clcclearsymsx%!义自变量xf=inline('x*exp(x)','x')%!义函数f(x)=x*exp(x),换函数时只需换该函数表达式即可f2=inline( '(4*exp(x) + x*exp(x)','x')%定义f(x)的四阶导数，输入程序1里求出的f2即可f3= '-(4*exp(x) + x*exp(x)'表达式带引号

8、，故取负号，一边求最大值e=5*10A(-8)a=1瞅分下限b=2瞅分上限x1=fminbnd(f3,1,2)大值点对应的x值for n=2:10000009因fminbnd()函数求的是表达式的最小值，且要求嘛度要求值%R负的四阶导数的最小值点，也就是求四阶导数的最说等分数nRn=-(b-a)/180*(b-a)/(2*n)F4*f2(x1)ifabs(Rn)<ebreakendendh=(b-a)/nSn1=0Sn2=0fork=0:n-1xk=a+k*hxk1=xk+h/2Sn1=Sn1+f(xk1)Sn2=Sn2+f(xk)end府算余项%用余项进行判断%符合要求时结束%求两组连

9、加和Sn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a)+f(b)%USn2力口了k=0时的值，故减去f(a)z=exp(2)R=Sn-z%R已知值与计算值的差fprintf('用Simpson公式计算的结果Sn=')disp(Sn)fprintf('等分数n=')disp(n)fprintf('已知值与计算值的误差R=')disp(R)运行结果为:E=2.7284e-08用三inpson公式l+苴的结果Sn-7.3091等分劫In=24已知值与计算值的误差Hf2.7284e-08A»2.4 数值分析误差分析：在上述计算中，若采

10、用复化梯形公式，则可以知，其与精确值的误差为2.8300Xe-8,等分数为n=7019;而采用复化Simpson公式知与精确值的误差为2.7284X18,等分数为n=24o故与复化梯形公式相比，复化Simpson公式误差相对较小。?收敛性分析：右妒N?=0?(?)=心？?(？?)？?？，复化Simpson公式的余?>项是：一？-?一4一(4)?=-询？?()(?)，？?？C?，？?4一、一.一？一.可以看出误差是h阶，实际上若f(x)ec(a,b),?1mo?=%?(?)?,因此复化Simpson公式是收敛的。稳定性分析：由于求积公式中A>0(i=0,1,.，n)则求积公式是稳定的

11、。3矩阵的LU分解3.1问题背景矩阵的LU分解主要用来求解线性方程组或者计算行列式。在使用初等行变换法求12-1解线性方程组的过程中，系数矩阵的变化情况如下：A=310经过日2(-3)、E3(1)、-1-1-212-1&(1/5)可得到0-53000-12/5由上可知：&(1/5)Ei3(1)Ei2(-3)A=U其中U就是上面矩阵A经过行变换后的上三角矩阵，Eij表示将i行元素与j行元素互换的初等矩阵；Eij(k)表示将i行元素的k倍加到j行上。因此：A=Ei2（3）Ei3（-1）E23（-1/5）12-11A= 310 = 3-1 - 1 - 2 - 1如果方阵A可以分解成单

12、位下三角矩阵 A的LU分解或三角分解。00 12- 110 0 - 53 =LU-1/5 100- 12/5L与上三角矩阵U的乘积，则式A=LU称为3.2 数学模型3.2.1 理论基础矩阵的LU分解在求解线性方程组时将十分简便。如对线性方程组Ax=b,设人=15其LU分解。我们先求解方程组Ly=b0由于L是下三角矩阵，则解向量y可以通过依次求出其分量yby2,yn而求出，再求解方程组Ux=y。解向量x可以通过该方程组依次求出分量Xn,Xn.1,X2,X1而快速得出。于是由两个方程组Ux=y,Ly=b的求解而给出LUx=Ly=b=A的解。若矩阵A非奇异，则A能分解为LU的充分必要条件是A的顺序主

13、子行列式不为00?1= ?11 " ?2=?11?21?11?1?3=?333?1?则存在惟一的主对角线上元素全为1的下三角阵L与惟一的上三角阵U,使得A=LU3.2.2 实例102030将矩阵204580进行LU分解。30801713.3 计算方法程序：clearallclcA=input('请输入一个方阵);输入一个n阶方阵n,n=size(A);L=zeros(n,n);U=zeros(n,n);fori=1:n%将L的主对角线元素赋值1L(i,i)=1;endforj=1:n晒矩阵U勺第一行元素U户A(1,j)；endfork=2:n晒矩阵L的第一列元素L(k,1)=

14、A(k,1)/U(1,1);endfori=2:n晒L、U矩阵元素forj=i:ns=0;fort=1:i-1s=s+L(i,t)*U(t,j);endU(i,j)=A(i,j)-s；endfork=i+1:nr=0;fort=1:i-1r=r+L(k,t)*U(t,i);endL(k,i)=(A(k,i)-r)/U(i,i);endend喻出矢!阵L、ULU输入一个方阵，输出结果如下:畲令行葡匚请输入一个方造10 20 30.20 45 80r3d S3 171L-二inoo20S03C-2(3.4 小组元的误差例如线性方程组Ax = b ,其中:如果系数矩阵被扰动成？=A = 110- 2

15、01；,b=0，可得理论解 x=-11 ；,可手算知：？= 10120°?=110 20120o01 - 10- 20若上述过程在计算机中进行，由浮点数运算规则可知, 则计算机中产生的矩阵为：两数相加时，大数吃掉小数,? = ? , ?=?，? = 10020-10- 20这时会发现？ ？' = 10-20L 11,且据？x=b解出的理论解x'=，,明显不再等于前面的理论解。这说明LU分解是稳定的，但是将LU分解用到解线性方程组上是不稳定的。究其原因，是因为？?中的第一个主元10-20太小，导致第二个主元中的1与值10-20相差悬殊，出现大数吃小数。为了避免上述危害，

16、引入一种选主元手段，即在消去的过程中，通过适当的选主元，避免放大数据误差。常用的选主元技术就是列选主元法（除此之外还有全选主元法、对角选主元法和随机选主元法等）：对m<n阶矩阵A,在确定第k个主元？?（？?？>k）时，先从该列自主元位置(?)(?)(kjk)至列尾的所有兀素中选择绝对俏最大的兀素与？?也换妹后将？?)1?IX,JrvAtsHJ"InJyI十Lu/、vHJy?J?O,I、，-I?+1,?,?化为零。继续这个过程，直至将矩阵A化为行阶梯形。4二分法求方程的根4.1 问题背景在科学研究与工程计算中，常遇到方程(组)求根问题。若干个世纪以来，工程师和数学家花了大量

17、时用于探索求解方程(组)，研究各种各样的方程求解方法。对于方程f(x)=0,当f(x)为线性函数时，称f(x)=0为线性方程；当f(x)为非线性函数时，称式f(x)=0为非线性方程。对于线性方程(组)的求解，理论与数值求法的成果丰富；对于非线性方程的求解，由于f(x)的多样性，尚无一般的解析解法。当f(x)为非线性函数时，若f(x)=0无解析解，但如果对任意的精度要求，设计迭代方程，数值计算出方程的近似解，则可以认为求根的计算问题已经解决，至少能够满足实际要求。4.2 数学模型使用二分法求方程x”+x-1=0在0,1内的近似根(误差<10A-5)。二分法：二分法是最简单的求根方法，它是利

18、用连续函数的零点定理，将含根区间逐次减半缩小，取区间的中点构造收敛点列xk来逼近根x。4.3 计算方法二分法程序代码：functiony=erfen1(m,n,er)symsxxka=m;b=n;k=0;ff=xA3+x-1;whileb-a>erxk=(a+b)/2;fx=subs(ff,x,xk);fa=subs(ff,x,a);k=k+1;iffx=0y(k)=xk;break;elseiffa*fx<0b=xk;elsea=xk;endy(k)=xk;endplot(y,'.-');gridon在命令窗口下执行:>>a-b-ezfenlC0T1,

19、1电-5):vpajatjS)实验结果如下：可以得到迭代区间中点数列分布及图像，数值如下：ans=0.5,0.75,0.625,0.6875,0.65625,0.671875,0.6796875。0.68359375,0.68164062,0.68261719,0.68212891,0.68237305,0.68225098,0.68231201,0.68234253,0.68232727,0.6823349依据题目要求的精度，则需彳十七次，由实验数据知x=0.6823349即为所求的根。4.4 二分法的收敛性二分法算法简单，容易理解，且总是收敛的；但是二分法收敛速度太慢，浪费时间，并且不能求

20、复根跟偶数重根。5雅可比迭代求解方程组5.1 问题背景在自然科学和工程技术中有很多问题的解决常常归结为解线性方程组，例如电学中的网络问题，化学中的配平方程式模型问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，非线性方程组求解问题，用差分法或者有限元法解常微分方程、偏微分方程边值问题等都导致求解线性方程组。在实践中，通常采用数值方法来讨论线性方程组Ax=b的解，其中迭代法是一种重要方法。5.2 数学模型5.2.1 理论基础雅可比迭代法：首先将方程组中的系数矩阵A分解成三部分，即：A=L+D+U其中D为对角阵，L为下三角矩阵，U为上三角矩阵。之后确定迭代格式，XA(

21、k+1)=B*XA(k)+f,(这里A表示的是上标，括号内数字即迭代次数)，其中B称为迭代矩阵，雅克比迭代法中一般记为J(k=0,1,)，再选取初始迭代向量XA(0),开始逐次迭代。设Ax=b,其中A=D+L+UZ非奇异矩阵，且对角阵D也非奇异，则当迭代矩阵J的谱半径p(J)<1时，雅克比迭代法收敛对于Ax=b, 0 ?91A= .7?1?2?3+0?11?220?12?13？1?0?22？2?+0？?3? r?, 0 =L+D+U因为？?？ w0(Jacob 假设)所以D可逆。故有：(L+D+U)x=bDx=-(L+U)x+bx=-D-1 (L+U)+D-1b5.2.2 实例用雅可比迭

22、代法解方程组:430?12435- 1 ?2= 30 0 - 14?3- 245.3计算方法雅可比迭代法程序:functionx,k,index=Jacobi(A,b,ep,it_max)ifnargin<4it_max=100;endifnargin<3ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;whilek<=it_maxfori=1:nifabs(A(i,i)<1e-10index=0;return;endy(i)=(b(i)-A(i,1:n)*x(1:n)+A(i,i)*x(i)/A(

23、i,i);endifnorm(y-x,inf)<epbreak;endk=k+1;x=y;end在命令窗口输入的命令和结果如下图：>>A二30;35-I;0-14;b=243。-24J;ep=le_5;it_nai=l00;kjinlfixI=Jacobi(A3bjit_JiLSM1x-4.20002.4000-b.400039indei=15.4收敛性分析由上面运行的结果可知方程组的解为4.2000,2.4000,-5.4000,迭代次数为39,由index=1表示计算成功。雅克比迭代法的优点明显，计算公式简单，每迭代一次只需计算一次矩阵和向量的乘法，且计算过程中原始矩阵A

24、始终不变，比较容易并行计算。然而这种迭代方式收敛速度较慢，而且占据的存储空间较大，所以工程中一般不直接用雅克比迭代法，而用其改进方法。与逐次超松弛迭代法相比，雅可比迭代法收敛速度相对较慢，逐次超松弛迭代法收敛速度较快。由逐次超松弛迭代法求出的方程组的数值解与该方程组的精确解十分接近。并且离散化后线性方程组的逐次超松弛迭代法的精确性较高，故相对于雅可比迭代法，逐次超松弛迭代法更加广泛地应用于实际，可以用逐次超松弛迭代法求解高阶稀疏线性方程组。6Romberg求积法6.1 问题背景复化求积方法对于提高精度是行之有效的方法，但复化公式的一个主要缺点在于要事先估计出部长。若步长过大，则精度难于保证；若

25、步长过小，则计算量又不会太大。而用复化公式的截断误差来估计步长，其结果是步长往往过小，而且f(x)和f(4)(x)在区间a刈上的上界M的估计是较为困难的。在实际计算中通常采用变步长的方法，即把步长逐次分半(也就是把步长二等分)，直到达到某种精度为止，这种方法就是Romberg积分法的思想。6.2 数学模型：6.2.1 理论基础记：TN0)TN,将区间N等分的梯形值。TN1)SN,将区间N等分的Simpson；T2)CN,将区间N等分的CotesoTN3)RN,将区间N等分的Romberg由其可构造一个序列TN,次序列称为Romberg字列，并满足如下递推关系:f(a (2k 1),2N k12

26、NT1詈f(a)f(b),T2(N)2TN(0)TNkk (k 1) (k 1)4 12N 1 N4k 1,k 1,2,L以上递推公式就是Romberg积分递推公式。6.2.2 实例.一一.1,2一4八利用Romberg积分法的思想，用龙贝格公式数值积分求函数乙,2?的积分6.3 计算方法龙贝格公式积分程序：functionI,step=romberg(f,a,b,eps)%romberg.m为用龙贝格求积分%f为被积函数%a.b为积分区间的上下限%eps为积分结果精度%I为积分结果；step为积分的子区间数if(nargin=3)eps=1.0e-4;end;M=1;tol=10;k=0;T

27、=zeros(1,1);h=b-a;T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b);whiletol>epsk=k+1;h=h/2;Q=0;fori=1:Mx=a+h*(2*i-1);Q=Q+subs(sym(f),findsym(sym(f),x);endT(k+1,1)=T(k,1)/2+h*Q;M=2*M;forj=1:kT(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j)/(4Aj-1);endtol=abs(T(k+1,j+1)-T(k,j);endI=T(k

28、+1,k+1);step=k;输入：»cLearall;a=-L2：b=l<2：eps=le-fi：I>step=rcTnberg1Jx"ta?bjeps)可以得到结果：Q,995323000000000step-6.4 误差分析Romberg求积公式是具有八阶精度的算法，收敛且稳定，比梯形公式，Simpson公式以及Cotes公式收敛的快。龙贝格公式的余项为：?ccc?,?(?)=/?(?)?-?,?=-+1)(?+2?)2?(?-?)一+?"+(?)?2Romberg积分公式高速有效，易于编程，适合于计算机计算。但他有一个主要的缺点是，每当把积分区

29、间对分后，就要对被积函数？?(？?)计算它在新分点处的值，而这些函数值的个数是成倍增加的。7曷法7.1 问题背景工程及物理中的许多实际问题需要计算矩阵的特征值及特征向量，对于给定的n阶实矩阵A,当n很小时用传统的方法是可以的，但当n稍大时，计算工作量将以惊人的速度增大，并且由于计算存在误差，用方程(入I-A)x=0求解十分困难。在实际应用中，有的时候不一定需要求出矩阵的全部特征值和特征向量，例如只需要求出按模最大的或最小的特征值。求这类特征值的方法，通常采用迭代法，其中两种是幕法和反幕法。7.2 数学模型7.2.1 理论基础设An有n个线性相关的特征向量v1,v2,，vn,对应的特征值i,2,

30、n,满足|i|>|2|n|,因为v1,v2,，vn为Cn的一组基，所以任给x0,(0)XaiVii1,所以有:n/k(0)«k/、AxA(aiVi)i1nkaiiVki1nkaiAVii1nkik1aV1()aiVii21若a10,则因1知，当k充分大时Ak)x(0)1ka1V1=CV1属人1的特征向量另一方面,记max(x)=Xi,其中|Xi|二|x|，则当k充分大时,max(Akx(0)max(1、心)1kmax(a1V1)/Ak1、.(叭zk-kzT1max(Ax)max(1aw)1max(aM)若d=0,则因舍入误差的影响，会有某次迭代向量在v1方向上的分量不为0,迭代

31、下去可求得入1及对应特征向量的近似值。7.2.2 实例用幕法计算矩阵3-10?=-132023绝对值最大的特征值及相应的特征向量，取精度要求为c=10-407.3 计算方法幕法Matlab程序：functionm,u,index=pow(A,ep,it_max)ifnargin<3it_max=100;endifnargin<2ep=1e-5;endn=length(A);u=ones(n,1);index=0;k=0;m1=0;whilek<=it_maxv=A*u;,i=max(abs(v);m=v(i);u=v/m;ifabs(m-m1)<epindex=1;br

32、eak;endm1=m;k=k+1;end在命令窗口输入和输出如下图所示：>、汽=3-I0-132:0231：mjufind«xj-po¥(At1b1)nt机2301-0,4136L0C0Q0,9112index"7.4误差分析幕法程序可以用来计算矩阵绝对值最大的特征值及相应的特征向量。幕法的缺点是开始的时候并不知道矩阵是否有单一的主特征值，也不知道如何选择Vo以保证它关于矩阵特征向量的表达中包含一个与主特征值相关的非零特征向量。采用幕法和反幕法，求矩阵的最大和最小特征值，从原理上看，这两种方法都是迭代法，因此迭代初始向量的选择对计算结果会产生一定影响，主要表现在收敛速度上。同时，原点位移m的选取也影响收敛的速度，但原点位移mo的