数值分析实验报告

1、数值分析实验报告姓名：张献鹏学号：173511038专业：冶金工程班级：重冶二班1.1问题背景11.2数学模型11.3 计算方法11.4数值分析2复化辛普森求积公式22.1问题背景22.2数学模型32.3 计算方法32.4数值分析5矩阵的LU分解53.1 问题背景53.2数学模型63.2.1理论基础. 63.2.2 实例63.3计算方法63.4 小组兀的误差 8一分法求方程的根 94.1 问题背景94.2数学模型94.3 计算方法94.4二分法的收敛性. 10雅可比迭代求解方程组. 105.1 问题背景105.2数学模型115.2.1理论基础. 115.2.2 实例11拉格朗日插值.12345

2、.3 计算方法 115.4 收敛性分析 136 Romberg求积法 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 问题背景1 f (X) 2£ / 对于函数1 X

3、, 5 x 5求拉格朗日插值。n 10，把f(x)和插值多项式的曲线画在同一张图上进行比较，观察数值积分中的Lagrange插值。1.2 数学模型取等距差值节点?=-5+10?/n, ?=0, 1,.,n,构造n次lagrange插值多项式：? 1?+1(?)?=乙22/?=01 + ?(?- ?)?"1(?)当n=10时，十次插值多项式L1o(x)以及函数f (x)的图像可以由Matlab画出1.3 计算方法f.m :fun cti onf= f( x )f=1./(1+x.A2);endLagra nge.mfun cti ony=Lagra nge(x0,y0,x);n=len

4、 gth(x0);m=le ngth(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1: nif j=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;yO=Lagra nge(x,y,xO);y1=1./(1+x0.A2);plot(xO,yO, 'b')hold onplot(x0,y1, 'r')运行这个文件可以得到如下拉格朗日图形曲线:1.4 数值分析

5、Le( x)的断误差Ro( x)= f (x)- Le(x)在区间-5 ， 5的两端非常大。例如， L1o(4.8)=1.8O438，而f (4.8)=0.04160。这种现象称之为龙格现象。不管n取多大，Runge 现象依然存在。因此，对函数作插值多项式时，必须小心处理，不能认为差值节点取得越多，差值 余项就越小。此外，当节点增多时，舍入误差的影响不能低估。为了克服高次插值的不 足，应采用分段低次插值。2复化辛普森求积公式2.1 问题背景2用复化Simpson公式计算定积分？?2 =彳？?？?的近似值，要求误差限E=1/2 X 10-7，利用其余项对算法做出步长的事前估计；并将计算结果与精确

6、值进行比较2.2 数学模型将积分区间a,b分为n等分，h=(b-a)/ n, Xk=a+ k h , k=0, 1,n。在每个子区 间Xk, Xk+i上用Simpson公式可得：?-1/ ?(x)dx=刀?(?)?=0?- 1?"百 刀?(?) + 4?(?+1)+ ?(?+1) ?= 02其中 Xk+1/2=Xk+1/2h。?_ 1?Sn ? = -6 刀?(?) + 4?(?+6 ?= 0?2)+?(?7?+1) = _?(?)?- 1?- 1+ 4 刀?(? 1)+2 刀?(?) + ?(?)+ 2?= 0 2 ?= 1此式即为复化Simpson公式。设f (x) C4a,b，

7、由Simpson公式的误差有5?= ?- ?= E?=1- 90(?) ?(4)(?) , ? ?,?+1。则复化Simpson公式的余项为：?- ?4?=-镒？ ？?( "(？?), ? ?, ?复化Simpson公式四阶收敛2.3 计算方法程序1(求f(x)的n阶导数)：syms xf=x*exp(x)%t义函数 f (x)n=input('输入所求导数阶数:')f2=diff(f,x,n)%求 f(x)的 n阶导数程序2:clcclearsyms x%定义自变量xf=inline( 'x*exp(x)' , 'x' )%定义函数f

8、(x)=x*exp(x)，换函数时只需换该函数表达式即可f2=i nlin e('(4*exp(x) + x*exp(x)','x' )%定义f(x)的四阶导数，输入程序1里求出的f2即可f3= '-(4*exp(x) + x*exp(x)'表达式带引号，故取负号，一边求最大值e=5*10A(-8)a=1%积分下限b=2%积分上限x1= fmi nbn d(f3,1,2)大值点对应的x值for n=2:1000000%因 fminbnd ()函数求的是表达式的最小值，且要求%精度要求值%求负的四阶导数的最小值点，也就是求四阶导数的最%求等分数nR

9、n=-(b-a)/180*(b-a)/(2* n)F4*f2(x1) if abs(Rn)<ebreakendendh=(b-a)/nSn 1=0Sn 2=0for k=0:n-1xk=a+k*hxk仁xk+h/2Sn仁Sn 1+f(xk1)Sn 2=S n2+f(xk)end%计算余项%用余项进行判断%符合要求时结束%求 h%求两组连加和Sn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a)+f(b)御 Sn2多加了 k=0时的值，故减去 f(a)z=exp(2)R=S n-z已知值与计算值的差fprintf('用Simpson公式计算的结果Sn=')disp(S

10、 n)fprintf('等分数 n=')disp (n)fprintf('已知值与计算值的误差R=')disp(R)运行结果为:2. 723J=e-Q8用匚impEwii式l+宜的谿果二 7. 3S91筹分埶炉 24已知值与讣負值的漠堇R= 2. 7281S-08A»2.4 数值分析误差分析：在上述计算中，若采用复化梯形公式，则可以知，其与精确值的误差为 2.8300 Xe-8，等分数为n=7019;而采用复化Simpson公式知与精确值的误差为2.7284 Xe" 8,等分数为n=24。故与复化梯形公式相比，复化 Simpson公式误差相对

11、较小。? ?收敛性分析：若 lim 刀?=o?(?) = %?(?)?，复化 Simpson公式的余项是：?_ ?4(4)?=-歸？ ?()(?)，? ?，?4?可以看出误差是h阶，实际上若f(x) GC(a,b)，?io?= £?(?)?，因此 复化Simpson公式是收敛的。稳定性分析：由于求积公式中 A>0(i=0, 1，.，n)则求积公式是稳定的。3矩阵的LU分解3.1 问题背景矩阵的LU分解主要用来求解线性方程组或者计算行列式。在使用初等行变换法求-10 经过 E12(_3)、E13(1)、-21 2解线性方程组的过程中，系数矩阵的变化情况如下：A= 31-1 - 1

12、1 2 - 1巳(1/5)可得到0 - 5300- 12/5由上可知：E23(1/5)日3(1) E12(_3)A=U其中U就是上面矩阵A经过行变换后的上三角矩阵，Eij表示将i行元素与j行元 素互换的初等矩阵；Eij (k)表示将i行元素的k倍加到j行上。因此：A壬12(3)巳3(_1) E23(_1/5)12-11A= 310 = 3-1 - 1 - 2 - 10 0 1 2 - 110 0 - 53 =LU如果方阵A可以分解成单位下三角矩阵L与上三角矩阵U的乘积，则式A=LU称为-1/5 100- 12/5A的LU分解或三角分解。3.2 数学模型3.2.1 理论基础矩阵的LU分解在求解线

13、性方程组时将十分简便。 如对线性方程组Ax=b，设A=LU是 其LU分解。我们先求解方程组Ly=b。由于L是下三角矩阵，则解向量y可以通过依次 求出其分量y1, y2, , yn而求出，再求解方程组Ux=y。解向量x可以通过该方程组 依次求出分量Xn，Xn-1， ，X2, X1而快速得出。于是由两个方程组 Ux=y，Ly=b的求解 而给出LUx=Ly=b=Ax的解。若矩阵A非奇异，则A能分解为LU的充分必要条件是A的顺序主子行列式不为0。?11?，?3= ?1 ?1= ?11 工 0,?2= ?11?21?1?则存在惟一的主对角线上元素全为1的下三角阵L与惟一的上三角阵U,使得A=LU3.2.

14、2 实例10 20将矩阵20 4530 803080进行LU分解1713.3 计算方法程序：clear allclcA=input('请输入一个方阵);输入一个n阶方阵n,n=size(A);L=zeros(n,n);U=zeros(n,n);for i=1:n % 将L的主对角线元素赋值1L(i,i)=1;endfor j=1:n%求矩阵U勺第一行元素U(1,j)=A(1,j);endfor k=2:n%求矩阵L的第一列元素L(k,1)=A(k,1)/U(1,1);endfor i=2:n%求 L、U矩阵元素for j=i:ns=0;for t=1:i-1s=s+L(i,t)*U(t

15、,j);endU(i,j)=A(i,j)-s;endfor k=i+1:nr=0;for t=1:i-1e+L(k,t)*U(t,i);endL(k,i)=(A(k,i)-r)/U(i,i);endend%输出矩阵L、ULU输入一个方阵，输出结果如下:缽口询输入一个方畸10 20 30,20 45 S0r3Q 3? H叮L =00210341匸=1020300S200013.4 小组元的误差例如线性方程组Ax = b，其中:如果系数矩阵被扰动成?=A = 010- 201；，b=Q，可得理论解X=-11 ；，可手算知：?= 蔦。0?=110- 201 0 1 - 10- 2°。若上述

16、过程在计算机中进行，由浮点数运算规则可知,则计算机中产生的矩阵为：两数相加时，大数吃掉小数,? = ?，?=?，? = 10020-10-2o?化为零。继续这个过程，直至将矩阵A化为行阶梯形。，d n- 20这时会发现? ? = 101；，且据? ? x=b解出的理论解x'=【0】，明显不再等于前面的理论解。这说明LU分解是稳定的，但是将LU分解用到解线性方程组上是不稳定的。究其原 因，是因为？?中的第一个主元10-20太小，导致第二个主元中的1与值10-20相差悬殊，出 现大数吃小数。为了避免上述危害，引入一种选主元手段，即在消去的过程中，通过适当的选主元， 避免放大数据误差。常用的

17、选主元技术就是列选主元法 （除此之外还有全选主元法、对 角选主元法和随机选主元法等）：对mK n阶矩阵A,在确定第k个主元？?（？?？ >k）时，先从该列自主元位置（k,jk）至列尾的所有兀素中选择绝对值最大的兀素，与?交换，然后将？?+',?，4二分法求方程的根4.1 问题背景在科学研究与工程计算中，常遇到方程(组)求根问题。若干个世纪以来，工程师和数学家花了大量时用于探索求解方程(组)，研究各种各样的方程求解方法。对于方程f(x)=O，当f(x)为线性函数时，称f(x)=O为线性方程；当f(x)为非线性函数时，称式 f(x)=O为非线性方程。对于线性方程(组)的求解，理论与数

18、值求法的成果丰富；对于 非线性方程的求解，由于f (x)的多样性，尚无一般的解析解法。当 f (x)为非线性函数 时，若f(x)=O无解析解，但如果对任意的精度要求，设计迭代方程，数值计算出方程 的近似解，则可以认为求根的计算问题已经解决，至少能够满足实际要求。4.2 数学模型使用二分法求方程xA3+x-1=0在0,1内的近似根(误差<10A-5 )。二分法：二分法是最简单的求根方法，它是利用连续函数的零点定理，将含根区间 逐次减半缩小，取区间的中点构造收敛点列Xk来逼近根x。4.3 计算方法二分法程序代码：fun cti on y=erfe n1(m, n,er) syms x xk

19、a=m;b=n ;k=0;ff=xA3+x-1;while b-a>erxk=(a+b)/2;fx=subs(ff,x,xk);fa=subs(ff,x,a);k=k+1;if fx=Oy(k)=xk;break;elseif fa*fx<0b=xk;elsea=xk;endy(k)=xk;endplot(y, '.-')；grid on在命令窗口下执行：>> abezfenl lf le-5):实验结果如下：可以得到迭代区间中点数列分布及图像，数值如下：ans =0.5，0.75，0.625，0.6875，0.65625，0.671875，0.6796

20、875。0.68359375， 0.68164062，0.68261719，0.68212891，0.68237305，0.68225098，0.68231201， 0.68234253，0.68232727，0.6823349T Figure 1-2 I寸曰*立粹田郝削D 却者聞 JS人(D 工也I：空直财 *(»)群城11J亠Ljdda丸耳背舉恵回目 依据题目要求的精度，则需做十七次，由实验数据知 x=0.6823349即为所求的 根。4.4 二分法的收敛性二分法算法简单，容易理解，且总是收敛的；但是二分法收敛速度太慢，浪费时间， 并且不能求复根跟偶数重根。5雅可比迭代求解方程组

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

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

23、= 304?3- 245.3 计算方法雅可比迭代法程序:fun cti onx,k,i ndex=Jacobi(A,b,ep,it_max)if nargin<4it_max=100;endif nargin<3ep=1e-5;endn=le ngth(A);k=O;x=zeros( n,1);y=zeros( n,1);i ndex=1;while k<=it_maxfor i=1:nif abs(A(i,i)<1e-10in dex=0;return ;endy(i)=(b(i)-A(i,1: n)*x(1: n)+A(i,i)*x(i)/A(i,i);endif

24、norm(y-x,inf)<epbreak;endk=k+1;x=y;end在命令窗口输入的命令和结果如下图：>>3 0;3 5 -I;0 -1 4 ;b=24 刃-24' :ep=le5;it_nai=l00;kj inlfi x I = J ac ob 1 (A3 bj 电p, i七Jtlqm 1x -4. 20002.4000-b.400039index =15.4 收敛性分析由上面运行的结果可知方程组的解为4.2000 , 2.4000 , -5.4000,迭代次数为39, 由index=1表示计算成功。雅克比迭代法的优点明显，计算公式简单，每迭代一次只需计算

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

26、主要缺点在于要 事先估计出部长。若步长过大，则精度难于保证；若步长过小，则计算量又不会太大。 而用复化公式的截断误差来估计步长， 其结果是步长往往过小，而且f（X）和f"）（x）在 区间a,b上的上界M的估计是较为困难的。在实际计算中通常采用变步长的方法，即 把步长逐次分半（也就是把步长二等分），直到达到某种精度为止，这种方法就是Romberg 积分法的思想。6.2 数学模型:6.2.1 理论基础记：Tf Tn，将区间N等分的梯形值。tN）Sn，将区间N等分的Simpson； t/ cn，将区间N等分的Cotes。tN® Rn，将区间N等分的Romberg 由其可构造一个序

27、列Tn（），次序列称为Romberg序列，并满足如下递推关系:T1(0)f(a)f(b),T2(N)b a N2N k 1f (a (2 k 1)b a2N),4kT(k 1) T(k 1)TnJ2；n ,k 1,2 丄41以上递推公式就是6.2.2 实例Romberg积分递推公式。124利用Romberg积分法的思想，用龙贝格公式数值积分求函数/仁2? ?的积分6.3 计算方法龙贝格公式积分程序：fun cti on l,step=romberg(f,a,b,eps)%romberg.m为用龙贝格求积分%f为被积函数%a.b为积分区间的上下限%eps为积分结果精度%I为积分结果；step为积

28、分的子区间数if (nargin=3)eps=1.0e-4;end;M=1;tol=10;k=0;T=zeros(1,1);h=b-a;T(1,1)=(h/2)*(subs(sym(f),fi ndsym(sym(f),a)+subs(sym(f),fi ndsym(sym(f),b);while tol>epsk=k+1;h=h/2;Q=0;for i=1:Mx=a+h*(2*i-1);Q=Q+subs(sym(f),fi ndsym(sym(f),x);endT(k+1,1)=T(k,1)/2+h*Q;M=2*M;for j=1:kT(k+1,j+1)=T(k+1,j)+(T (k+

29、1,j)-T(k,j)/(4Aj-1);endtol=abs(T(k+1,j+1)-T(k,j);endI=T(k+1,k+1);step=k;输入:» clear all;a=-l 2：b=l 2：eps=lt-6:I> step-rcniberg ；，:<”/ t beps)可以得到结果:0. 995326000000000st ep -6.4 误差分析Romberg求积公式是具有八阶精度的算法，收敛且稳定，比梯形公式，Simpson公式以及Cotes公式收敛的快。龙贝格公式的余项为：?(?) = ?(?)?- ?=-?2?+ 22(?+1)( ?+2?)2?!(?-

30、 ?)2?+3?(2?+2)(?)Romberg积分公式高速有效，易于编程，适合于计算机计算。但他有一个主要的缺 点是，每当把积分区间对分后，就要对被积函数？?（？?）计算它在新分点处的值，而这些函数值的个数是成倍增加的。7 幂法7.1 问题背景工程及物理中的许多实际问题需要计算矩阵的特征值及特征向量，对于给定的n阶实矩阵A,当n很小时用传统的方法是可以的，但当 n稍大时，计算工作量将以惊人的速度增大，并且由于计算存在误差，用方程（入l-A）x=0求解十分困难。在实际应用中， 有的时候不一定需要求出矩阵的全部特征值和特征向量， 例如只需要求出按模最大的或 最小的特征值。求这类特征值的方法，通常

31、采用迭代法，其中两种是幕法和反幕法。7.2 数学模型7.2.1 理论基础满足| i| > | 2|设An有n个线性相关的特征向量vi, v2,，vn，对应的特征值i, 2,，n, n| ,因为vi , v2,，vn为Cn的一组基，所以任给x(0) 0,n(0)xaivii i ，所以有:若ai0,贝U因Akx(0)nk ,、A (aM)i inkai i vki in3i Ai ii aiVi当k充分大时A/kVin k()aMi 2 iikaiVi = cvi属入i的特征向量另一方面，记max(x)=Xi，其中| Xi| = |x|，则当k充分大时,max(Ak x(0)max(Ak

32、ix(0)kkmax( i aivi)i max(aivi)FTmax( i aivi)iimax(aiVi)会有某次迭代向量在 下去可求得入i及对应特征向量的近似值。若ai=0,则因舍入误差的影响,vi方向上的分量不为0，迭代7.2.2 实例用幕法计算矩阵0233?= - i0绝对值最大的特征值及相应的特征向量，取精度要求为£ =i0-4。7.3 计算方法幕法Matlab程序：fun cti onm,u,i ndex=pow(A,ep,it_max)if nargin<3 it_max=i00; end if nargin<2 ep=ie-5; end n=len gt

33、h(A);u=ones(n ,i);in dex=0;k=0;mi=0;while k<=it_maxv=A*u;,i=max(abs(v); m=v(i);u=v/m;if abs(m-m1)<epindex=1;break ;endm1=m;k=k+1;end在命令窗口输入和输出如下图所示：» A=f3 -1 0;-1 3 2;0 2 3:口index(A, 1m th 2301-0,4136LOCOQ0,9112index 7.4 误差分析幕法程序可以用来计算矩阵绝对值最大的特征值及相应的特征向量。 幕法的缺点是 开始的时候并不知道矩阵是否有单一的主特征值，也不知道如何选择V。以保证它关于矩 阵特征向量的表达中包含一个与主特征值相关的非零特征向量。采用幕法和反幕法，求矩阵的最大和最小特征值，从原理上看，这两种方法都是迭 代法，因此迭代初始向量的选择对计算结果会产生一定影响，主要表现在收敛速度上。同时，原点位移m的选取也影响收敛的速度，但原点位移m的适当选取