数值分析报告实验报告材料

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 问题

2、背景 101.2 数学模型 111.2.1 理论基础 111.2.2 实例 111.3 计算方法 111.4 收敛性分析 126 Romber球积法 136.1 问题背景 136.2 数学模型： 136.2.1 理论基础 136.2.2 实例 146.3 计算方法 146.4 误差分析 157 幕法 157.1 问题背景 157.2 数学模型 157.2.1 理论基础 157.2.2 实例 167.3 计算方法 167.4 误差分析 178 改进欧拉法 178.1 问题背景 178.2 数学模型 178.2.1 理论基础 178.2.2 实例 188.3 数学模型 188.4 误差分析 19文

3、案大全1拉格朗日插值1.1 问题背景1f (x) 2对于函数 1 x , 5 X 5求拉格朗日插值。n 10,把f(x)和插值多项式的曲线画在同一张图上进行比较，观察数值积分中的Lagrange插值。1.2 数学模型取等距差值节点？?=-5+10?n, ?=0, 1, . , n,构造n次lagrange插值多项式:?= Z2?=01?“(？)1 + ?2?(? ?街?+1(?当n=10时，十次插值多项式L0(x)以及函数f(x)的图像可以由Matlab画出1.3 计算方法f.m :function f= f( x ) f=1./(1+x.A2);endLagrange.mfunction y

4、=Lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:n if j=k p=p*(z-x0(j)/(x0(k)-x0Q);endends=p*y0(k)+s;endy(i)=s；End拉格朗日插值的曲线:x=-5:1:5;y=1./(1+x.A2);xO=-5:0.001:5;y0=Lagrange(x,y,x0);y1=1./(1+x0A2);plot(x0,y0, hold on'b')plot(x0,y1, 'r')运行这个文件可以得

5、到如下拉格朗日图形曲线:1.4 数值分析Lio(x)的断误差Rio(x)= f(x)-Lio(x)在区间卜5 , 5的两端非常大。例如， Lio(4.8)=1.8O438，而f (4.8)=0.04160。这种现象称之为龙格现象。不管n取多大，Runge 现象依然存在。因此，对函数作插值多项式时，必须小心处理，不能认为差值节点取得越多，差值 余项就越小。此外，当节点增多时，舍入误差的影响不能低估。为了克服高次插值的不 足，应采用分段低次插值。2复化辛普森求积公式2.1 问题背景用复化Simpson公式计算定积分？ = /?2?我近似值，要求误差限己=1/2X10- 7,利用其余项对算法做出步长

6、的事前估计；并将计算结果与精确值进行比较。2.2 数学模型将积分区间a,b分为n等分，h=(b-a)/ n, Xk=a+ k h , k=0, 1,门。在每个子区 问Xk, Xk+i上用Simpson公式可得：?-i?-i?(?+i)? ?x)dx = E ?(?)? E ?(?) + 4?7?j) +?=06 ?=02其中 Xk+1/2=Xk+1/2h。?-1?-1?-1Sn? = yE ? + 4?(?,?+1) + ?+1) =8? + 4 汇？7?g) + 2 汇？?钩 + ?(?) ?=0?=0?=1此式即为复化Simpson公式。设f(x) C4a,b,由Simpson公式的误差有

7、一 1? 5?= ? ?= U?=1o- 90(7) ?夕)(?为，?,?+1o则复化Simpson公式的余项为：一 一?-? 一 4 一?=- 288880-? ?夕)(?九? ? ?复化Simpson公式四阶收敛。2.3 计算方法程序1(求f(X)的n阶导数)：syms xf=X*eXp(X)%定义函数 f (x)n=input('输入所求导数阶数:')f2=diff(f,x,n)%t f(x)的 n阶导数程序2:clcclearsyms x%!义自变量xf=inline( 'x*exp(x)' , 'x' )%!义函数f(x)=x*exp(

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

9、*n)A4*f2(x1)府算余项endif abs(Rn)<e%用余项进行判断endh=(b-a)/nSn1=0Sn2=0break%符合要求时结束%R hfor k=0:n-1%求两组连加和xk=a+k*hxk1=xk+h/2Sn1=Sn1+f(xk1)Sn2=Sn2+f(xk)endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a)+f(b)%H Sn2力口了 k=0时的值，故减去f（a）z=exp(2)R=Sn-z%R已知值与计算值的差fprintf（'用Simpson公式计算的结果Sn='）disp(Sn)fprintf('等分数 n='

10、)disp(n)fprintf（'已知值与计算值的误差 R='）disp(R)运行结果为:E =2. 7284= e-0®用三inpwon公式l+算的结果£n=7. 3S91等分薮“24已知值与计算值的误差R= 2.7281e-08A»2.4数值分析误差分析：在上述计算中，若采用复化梯形公式，则可以知，其与精确值的误差为 2.8300 Xe8,等分数为n=7019;而采用复化Simpson公式知与精确值的误差为 2.7284 Xe-8,等分数为n=24o故与复化梯形公式相比，复化 Simpson公式误差相对较小。收敛性分析：若lim Zj?=d?钞

11、=/：?复化Simpson公式的余项是： ?° ,?-?4?=-痂？ ？夕)(?冽，？?？? ?可以看出误差是 h4阶，实际上若f(x) CC(a,b) , lim ?=4?(?)?？此复化 ?今Simpson公式是收敛的。稳定性分析：由于求积公式中 A>0(i=0, 1,.，n)则求积公式是稳定的。3矩阵的LU分解3.1 问题背景矩阵的LU分解主要用来求解线性方程组或者计算行列式。在使用初等行变换法求2-110经过E12（-3）、-1-2表示将i行元素与j行元1解线性方程组的过程中，系数矩阵的变化情况如下：A= 3-112-1日3(1)、E3(1/5)可得到0-53。00-1

12、2/5由上可知：E23(1/5) E13(1) E12(-3)A=U其中U就是上面矩阵A经过行变换后的上三角矩阵，Eij素互换的初等矩阵；Eij (k)表示将i行元素的k倍加到j行上 因此：A=E12(3) E13(-1) &(-1/5)12-110A = 310 = 31-1-1-2-1-1/50 12-10 0 -53 =LU1 00-12/5如果方阵A可以分解成单位下三角矩阵 L与上三角矩阵U的乘积，则式A=LU称为A的LU分解或三角分解。3.2 数学模型3.2.1 理论基础矩阵的LU分解在求解线性方程组时将十分简便。 如对线性方程组Ax=b,设人=15 其LU分解。我们先求解方

13、程组Ly=b0由于L是下三角矩阵，则解向量y可以通过依次 求出其分量yb y2, , yn而求出，再求解方程组Ux=y。解向量x可以通过该方程组 依次求出分量Xn, Xn-1, , X2, X1而快速得出。于是由两个方程组 Ux=y, Ly=b的求解 而给出LUx=Ly=b=A的解。?若矩阵A非奇异，则A能分解为LU的充分必要条件是A的顺序主子行列式不为00则存在惟一的主对角线上元素全为1的下三角阵L与惟一的上三角阵U,使得A=LU10将矩阵20302045803.2.2 实例3080 进彳亍LU分解。1713.3 计算方法程序：clear allclcA=input('请输入一个方阵

14、);输入一个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户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,j);endU(i,j)=A(i,j)-s；endfor k=i+1:nr=0;for t=1:i-1r=r+L(k,t)*U(t,i);endL(k,i)=(A(k

15、,i)-r)/U(i,i);endend喻出矢!阵L、ULU输入一个方阵，输出结果如下:畲令行葡匚请输入一个方晒10 20 30.20 45 80.3(1 盯 171L 二00210341L020300S200D13.4 小组元的误差例如线性方程组Ax = b ,其中：A= 0 1, b=1,可得理论解x=-1 。 1101如果系数矩阵被扰动成？= 10：01,可手算知：？N -20 0, ?=11101 10-20101 - 10-20 0若上述过程在计算机中进行，由浮点数运算规则可知，两数相加时，大数吃掉小数， 则计算机中产生的矩阵为：?77 = ?7 ?"? = r10 201

16、1'一 .'0 0-10 -20 这时会发现？ = 101200,且据？ ？ x=b解出的理论解x =0,明显不再等于前面的理论解。这说明LU分解是稳定的，但是将LU分解用到解线性方程组上是不稳定的。究其原 因，是因为？中的第一个主元10-20太小，导致第二个主元中的1与值10-20相差悬殊，出 现大数吃小数。为了避免上述危害，引入一种选主元手段，即在消去的过程中，通过适当的选主元， 避免放大数据误差。常用的选主元技术就是列选主元法 （除此之外还有全选主元法、对 角选主元法和随机选主元法等）：对“< n阶矩阵A,在确定第k个主元？?）?照>k）时，先从该列自主元位置

17、（k,jk）至列尾的所有元素中选择绝对值最大的元素，与 蜜?衮换，然后将？岸？?"?？?隽为 零。继续这个过程，直至将矩阵 A化为行阶梯形。4二分法求方程的根4.1 问题背景在科学研究与工程计算中，常遇到方程(组)求根问题。若干个世纪以来，工程师 和数学家花了大量时用于探索求解方程 (组)，研究各种各样的方程求解方法。对于方程 f(x)=0,当f(x)为线性函数时，称f(x)=0为线性方程；当f(x)为非线性函数时，称式 f(x)=0为非线性方程。对于线性方程(组)的求解，理论与数值求法的成果丰富；对于 非线性方程的求解，由于f(x)的多样性，尚无一般的解析解法。当 f(x)为非线性

18、函数 时，若f(x)=0无解析解，但如果对任意的精度要求，设计迭代方程，数值计算出方程 的近似解，则可以认为求根的计算问题已经解决，至少能够满足实际要求。4.2 数学模型使用二分法求方程x”+x-1=0在0,1内的近似根(误差10A-5)。二分法：二分法是最简单的求根方法，它是利用连续函数的零点定理，将含根区间 逐次减半缩小，取区间的中点构造收敛点列 xk来逼近根x。4.3 计算方法二分法程序代码：function y=erfen1(m,n,er) syms x xka=m;b=n;k=0;ff=xA3+x-1;while b-a>erxk=(a+b)/2;fx=subs(ff,x,xk

19、);fa=subs(ff,x,a);k=k+1;if fx=0y(k)=xk;break;elseif fa*fx<0b=xk;elsea=xk;endy(k)=xk;endplot(y, '.-');grid on在命令窗口下执行:>> ab=ezfenl 1, Le-5): vpKabnS)实验结果如下：可以得到迭代区间中点数列分布及图像，数值如下：ans =0.5 , 0.75, 0.625, 0.6875, 0.65625, 0.671875, 0.6796875。0.68359375, 0.68164062, 0.68261719, 0.682128

20、91, 0.68237305, 0.68225098, 0.68231201, 0.68234253, 0.68232727, 0.6823349依据题目要求的精度，则需彳十七次，由实验数据知x=0.6823349即为所求的根。4.4 二分法的收敛性二分法算法简单，容易理解，且总是收敛的；但是二分法收敛速度太慢，浪费时间， 并且不能求复根跟偶数重根。5雅可比迭代求解方程组5.1 问题背景在自然科学和工程技术中有很多问题的解决常常归结为解线性方程组，例如电学中的网络问题，化学中的配平方程式模型问题，船体数学放样中建立三次样条函数问题，用最小二乘法求实验数据的曲线拟合问题，非线性方程组求解问题，用

21、差分法或者有限 元法解常微分方程、偏微分方程边值问题等都导致求解线性方程组。在实践中，通常采 用数值方法来讨论线性方程组 Ax=b的解，其中迭代法是一种重要方法。5.2 数学模型5.2.1 理论基础雅可比迭代法：首先将方程组中的系数矩阵 A分解成三部分，即：A = L+D+U其中D为对角阵，L为下三角矩阵，U为上三角矩阵。之后确定迭代格式，XA(k+1) = B*XA(k)+f,(这里A表示的是上标，括号内数字即迭代次数)，其中B称为迭代矩阵，雅克比迭 代法中一般记为J (k= 0,1 ,)，再选取初始迭代向量 XA(0),开始逐次迭代。设Ax= b,其中A=D+L+UJ非奇异矩阵，且对角阵

22、D也非奇异，则当迭代矩阵J的谱半径p (J)<1时，雅克比迭代法收敛0?2?3?0?2？?+0？我?10 对于Ax=b,A=?01? ?+ ?1 ?2?1?2 ?3 0=L+D+U因为？?w 0(Jacob 假设)所以D可逆。故有：(L+D+U)x=bDx=-(L+U)x+bx=-D-1 (L+U)+D-1b5.2.2 实例用雅可比迭代法解方程组：430?2435-1 ?= 30 0 -14? -245.3 计算方法雅可比迭代法程序：function x,k,index=Jacobi(A,b,ep,it_max)if nargin<4it_max=100;end if nargin

23、<3 ep=1e-5;endn=length(A);k=0;x=zeros(n,1);y=zeros(n,1);index=1;while k<=it_max for i=1:n if abs(A(i,i)<1e-10index=0;return ; endy(i)=(b(i)-A(i,1:n)*x(1:n)+A(i,i)*x(i)/A(i,i); end if norm(y-x,inf)<epbreak; endk=k+1;x=y; end在命令窗口输入的命令和结果如下图：>> A=H 35 -I;0 -1 4 ;b=24 33 -24J ;ep=le-5

24、;it_wiaz=l00;xT kj indss x = J ac ob i (Aj b, 0p,'x 二4. 2000 2.4000 -b. 40M39index =15.4收敛性分析由上面运行的结果可知方程组的解为4.2000 , 2.4000 , -5.4000,迭代次数为39, 由index=1表示计算成功。雅克比迭代法的优点明显，计算公式简单，每迭代一次只需计算一次矩阵和向量的乘法，且计算过程中原始矩阵 A始终不变，比较容易并行计算。然而这种迭代方式收敛 速度较慢，而且占据的存储空间较大，所以工程中一般不直接用雅克比迭代法，而用其 改进方法。与逐次超松弛迭代法相比，雅可比迭代

25、法收敛速度相对较慢，逐次超松弛迭代法收 敛速度较快。由逐次超松弛迭代法求出的方程组的数值解与该方程组的精确解十分接近。 并且离散化后线性方程组的逐次超松弛迭代法的精确性较高，故相对于雅可比迭代法， 逐次超松弛迭代法更加广泛地应用于实际， 可以用逐次超松弛迭代法求解高阶稀疏线性 方程组。6 Romberg求积法6.1 问题背景复化求积方法对于提高精度是行之有效的方法，但复化公式的一个主要缺点在于要 事先估计出部长。若步长过大，则精度难于保证；若步长过小，则计算量又不会太大。 而用复化公式的截断误差来估计步长， 其结果是步长往往过小，而且f (刈和f(刈在 区间a，b上的上界M事估计是较为困难的。

26、在实际计算中通常采用变步长的方法，即 把步长逐次分半(也就是把步长二等分)，直到达到某种精度为止，这种方法就是Romberg 积分法的思想。6.2 数学模型：6.2.1 理论基础记：TN0) TN ,将区间N等分的梯形值。T* SN ,将区间N等分的Simpson；TN2) CN ,将区间N等分的Coteso TN3)RN ,将区间N等分的Romberg(k)由其可构造一个序列'N ',次序列称为Romberg字列，并满足如下递推关系:TiMm 3Nb af(ak 12N(2 k以上递推公式就是4kT2(N1)TN4k 1k 1)一,k1,2,LRomberg积分递推公式6.2

27、.2 实例利用Romberg积分法的思想，用龙贝格公式数值积分求函数L' ?密积分-1.26.3 计算方法龙贝格公式积分程序：function I,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=zeros(1,1);h=b-a;T(1,1)=(h/2)*(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(

28、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,findsym(sym),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+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;输入：» clea.r all;a=-L, 2：b=l. 2：fps=le-6:s±ep=rcinberg 1曰 jb,epm

29、)可以得到结果：I =0. 995326000000000st ep 二6.4 误差分析Romberg求积公式是具有八阶精度的算法，收敛且稳定，比梯形公式，Simpson公式以及Cotes公式收敛的快。龙贝格公式的余项为:?,?(?= / ?(?)?/,?=- ?+22(?+1 )( ?+2?)2? !(? ?2?+3 ?夕?+2 )(?)Romberg积分公式高速有效，易于编程，适合于计算机计算。但他有一个主要的缺点是，每当把积分区间对分后，就要对被积函数 ？?(?野算它在新分点处的值，而这些函 数值的个数是成倍增加的7曷法7.1 问题背景工程及物理中的许多实际问题需要计算矩阵的特征值及特征

30、向量，对于给定的n阶实矩阵A,当n很小时用传统的方法是可以的，但当 n稍大时，计算工作量将以惊人的 速度增大，并且由于计算存在误差，用方程(入I-A)x=0求解十分困难。在实际应用中， 有的时候不一定需要求出矩阵的全部特征值和特征向量，例如只需要求出按模最大的或最小的特征值。求这类特征值的方法，通常采用迭代法，其中两种是幕法和反幕法。7.2 数学模型7.2.1 理论基础设An有n个线性相关的特征向量v1, v2,，vn,对应的特征值1, 2,n, 满足| i| > | 2| n| ,因为v1 , v2,，vn为Cn的一组基，所以任给x(0) 0,x(0) aivii 1,所以有:n /

31、k (0)«k /、A x A (aivi)i 1nkai iVki 1nkai A vii 1n：觇1()kaivii 21若a10,则因1 知，当k充分大时Ak)x(0)1&V1 = CV1属入1的特征向量另一方面,记max(x)=x,其中|Xi| 二 | x| ，则当k充分大时,max(Akx(0)max( 1ka1v1)1k max(a1v1)Ak 1.(0)kT"k 11max(A x ) max( 1 a1v1)1 max(a1v1)若a0,则因舍入误差的影响，会有某次迭代向量在v1方向上的分量不为0,迭代 下去可求得入1及对应特征向量的近似值。7.2.

32、2 实例用幕法计算矩阵3-10?= -132023绝对值最大的特征值及相应的特征向量，取精度要求为e=10-4 07.3 计算方法幕法Matlab程序：function m,u,index=pow(A,ep,it_max) if nargin<3 it_max=100; end if nargin<2 ep=1e-5; end n=length(A);u=ones(n,1);index=0;k=0;m1=0;while k<=it_maxv=A*u;,i=max(abs(v);m=v(i);u=v/m;if abs(m-m1)<ep index=1; break ; e

33、ndm1=m;k=k+1;end在命令窗口输入和输出如下图所示:» -l 0 -1 3 2:0 2 31:Mj uf ind«xj -po¥ (A, 11)U1L =5,-0,4136 LOCCO 0,9113index 7.4 误差分析幕法程序可以用来计算矩阵绝对值最大的特征值及相应的特征向量。 幕法的缺点是 开始的时候并不知道矩阵是否有单一的主特征值，也不知道如何选择Vo以保证它关于矩 阵特征向量的表达中包含一个与主特征值相关的非零特征向量。采用幕法和反幕法，求矩阵的最大和最小特征值，从原理上看，这两种方法都是迭 代法，因此迭代初始向量的选择对计算结果会产生一定影响，主要表现在收敛速度上。同时，原点位移m的选取也影响收敛的速度，但原点位移mo的适当选取依赖于对矩 阵A的大致了解。8改