数值分析实验报告

1、天津工业大学数值分析实验报告电子与信息工程学院例题一用Gauss消去法来解方程组?1+ ?2+ ?3= 6 12?1- 3?2+ 3?3= 15-18?1 + 3?2- ?3= -15模型：高斯(GausS消去法高斯消去法的基本思想是： 先逐次消去变量， 将方程组化成同解的上三角形 方程组，此过程成为消元过程。 然后按方程相反顺序求解上三角方程组， 得到原 方程组的解，此过程称为回代过程。这种方法称为Gauss消去法。将方程组改成如下形式?1)?1+ ?2+ +? ?1)?坊)？1+ ?2)?2f +? ?1 :：(2- 1) ?；?1+ ?2?片 +?；? ?消元过程 ：第一步：设？刃工0,

2、记昭二？?/?；)(? 2,3, , ?)方程中第i个方程减去第一个方程乘以 ?1，完成第一次消元。将上述方程式可以改写为?(111)?1+ ?(1 12 ) ?2+?1(?1?)?=? ?1(?1)?2?(22) ?2+?2(2?)?=? ?2?(2)22 2? ： 2 ( 2-2 )?(?22)?2+?(2?)?=? ?(2?)其中 ?(?2?)?= ?(1?)?- ?1 1?(1?)? ， ?(2?) = ?(1?) - ?1?1?(1)(?,?= 2,3， ? 。)方程组(2-2)简记为 ?2?= ?2.第二步：设?2)工0,记？?2= ?22/?2)(?/= 2,3, , ?)方程中

3、第i个方程减去第一个方程乘以 ?1， 完成第二次消元 第 n-1 步：同上,得如下方程组?(111)?1+ ?(112)?2+?(11?)?=? ?(1?1 )?(222)?2+?(22?)?=? ?(22)? (?2?)? ?=? ?(?2?)回代过程：按变量的逆次序逐步回代得到方程组的解(k = n-1,n-2 ,1)?= ?/?(?)?= (?；?- 2?=?+1?)1/?；?算法：1. 输入 A=(?), b=(? ?), n，置 k=102. 若?k丰0,转3;否则输出失败信息，停机。3. 对，置i= k + 1 ，-n 置aik/a kk ? bi 对bj-?；?4. 若 k=n

4、-1,转 5;否则 k+1? k， 转 2,05. 若？需=0 ,输出失败信息，停机；否则，置 bn/a nn ? bn。6. 对 k = n-1,n-2 , 1 置bk - EJU+i akl bi)/a? bk7. 输出x=b停机代码： 步骤一：编写高斯消去法函数 function X =Fgauss(A ,b) zengguang =A b;n=length(b); ra=rank(A);rz=rank(zengguang); temp1=rz - ra;if temp10,disp( 错误 )returnendif ra=rzif ra=n;X=zeros(n,1); C=zeros(

5、1,n+1);for p= 1:n-1for k=p+1:n m=zengguang(k,p)/zengguang(p,p); zengguang(k,p:n+1)=zengguang(k,p:n+1)- m*zengguang(p,p:n+1);endendb=zengguang(1:n,n+1);A=zengguang(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)/A(q,q);endelsedisp( 错误 )endend步骤二：编写方程组信息a=1 1 112 -3 3-18 3 -

6、1;b=6 15 -15; x=Fgauss(a,b )步骤三：运行主函数 得到结果：X1=1.000X2=2.000X3=3.000结果分析：高斯消去法简单易行，但其计算过程中要求？?k工0,因而适用范围小，只能 适用于从1到n-1阶顺序主子式不为零的矩阵A,计算实践表明高斯消去法数值稳 定性差，当出现小主元素时，会严重影响计算精度，甚至导致错误的结果。例题二用雅克比（Jacob）迭代法求解下列方程组10x1 - x2 - 2x3 = 72-x1 + 10x2 - 2x3 = 83-x1 - x2 + 5x3 = 42模型：雅克比（Jacob）迭代法因为方程组aiixl + ai2x2+ +

7、a=nX n = bia2i x1 + a22x2 + +a2nXn = b2:（3- 1 ） an1 x1 + an2 x2 + +a nn xn = bn记的系数矩阵A非奇异，不妨设aii工0（i=1,2 ，n）。将方程组（3-1）变形为x2 = b21x1+ + 2nX n + g2（3-2 ） xn = bn1X1 + bn2X2 + + nn-1 Xn -1 + gn其中bij =-?，（i 工 j，i，j= 1，2，aii、bi，n），gi = （i =aii=1， 2，n)。若x1 =b12x2 + + 1nx n + g10bi2b13 1n-1 b1nbb220b3 2n-1

8、 b2nbB =bn1 bn2bn3 nn-10)g = g1 g2 nTg则方程组(3-2)可以简记为x = Bx + g (3-3)选取初始向量？?)代入迭代公式x(k+1) = Bx(k) + g(k = 0,1,2, ( 3-4)产生向量序列?%),由上述计算过程所给出的迭代法称为Jacobi迭代法，又叫简单迭代法。式(34)为它的迭代公式，迭代矩阵为B。这是Jacobi迭代法的计算公式。如果用系数矩阵A来表示，则有B = I - D-1 A，g= D-1 b (3-5)111其中 D= diag (an，a?2，ann )，D-1 = diag(-，-，a )。a11a22ann算法

9、：1输入 A=(aij), b=(b1，b2，bn)，维数 n, x = (x(0)，x；0， ，xf)，E,最大容许迭代次数N。2. 置 k=13. 对 i=1, 2, ，n xi = (bi -耳1 aij xj)/a ijj勺4. 若llx x(0) | ?输出x,停机；否则转5。5. 若 k=tolx0=x;x=B*x0+f;k=k+1;if (k=max1)disp( 错误 );return endk ,xend步骤二：编写方程组信息a=10 -1 -2;-1 10 -2;-1 -1 5; b=72 83 42;x0=0 0 0;x,k=Fjacobi(a,b,x0,1e-7) 步骤

10、三：运行主函数得到结果 jwwbi1 1. 000012. 000013. 0000结果分析：由表所示方程组的精确解为x= (11 12 13 )T当迭代次数增加时，迭代结果越来越 接近精确解。在迭代到11次后数值不变2.00009.710010.700011.5000|3.000010.570011.571012.48204.000010.853511.853412.82825.000010.951011.951012.94146.000010.983411.983412.98047.000010.994411.994412.99338.000010.998111.998112.99789.

11、000010.999411.999412.999210.000010.999811.999812.9997|11.000010.999911.999912.9999|12.000011.000012.000013.0000|13.000011.000012.000013.0000|14.000011.000012.000013.0000|15.000011.000012.000013.0000116.000011.000012.000013.0000|17.000011.000012.000013.0000118.000011.000012.000013.0000|?19.011.000012

12、.000013.0000例题三：已知函数表如下:X-1 0 1yi2X0.512利用线性插值和二次插值计算20.3的近似值模型：Lagrange插值函数y f x在区间a,b上n 1个互异节点x0,为，,xn的函数值y f(x), i 0,1,L , n。求一个至多 n次多项式 n(x) a。aix Lanxn使其在给定点处与f x同值，即满足插值条件。线性插值，即只有两个节点xo, xi（n 1）的插值多项式。插值多项式设为i（x） ao aix,且满足插值条件1 (xo)1(X1)a0 a1x0 y0f(x0)，解此方程组得ao CX1 y1f(xja0, a1X0x!也丄，所以两点的一次

13、插值多项式为X 人1（X）xX,y。-XX1xX0X1X(3-1)记 I（X）3l1(x)X0 X1，则（3-1 ）可写成1（x）X1 X0yl(x) y1h(x)这是用两点（x, y ）,（X1, yj的直线y1（x）近似曲线y f x ,故这种插值又称线性插值。由此得到启发，当节点增加到n 1时，可以先构造n次多项式li(x) (i 0,1,L , n)，它们满足 li(Xj)Xjj 0 XiXj插值基函数为 h(X) (X X0)L(X x1)(x Xi1)L(X Xn)(Xi X0)L (Xi X 1)(Xi Xi 1)L (Xi Xn)nnn X X插值多项式为 n(x)yili(x

14、)yi(-) (3-2)i 0i 0 j 0 Xi Xjj i式(3-2)称为n次Lagrange插值多项式，为了以后便于区别，常用Ln(x)代替n(x)以突出表示这是由Lagrange插值所得到的插值多项式，即(3-3)nLn(x)yili(x)i 0特别地，n=1时，两点一次Lagra nge插值(也叫线性插值)多项式为x x1 YoJ(x) yl(x)川&)X。 X1x X0-Y1 ； X1X0n=2时称为二次插值，也叫抛物线插值，插值多项式为L2(x) yl(x)%li(x)y2(x)x x-1x x?yo -xoX1x。X2X X0Xy1XI X0X-1x2X2x x0X x1y2-

15、X2 X0X2X1算法:(1)输入 x, xi, yi (i 0,1 丄，n)。(2)n X xj对 i 0,1,L , nljj 0 x Xj(3)(4)输出L f(x),停机。nL Yili-i 0程序:fun cti on yp=mlagr(x,y,xp)n=len gth(x);m=le ngth(xp); yp=zeros(1,m); c1= on es( n-1,1); c2=on es(1,m); for i=1: nxb=x(1:i-1,i+1: n);yp=yp+y(i)*prod(c1*xp-xb*c2)./(x(i)-xb*c2); end结果分析：在命令窗口输入x=0

16、1;y=1 2; xx=0.3;yy仁mlagr(x,y,xx)输出结果为yyi =1.300000000000000这是线性插值的求得的203的近似值。在命令窗口输入x=-1 0 1;y=0.5 1 2; xx=0.3;yy2=mlagr(x,y,xx)输出结果为yy2 =1.247500000000000这是抛物线插值求得的20.3的近似值。20.3的准确解为1.2311，可见抛物线插值比线性插值的误差要小。这是因为抛物线插值选取了比线性插值更多的节点求解近似值。由插值余项f(n1)() 亠 一、RJx) - - n 1(x)可知，Lagrange插值法要提咼插值多项式的函数逼近程度，(n

17、 1)!要增加节点数，提高多项式的次数，但这样往往不能达到预想的效果。为了克服 Runge现象，我们常用分段线性插值来改进。例题四:分别用复化梯形公式和复化Simpson公式计算下列积分的近似值01xTdx nx模型：复化梯形公式和复化Simpson公式（1）复化梯形公式用分段线性插值函数近似被积函数，等于把积分区间分成若干小区间，在每 个小区间上以梯形的面积近似曲边梯形的面积。 即用梯形公式求小区间上积分的 近似值。如图1图1将积分区间a,b n 等分，记 h -_-， xk a kh(k 0,1,., n)nbfax dxn 1Xk1fk 0 xkn 1 hx dxf xkk 0 2f X

18、k 1整理得bhn 1f x dxfa f b2f xkTn（1）a2k 1式（1）称为复化梯形公式。在每个小区间Xk,Xk1（k0,1,., n 1）上用梯形公式求和，得如果f xC 2 a,b，在小区间X,xk 1上，梯形公式的截断误差为xk 1hr上h3f x dxf Xkf Xk 1fk ,kxk, xk 1xk212因此bh3 n1Rtff x dx Tnfka12k 0因为fx在a,b连续，有介值定理，存在a,b ,使得1n 1从而有bhb a Rt f f x dx Tnnfh2 fa, b（2）a 12 12这就是复化梯形公式的截断误差。（2）复化Simpson公式如果用分段而

19、插值函数近似被积函数，即在小区间上用Simpson公式计算积分近似值，就导出复化Simpson公式。将区间a,b分成n 2m等分，分点为xk a kh k 0,1,L ,n , h b_a，在n每个小区间x2k2,x2k上，用Simpson公式求积分，则有2kf xdxX2kX2k 2 fx2k 24 f x2k 1fX2k2k 26h3f x2k 24fx2k 1f x2k求和得bmx?km hf xdxf x dxf x2k 2 4f X2k 1f X2kak 1x2k 2k 1 3整理后得到bhm 1mfx dx1 1f afb 2f X2k4 f X2k 1（3）a3k 1k 1式（3

20、）称为复化Simpson公式。如果f x C 4 a,b ,复化Simpson公式的截断误差为农fbf Xadx - f a f3m 1mb 2 f X2k4 f X2k 1k 1k 15m2h f 4k kX2 k 2 , X2 kk 12880因为f 4 x连续，故存在a, b ,使得1 m4ffkm k 1代入上式得52h上4b aRs fmfh f,a,b2880180时，用复化Simpson式（4）表明，步长h越小，截断误差越小。当n 2m 公式所求的近似值收敛域积分值，而且算法具有数值稳定性。算法:a.复化梯形公式(1)输入f(x) , a, b,区间等分数nbhn 1(4)f x

21、 dxf af b 2 f xka2k 1(3)xka kh(k 0,1，n)Tn(5)停机。b.复化 Simps on公式(1) 输入f(x) ，a，b，区间等分数n，n=2m(2)h b a, nXka kh(k0,1,., n)bhm 1m(3)f x dxhf a f b2f X2k4 f X2ka3k 1k 1(4)停机。程序：a.复化梯形公式fun cti on s=mtrap(f,a,b ,n) h=(b-a)/n;x=li nspace(a,b ,n+1);y=feval(f,x);s=0.5*h*(y(1)+2*sum(y(2:n )+y( n+1);endb.复化 Simp

22、s on公式fun cti on s=msimp(f,a,b ,n)h=(b-a)/n;x=li nspace(a,b,2* n+1); y=feval(f,x);s=(h/6)*(y(1)+2*sum(y(3:2:2* n-1)+4*sum(y(2:2:2* n)+y( n+1); end结果分析：在MATLA工作窗口输入 format lo ng f=(x)x./(1+x.A2);s=mtrap(f,0,1,8) p=msimp(f,0,1,8)输出结果为s =0.345268948437223P =0.3444908975388951 xs和p分别是复化梯形公式、复化 Simpson公式

23、对一dx在n 8的近似值01 x例题五确定方程x3 - x+ 4 = 0的实根的分布情况，并用二分法求在开 区间(-2 , -1 )内的实根的近似值，要求精度为0.001.模型：二分法二分法也称对分区间法。它的基本思想是：先确定方程f(x)=0含根的区间(a,b)，再把区间逐次二等分。设f(x)在区间a,b上连续，且有f( a) f( b) 0则由连续函数介值定理，f(x)在(a,b)内必有零点，称(a,b)为方程 f(x)=0的有根区间。不妨设f(a) 0,取x0 =号，若f(x0)= 0,则x = x0就是方程f(x)=0的解。否则，若f(x0) 0，取ai = a，bi = x，则有ai

24、, bi ?a,b， bi - aib-a2且f(x)在色，6上连续，满足f(ai)f(bi) 0。重复上述过程又b a得到区间a?，满足a2,b2?印，b? - a?二亍，且f(a2)f(b2) 0。如此继续下去，得到一个区间序列a,b ?ai，bi?a2，b2】？ an，bn?满足f(an)f(bn) 0，bn -為二第因而满足an,bn(n=i,2, ?)均为方程f(x)=0的有根区间。由区间套 定理，存在x?a,b，使得lim an = lim bn = x?nn qx且x=x?是方程f(X)=O的根。若取xn = an；b n作为根X?的近似值，则误 差为|X?- Xnl 呼=籍(1

25、-1)按上述方法求非线性方程f(x)=0的近似解称为对分区间法，式(1-1 ) 表明，只要f(x)连续，对分区间法总是收敛的。算法：(1)(2)置 n=0。(3)(4)(5)a+bXn =。若f(Xn) = 0或牛 ?则输出X?= Xn，停机；否则转若nN,转6;否则输出失败信息，停机。5。输入f(x) , a, b,误差限,最大容许迭代次数N。(6) 若f(a)f(Xn) 0,disp(注意：ya*yb0 ,请重新调整区间端点a和b),return end max1=-1+ceil(log(b-a)-log(abtol)/log(2); for k=1:max1+1ya=fun(a);yb=

26、fun(b);x=(a+b)/2;yx=fun(x);wuca=abs(b-a)/2;k=k-1;k,a,b,x,wuca,ya,yb,yxif yx=0a=x;b=x;elseif yb*yx0b=x;yb=yx ；elsea=x;ya=yx;endif b-aabtolreturnendend k=max1; x;wuca;yx=fu n( x);end k,x,wuca,yx=erfe n(-2,-1,0.001)结果分析：首先先确定方程x3 - x+ 4 = 0的根的分布情况，在MATLAB工作窗口输入程序x=-4:0.1:4;y=x.A3-x+4;Plot(x,y)grid,gtex

27、t(y二x.A3-x+4)画出函数f(x)=x3 - x+ 4的图像，如图1-1所示，从图像中可以看出， 此曲线有两个驻点士血都在x轴的上方，在(-2 , -1 )内曲线与x轴3只有一个交点，则该方程在(-2 , -1 )内唯一的一个实根其次，再在区间(-2 , -1 )内求出实根的近似值，在MATLAB的工作窗口输入程序：k,x,wuca,yx二erfe n -2,-1,0.001运行后屏幕上显示二分法的计算结果，如表1-2所示：表1-2次数k左端点ak右端点bk中点Xkbn - an| 2 |函数值伽)函数值f(bk)函数值f(Xk)0-2.0000-1.0000-1.50000.5000

28、-2.00004.00002.12501-2.0000-1.5000-1.75000.2500-2.00002.12500.39062-2.0000-1.7500-1.87500.1250-2.0000.3906-0.71683-1.8750-1.7500-1.81250.0625-0.71680.3906-0.14184-1.8125-1.7500-1.78130.0313-0.14180.39060.12965-1.8125-1.7813-1.79690.0156-0.14180.1296-0.00486-1.7969-1.7813-1.78910.0078-0.00480.12960.0

29、6277-1.7969-1.7891-1.79300.0039-0.00480.06270.02908-1.7969-1.7930-1.79490.0020-0.00480.02900.01219-1.7969-1.7949-1.79590.0010-0.00480.01210.0037经过计算，我们可以得出方程根的近似值为x=-1.7959,二分区间的次数为k=9，误差为wuca=9.7656e004。可见，二分法的优点是对函数的要求低，方法简单、可靠、程序 设计容易，可以容易估计计算次数，收敛速度恒定；缺点是不能求出 方程的偶重根，收敛速度慢。例题六 用迭代法求方程f(x) = x2 +

30、2x - 10的一个正根模型：迭代法迭代法的基本思想是利用某种迭代公式，使某个近似根逐步精确 化，直到满足精度要求的近似根为止。把给定的方程f(x)=0改写成x = (x)(其中(X)称为迭代函数)， 选择合适的初始值X。，代入林X)算得的结果记为xi =林Xo)，般 x1工X。；把乂1代入(x)算得的结果记为x2 = (x1), ?.若把xk代入 X)，般有迭代公式Xk+1=Xk)(k=0,1,2,?)(2-1)由此得到序列Xo , x1 , ? ,xk ,?.若迭代序列xk收敛到X?(及 lim Xk = x?),则当函数(x)连续时，由(2-1)得klim xk = limkk(Xk)=

31、 (怀 Xk)X? = (x?)，称为的不动点，即x?为原方程f(x)=0的根。一般计算 有限步，即可得到某种精度的近似根Xk。这种求方程根的方法称为简 单迭代法。算法：(1) 输入(x)，初始值X。，误差限,最大容许迭代次数N。(2) 置 n=1。(3) 计算 x= (x0)。(4) 若|x- Xo| 1)&(xdpiancha0.5)&(k3)disp(请用户注意：此迭代序列发散，请重新输入新的迭代公式)return;end if(piancha0.001)&(xdpiancha3) disp(祝贺您！此迭代序列收敛，且收敛速度较快) return;endp=(i-1) piancha xdpiancha xk;结果分析:我们容易知道方程f(x) = x2 + 2x - 10的一个正根x?(2,3)，当然 利用二次方程求根公式不难得到X? =2.316624。下面我们构造不同的 迭代函数(X)，再以Xo =2为初始值，分别用下面三种迭代法求解。x=奶(x) = (10 - x2)?2;x= 2(x) = 10?(x + 2);x= 3(x) = x- (x2 + 2x- 10)?(2x+ 2).(1) 当奶(x) = (10 - x2)?2时我们在MATLAB窗口输入程序k,pia ncha,xdpia