计算方法试验报告

1、实验报告一、求方程 f(x)=xA3-sinx-12x+1 的全部根,=le-61、用一般迭代法;2、用牛顿迭代法；并比较两种迭代的收敛速度。一、首先，由题可求得：/W =3x2-cosx-12 .其次，分析得到苴根所在的区间。1令f(x)-0,可得到%3-12x4-1 = sinx.2用一阶导数分析得到X3-12X+ 1和sinx两个函数的增减区间：再用二阶导数分析得到两个函数的拐点以及凹凸区间.3在直角坐标轴上描墓出X3-12X+1= 0和sinx = 0的图，在图上可以看到他们的交点, 然后估计交点所在的区间， 即是所要求的根的区间。 经过估计， 得到根所在的区间为-4-3,0,1和3,

2、4.1、一般迭代法(1)算法步骤：设e为给定的允许精度，迭代法的计算步骤为：1选定初值A().|h /(A) = 0确定函数g(x),得等价形式x = g(x).2计算gg).由迭代公式得X|=g(x。).3如果|x,-x0|,则迭代结束，取“为解的近似值；否则，用“代替儿，重 复步骤和步骤.(2)程丿了:代码：1在区间-4,-3内， 代码：clcx0=-35;$初值心iter_max=100;$迭代的最大次数ep=le-6;$允许精度k=0;while k=iter_max$从0开始至ljiter_max循坯；xl=(sin(x0)+12*x0-l) .(1/3);貂弋入儿，算出旺的值if

3、abs (xl-xO) ep %|x,-x0与允许精度作比较break;匕条件卜|一对5成立，跳出循环endxO=xl;务条件|Xj-XQ 不成立,用“代替兀0k=k+l;%kiniendx_star=xl, iter=k为解的近似值，iter为迭代次数运行结果：x_star = -3.4101 : iter =142在区间0,1内，代码：clcxO=O . 5;$初值旺iter_max=100;$迭代的最大次数ep=le-6;$允许精度gk=0;while k=iter_max$1从0开始至ljiter_max循坯;xl=(l/12) * (x0.3-sin(x0) +1);貂弋入心,算岀“

4、的值if abs (xl-xO) ep%|xj -x0|与允许精度作比较break;务条件x-x()| s成立，跳出循环endxO=xl;匕条件|x,-XQ不成立，用“代替儿k=k+l;Ek力mendx_star=xl, iter=k为解的近似值，iter为迭代次数结果:x_star = 0.07696；iter =63在区间3,4内，代码：clcx0=3.5;&初值心iter_max=100;$迭代的最大次数ep=le-6;$允许精度k=0;while k=iter_max靱从0开始至ljiter_max循环xl=(sin(x0)+12*x0-l)(1/3);訂弋入心，算出“的值if

5、 abs(xl-xO)ep引召-对与允许精度作比较break;匕条件-x()| 成立，跳岀循环endxO=xl;$条件|Xj -x0| 不成立，用X代替X。k=k+l;力mendx_star=xl, iter=k为解的近似值，iter为迭代次数运行结果：x_star =3.4101；iter =102、牛顿迭代法（1）算法步骤：1选泄初值计算/（“），/ （%（）2按公式心+|=兀一丄2迭代，得新的近似值兀+并计算/（林+J，/（A.+l）.fxk）3对于给立的允许精度 S 如果卜阳一兀卜,则终止迭代，取“ SJ；否则,k = k + ,在转回步骤计算.(2)程丿了代码：1在区间-4-3内，c

6、lcxl=-3.5;$初值“k=0;while k=100%k从0开始到100循环x0=xl;$将初值勺赋给儿f0=x0 .Zv3-sin(x0)-12*x0 + l;&计算/(心)fl=3*x0.2-cos (x0) -12;务计算/(x0) xl=xO-fO/fl;$计算得到新的近似值“if abs (xl-xO) 1. Oe-6 %| -x0|与允许精度作比较break;匕条件x-x0成立，跳出循环endk=k+l;务条件1%,-XQ不成立，k加1 endx_star=xl, iter=k%坷为解的近似值，iter为迭代次数运行结果：x_star =-3.4911；iter =2

7、2在区间0,1内，clcxl=0.5;%初值Xk=0;while k=100从0开始到100循环x0=xl;&将初值坷赋给儿f 0=x0 .A3-sin (x0) -12*x0+l;$计算/(x0)fl=3*x0./x2-cos (xO) -12;%i|-M/(x0)xl=x0-f0/fl;$计算得到新的近似值“if abs (xl-xO) 1. Oe-6-x与允许精度作比较break;匕条件-x0| 成立，跳出循环endk=k+l;%条件|Xj-XQ不成立，k加1endx_star=xl, iter=k为解的近似值，iter为迭代次数运行结果：x_star =0.07696；iter

8、 =33在区间3,4内，clcxl=3.5;%初值“k=0;while k=100靱从0开始到100循环x0=xl;$将初值X赋给几f 0=x0八3-sin (x0) -12*x0+l;%计算/(x0)fl=3*x0./s2-cos (xO) -12;$计算/(x0) xl=xO-fO/fl;咎计算得到新的近似值“if abs (xl-xO) 1. Oe-6%|xj-x0与允许精度作比较break;E条件x-xQ成立，跳出循环endk=k+l;匕条件-x0| 不成立，k加1endx_star=xlziter=k*山为解的近似值，iter为迭代次数运行结果：x_star =3.4101；iter

9、 =33.运行结果：xe-4,-3X e o,lX e 3,4x_stariterx_stariterx_stariter一般迭代3.4101140.769663.41011()牛顿法3.491120.7069633.410134.结果分析:从这题的结果可以看出，牛顿迭代法的迭代速度比一般迭代法的速度要快，牛顿法的迭 代次数比较少。例如，求在-4-3的根时，一般迭代法迭代了14次仃ter =14),而牛顿法只用了2次(iter =2).总而言之，一般迭代法是一种逐次逼近的方法，具有原理简单、编制程序方便等优点, 但存在是否收敛和收敛速度的快慢等问题。牛顿迭代法是一种特姝的迭代法，用于求方程的单

10、根时具有2阶收敛。因此是一种求非 线性方程解的好方法，还可以用于求重根和复根，而且可以推广到求非线性方程组的解.二. 解方程组直接法421872A =483对矩阵A做LU分解。K杜利特尔分解法(直接三角分解法)：设方程组心=/2的系数矩阵A的各阶主子式det(AjO(A: = 1,2,则可知该方程组存在唯一的杜利特尔分解：A = LU,其中L =1,U =UliU12“Jll2n，211仃2hin-L% _我们记矩阵A的第i行第丿列元素为“厂 则矩阵U的第一行和厶的第一列可山公式(1.1)求出:(2.1)再由公式(2.2)求出(/的第k行(&=2,/)、厶的第列(& =2,-1

11、)元素:2 1 -71 2 1-5* * b =1 2 1-51 2_9_-5_A =51062012 6 112.用追赶法解下述方程组Ax = b9并给出n=10的结果，其中厶1 =J = 2,3,/.（1）算法步骤：由于系数矩阵A的各阶主子式det（人）H0（k = l,2,时，则可知该方程组存在唯一的杜利特尔分解：A = LU 1利用公式（2.1）和（2.2）求出U和厶的元素和2求解单位下三角形方程组Ly = b,即按公式儿-工/灯儿,2,3,“23）求出儿2,，儿。4求解上三角形方程组Ux=y,即按公式可求得方程组的解：占,/(2)程序代码：function LzU=LU(A)Anzn

12、=size (A) ;%定义A为/阶矩阵L=zeros (n, n) ;%矩阵厶初始化(2.2)(2.4)儿-m内；-+i%U=zeros (n, n) ;$矩阵U初始化for i=l:nL(izi)=l；$矩阵厶的对角线元素为1X】=11V1X2 endfor k=l:nfor j=k:nU(kzj)=A(kzj)-sum(L(kz1:k-1) .*U(1:k-1,j)1);匕求出u元素endfor i=k+l:nL(izk) = (A(i,k) -sum(L(i, 1 :k-l) .*U(l:k-lzk)1) ) /U(k,k);E求出厶的元素 gendendA=4 2 1 5; 8 7

13、2 10; 4 8 3 6;12 6 11 20 ;&矩阵AL,U=LU(A)A = LU2追赶法设所求方程组为bxx+cx2= ra2x+b2x2+c2x3= r2a3x2+ bg + c3x4= r3(1)算法步骤：1由方程组的第1个方程，将未知元“用七表示为,_ 八一 3i己M = V)=,得(2.5)+ bnXn=rn1以此类推，*2用心表示，无用林 4 表示。我们可以得到一般式（其中记上式右端恰好是式（2.6）的心。于是，可由式（2.7）求岀三对角方程组的解心，即（2）程序代码:clca=01111 1 11111;b=22222 2 2222；c= 11111 1 1110

14、；r=-1 -5 c5 -5-5 -5 -5 -5 -5 -5;u=00000 0 0000；v= 00000 0 0000；x=00000 0 000 0；u(l)=r(l)/b(l) ;&计算妁=v(l)=c(l)/b(l) ; E计算儿=乞for k=2:10u(k) = (r(k)-u(k-l)*a(k) )/ (b(k)-v(k-l)*a(k) ;&计算叭v(k)=c(k)/(b(k)-v(k-l) *a(k) ) ;$计算吸endx(10)=u (10); bk -vk=-，比=1,2,/-1Xg=k_vx如,k = l,2,，n_l(2.6)(2.7)2将畑=仏-

15、忖心代入方程组（2.5）的第n个方程,并解出耳得(2.8)for k=l:9x (k) =u (k) -v (k) *x (k+1) ;xk=uk,解出x如endx运行结果:x=-3.5000 -1.0000-3.0000 -1.6000-2.8333 -1.8571 -2.7500-2.0000-0.8182 -2.09093、运行结果：（1）杜利特尔分解法：5、2 15、2 13 0 0L =U =1 2 12 1 Jacobi（雅可比）迭代法：（1）算法步骤：1根据矩阵A ,算出矩阵-厶、D = diag(sd22,2计算雅可比迭代矩阵乞：Bj =Dl（L + U）.3计算矩阵fj:fj

16、=Db.4代入给出的初值*）=（却）,雋）,普）7,得到雅可比迭代公式（3.1）的矩阵 形式：xm+=BjXm+ fj ,加= 0,1,.（3.2）由式（3.2）就可以得到向量序列：Jo） J0 .M .（2）程序代码:clcA=3 -1/2 -1/4 0 0 0 0 0 0 0;-1/2 3 -1/2 -1/4 0 0 0 0 0 0; -1/4 -1/2 3 -1/2 -1/4 0000 0;0 -1/4 -1/23 -1/2 -1/4 0 0 0 0;0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0;0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0;,-U =

17、知03a23a2n00 .0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0;0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4;0 0 0 0 0 0 -1/4 -1/2 3 -l/2;0 0 0 0 0 0 0 -1/4 -1/2 3;x0=0 00000000 01;$ 初值x二叶),),乂o)7b=5 555555 5 5 51;L= 000000000 0;l/2 000000000;1/41/20000000 0;0 1/41/2000000 0;0 01/41/200000 0;0 001/41/20000 0;0 0001/41/2000 0;0 000

18、01/4 1/2 000;0 000001/4 1/2 0 0;0 0000001/4 1/2 0;U= 01/2 1/400 0 0 0 0 0;0 0 1/2 1/4 0 0 0 0 0 0;0 001/21/40000 0;0 0001/21/4000 0;0 00001/2 1/4 0 0 0;0 000001/21/40 0;0 000000 1/2 1/4 0;0 00000001/2 1/4;0 00000000 l/2;0 000000000;D=3000000 0 0 0;03 0000 000 0；003000 0000；00 0300 0000；00 0030 0000；

19、00 00030000；00 00003000；00 0000 0300；00 0000 003 0；00 0000 000 3；BJ=i nv (D )* (L +U ）；匕讣算雅可比迭代矩阵乞=Dl（L + U）FJ=inv(D) *b;$计算矩阵刀=DbN=1000;ep=le-10;% 误差 =10_1k=0;while k=iter_max%1从0开始至ljiter_max循环xl=BJ*xO + f J;% 计算x(n,+1)=+ fjif norm( (xl-xO) , *inf *)ep列卜与误差 做比较break;$满足卜)匸S ,跳出循环end x0=xl;&不满足

20、|Xj-xQ - 皿尸（1一co）D + a）u）其中，e 为松弛因子，ew(-s,+oo)。3计算矩阵几：fj = co(D - coL)lb.其中，e 为松弛因子，ty e (-co,4o)o4代入给岀的初值*)=(#),雋) -,卅) ， 得到超松弛迭代公式的矩阵形式:W)= %Z)+几，加= 0,1,.山式(3.2)就可以得到向量序列：X()M （2）程丿了代码：clc0a2a35a2l0 %a2na3l。320,-U =: 5+.5匕20L =Q =伽g(,22，”仏)。A=3 -1/2 -1/4 0 0 0 0 0 0 0;-1/2 3 -1/2 -1/4 0 0 0 0 0 0;

21、 -1/4 -1/2 3 -1/2 -1/4 0000 0;0 -1/4 -1/23 -1/2 -1/4 0 0 0 0;xl=Bw*xO+fw;0 0 -1/4 -1/2 3 -1/2 -1/4 0 0 0;0 0 0 -1/4 -1/2 3 -1/2 -1/4 0 0;0 0 0 0 -1/4 -1/2 3 -1/2 -1/4 0;0 0 0 0 0 -1/4 -1/2 3 -1/2 -1/4;0 0 0 0 0 0 -1/4 -1/2 3 -l/2;0 0 0 0 0 0 0 -1/4 -1/2 3; x0=0 00000000 01;0；0；0；0；0；0；3；w=l .3;&

22、超松弛因子血Bw=(inv(D-w*L) ) * ( (1-w) *D+w*U) ;&计算超松弛迭代矩阵Fw=w* (inv(D-w*L) ) *b;&计算矩阵iter_max=1000;ep=le-10;E误差= 10,0k=0;while k=iter_maxif norm( (xl-xO) , *inf1)ep务卜_兀|L与误差做比较break;endxO=xl;k=k+l;endx_star=xl, iter=kL= 00000 00 000；l/2 0000 00 00 0f1/41/20 00 00000;01/41/20 00 0000；0 01/41/20 00

23、000;0001/41/20 0000；0 0001/41/20000;0000 01/41/2000；0 0000 01/41/200;0000 00 01/41/20；U= 01/2 1/4 00 00000；0 01/21/40 00000；0 001/2 1/4 00000；0 0001/21/40000；0 0000 1/2 1/4000；0 0000 01/21/400；0 0000 00 1/21/40；0 0000 00 01/21/4;0 0000 00 001/2;0 0000 00 000fD=30000 00 000；b=5 5 5 5 5 5 5 5 5 51;030

24、0000000；0030000000；3、运行结果：(1)雅可比迭代法：x, =(2.40792.8663 3.16213.2558 3.2929 3.2929 3.2558 3.16212.86632.4079)T迭代次数：iter =31(2)超松弛迭代法：X, =(24079 2.8663 3.1621 3.2558 3.2929 3.2929 3.2558 3.1621 2.8663 2.4079)T迭代次数：iter =254、结果分析：雅可比迭代法和超松弛迭代法从迭代次数上比较，可以看出超松弛迭代比雅可比迭代的次数少，表明超松弛迭代法比雅可比迭代法的收敛速度快。四、插值逼近题目：/

25、 = -r-ve-5,51,将-5,5 10等分，作Lagrange插值，将1 + X插值函数的图形与y = /(x)的图形比较，并给出结论。nX X*=1;=1XkXj1、算法步骤：已知y = Jx 求出各个节点的拉格朗日插值基函数多项式:! = (x一X。X一J(x一兀_1 X一召+J(兀一x”)= A %-兀1U - “ X刃- %J(召-Vi XN-和)(兀-兀)1=1 “-兀k和厶M就称为以儿为节点的拉格朗日插值基多项式。 利用基函数的线性组合得到插值多项式：厶(x)=yj()(x)+儿/“(x)。2程序代码:function y=lagrange(xO,yO,x)n=length(

26、xO);m=length(x);for i=l:mz=x(i);s=0.0;for k=l:np=1.0;for j=l:nifj=k p=p*(z-xO(j)/(xO(k)-xO(j);%求出厶=fj _，p即为厶jt=OXj兀女k*iendends=p*yO(k)+s;%计算Ln(x)=儿/()(x)+ +儿/”(x), s即为厶”(x)endy(i)=s；endx=-5:l:5;y=l./(l+x.A2);%插值函数x0=-5:0.1:5;yO=lagrange(x,y,xO);yl = l./(l +xO.A2);%绘制图形plot(xO,yO,-r)%插值曲线hold onplot(

27、xO,yl/-b,)%原曲线3、运行结果:插值曲线为红色，原曲线为蓝色:4、结果分析：当“很大时，在3.63内乙能很好的逼近f（x）,而在这个区间外，误 差反而很大。由图像我们可以看出，在-5,5内部，厶（x）与/（X）偏差小；而在 端点x = 5附近，偏差大。这就是龙格现象。五、复化求积题目：分别用复化梯形公式、复化辛卜生公式计算其中1心-忑I = 0.5X10-7.（用区间逐步分半递推算法）1、复化梯形公式:（1）算法步骤： 将区间肚切分成” =2如=0,1,2,）等分,步长为分点为xk= +也(k = 0,1,2,2)其中,我们乂题意已知：sb =1,2。2建立7；心与的递推公式：厶Q H +九土他+(2k-叽)(m = 1,2,)。3给定误差 s 当Tr, -Tvs时，停止计算，取G为所给积分的近似值。(2)程丿了;代码：clca=l