计算方法实验报告

《计算方法实验报告》由会员分享，可在线阅读，更多相关《计算方法实验报告（22页珍藏版）》请在装配图网上搜索。

1、实验报告一、求方程f(x)=x3-sinx-12x+1的全部根, =1e-61、 用一般迭代法; 2、 用牛顿迭代法; 并比较两种迭代的收敛速度。一、 首先，由题可求得：.其次，分析得到其根所在的区间。 令，可得到. 用一阶导数分析得到和两个函数的增减区间；再用二阶导数分析得到两个函数的拐点以及凹凸区间. 在直角坐标轴上描摹出和的图，在图上可以看到他们的交点，然后估计交点所在的区间，即是所要求的根的区间。经过估计，得到根所在的区间为，和.1、 一般迭代法 （1）算法步骤：设为给定的允许精度，迭代法的计算步骤为： 选定初值.由确定函数，得等价形式. 计算.由迭代公式得. 如果，则迭代结束，取为解

2、的近似值；否则，用代替，重复步骤和步骤.（2）程序代码： 在区间内，代码： clcx0=-3.5; %初值iter_max=100; %迭代的最大次数ep=1e-6; %允许精度 k=0;while k=iter_max %k从0开始到iter_max循环 x1=(sin(x0)+12*x0-1).(1/3); %代入，算出的值 if abs(x1-x0)ep %与允许精度作比较 break; %条件成立，跳出循环 end x0=x1; %条件不成立，用代替 k=k+1; %k加1endx_star=x1, iter=k %为解的近似值，iter为迭代次数运行结果：x_star = -3.41

3、01 ；iter =14在区间内， 代码：clcx0=0.5; %初值iter_max=100; %迭代的最大次数ep=1e-6; %允许精度k=0;while k=iter_max %k从0开始到iter_max循环x1=(1/12)*(x0.3-sin(x0)+1); %代入，算出的值 if abs(x1-x0)ep %与允许精度作比较 break; %条件成立，跳出循环 end x0=x1; %条件不成立，用代替 k=k+1; %k加1endx_star=x1, iter=k %为解的近似值，iter为迭代次数结果：x_star = 0.07696；iter =6在区间内，代码：clcx

4、0=3.5; %初值iter_max=100; %迭代的最大次数ep=1e-6; %允许精度k=0;while k=iter_max %k从0开始到iter_max循环x1=(sin(x0)+12*x0-1).(1/3); %代入，算出的值 if abs(x1-x0)ep %与允许精度作比较 break; %条件成立，跳出循环 end x0=x1; %条件不成立，用代替 k=k+1; %k加1endx_star=x1, iter=k %为解的近似值，iter为迭代次数运行结果：x_star = 3.4101 ；iter =102、 牛顿迭代法（1）算法步骤： 选定初值，计算，. 按公式迭代，得

5、新的近似值，并计算，. 对于给定的允许精度，如果，则终止迭代，取；否则，在转回步骤计算.（2）程序代码：在区间内，clcx1=-3.5; %初值k=0;while k=100 %k从0开始到100循环x0=x1; %将初值赋给 f0=x0.3-sin(x0)-12*x0+1; %计算f1=3*x0.2-cos(x0)-12; %计算x1=x0-f0/f1; %计算得到新的近似值if abs(x1-x0) 1.0e-6 %与允许精度作比较break; %条件成立，跳出循环endk=k+1; %条件不成立，k加1end x_star=x1, iter=k %为解的近似值，iter为迭代次数运行结果

6、：x_star = -3.4911；iter =2在区间内，clcx1=0.5; %初值k=0;while k=100 %k从0开始到100循环x0=x1; %将初值赋给f0=x0.3-sin(x0)-12*x0+1; %计算f1=3*x0.2-cos(x0)-12; %计算x1=x0-f0/f1; %计算得到新的近似值 if abs(x1-x0) 1.0e-6 %与允许精度作比较break; %条件成立，跳出循环endk=k+1; %条件不成立，k加1end x_star=x1, iter=k %为解的近似值，iter为迭代次数运行结果：x_star =0.07696 ；iter =3在区间

7、内，clcx1=3.5; %初值k=0; while k=100 %k从0开始到100循环x0=x1; %将初值赋给f0=x0.3-sin(x0)-12*x0+1; %计算f1=3*x0.2-cos(x0)-12; %计算 x1=x0-f0/f1; %计算得到新的近似值if abs(x1-x0) 1.0e-6 %与允许精度作比较break; %条件成立，跳出循环endk=k+1; %条件不成立，k加1end x_star=x1, iter=k %为解的近似值，iter为迭代次数运行结果：x_star =3.4101；iter =33、运行结果：x_stariterx_stariterx_sta

8、riter一般迭代3.4101140.769663.410110牛顿法3.491120.7069633.410134、结果分析： 从这题的结果可以看出，牛顿迭代法的迭代速度比一般迭代法的速度要快，牛顿法的迭代次数比较少。例如，求在的根时，一般迭代法迭代了14次（iter =14），而牛顿法只用了2次（iter =2）。总而言之，一般迭代法是一种逐次逼近的方法，具有原理简单、编制程序方便等优点，但存在是否收敛和收敛速度的快慢等问题。牛顿迭代法是一种特殊的迭代法，用于求方程的单根时具有2阶收敛。因此是一种求非线性方程解的好方法，还可以用于求重根和复根，而且可以推广到求非线性方程组的解.二、解方程组

9、直接法1、已知对矩阵A做LU分解。2、用追赶法解下述方程组，并给出n=10的结果，其中，1、杜利特尔分解法（直接三角分解法）：设方程组的系数矩阵的各阶主子式，则可知该方程组存在唯一的杜利特尔分解：，其中，.我们记矩阵的第行第列元素为，则矩阵的第一行和的第一列可由公式（1.1）求出： （2.1）再由公式（2.2）求出的第行、的第列元素： （2.2）（1）算法步骤： 由于系数矩阵的各阶主子式，则可知该方程组存在唯一的杜利特尔分解：。 利用公式（2.1）和（2.2）求出和的元素和。 求解单位下三角形方程组，即按公式 （2.3） 求出。 求解上三角形方程组，即按公式 （2.4） 可求得方程组的解：。（

10、2）程序代码：function L,U=LU(A)An,n=size(A); %定义为阶矩阵L=zeros(n,n); %矩阵初始化U=zeros(n,n); %矩阵初始化for i=1:n L(i,i)=1; %矩阵的对角线元素为1endfor k=1:n for j=k:n U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j);%求出元素 end for i=k+1:n L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)/U(k,k); %求出的元素 endendA=4 2 1 5;8 7 2 10;4 8 3 6;12 6 11

11、 20; %矩阵L,U=LU(A) %2、追赶法 设所求方程组为， （2.5）（1）算法步骤： 由方程组的第1个方程，将未知元用表示为 记，得。 以此类推，用表示， 用表示。我们可以得到一般式（其中记）： （2.6） （2.7） 将代入方程组（2.5）的第n个方程，并解出得， 上式右端恰好是式（2.6）的。于是，可由式（2.7）求出三对角方程组的解，即， （2.8）（2）程序代码：clca=0 1 1 1 1 1 1 1 1 1;b=2 2 2 2 2 2 2 2 2 2;c=1 1 1 1 1 1 1 1 1 0;r=-7 -5 -5 -5 -5 -5 -5 -5 -5 -5;u=0 0 0

12、 0 0 0 0 0 0 0;v=0 0 0 0 0 0 0 0 0 0;x=0 0 0 0 0 0 0 0 0 0;u(1)=r(1)/b(1); %计算v(1)=c(1)/b(1); %计算for k=2:10 u(k)=(r(k)-u(k-1)*a(k)/(b(k)-v(k-1)*a(k); %计算 v(k)=c(k)/(b(k)-v(k-1)*a(k); %计算endx(10)=u(10); for k=1:9 x(k)=u(k)-v(k)*x(k+1); %由，解出end x运行结果：x=-3.5000 -1.0000 -3.0000 -1.6000 -2.8333 -1.8571

13、-2.7500 -2.0000 -0.8182 -2.09093、运行结果：（1）杜利特尔分解法：（2）追赶法： 。4、结果分析：（1）杜利特尔分解法：在实现分解后，解具有相同系数矩阵的方程组，非常方便，每解一个方程只需用式（2.3）和式（2.4）解两个三角形方程组即可，大大减少了运算工作量。另外，矩阵和的元素可采用紧凑格式存放在系数矩阵的相应元素位置上，节省了存储单元。（2）追赶法：追赶法仅适用于三对角方程组的求解。我们由式（2.6）和式（2.8）可以算出，追赶法的乘、除运算次数仅为，加、减运算次数仅为。故追赶法具有运算量小和存储量小的优势。三、解方程组迭代法1、用迭代法解Ax=b，其中b=

14、(5，5，5)T，给定误差，用Jacobi和SOR两种迭代法计算，并给出n=10的结果。对于阶方程组，假定系数矩阵的对角元时，可得雅可比迭代格式为：， （3.1）记为方程组的系数矩阵的对角元组成的对角阵，和分别为由的元素组成的严格下三角阵和上三角阵，则系数矩阵可分解为：.1、Jacobi（雅可比）迭代法：（1）算法步骤： 根据矩阵，算出矩阵、和，即，。 计算雅可比迭代矩阵：. 计算矩阵：. 代入给出的初值，得到雅可比迭代公式（3.1）的矩阵形式：. (3.2)由式（3.2）就可以得到向量序列：.（2）程序代码： clcA=3 -1/2 -1/4 0 0 0 0 0 0 0;-1/2 3 -1/

15、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 -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;0 0 0 0 0 0 0 -1/4 -1/2 3;x0=0 0 0 0 0 0 0 0 0 0; %初值 b=

16、5 5 5 5 5 5 5 5 5 5;L= 0 0 0 0 0 0 0 0 0 0;1/2 0 0 0 0 0 0 0 0 0; 1/4 1/2 0 0 0 0 0 0 0 0;0 1/4 1/2 0 0 0 0 0 0 0; 0 0 1/4 1/2 0 0 0 0 0 0;0 0 0 1/4 1/2 0 0 0 0 0; 0 0 0 0 1/4 1/2 0 0 0 0;0 0 0 0 0 1/4 1/2 0 0 0; 0 0 0 0 0 0 1/4 1/2 0 0;0 0 0 0 0 0 0 1/4 1/2 0;U= 0 1/2 1/4 0 0 0 0 0 0 0;0 0 1/2 1/4

17、0 0 0 0 0 0; 0 0 0 1/2 1/4 0 0 0 0 0;0 0 0 0 1/2 1/4 0 0 0 0; 0 0 0 0 0 1/2 1/4 0 0 0;0 0 0 0 0 0 1/2 1/4 0 0; 0 0 0 0 0 0 0 1/2 1/4 0;0 0 0 0 0 0 0 0 1/2 1/4; 0 0 0 0 0 0 0 0 0 1/2;0 0 0 0 0 0 0 0 0 0;D=3 0 0 0 0 0 0 0 0 0; 0 3 0 0 0 0 0 0 0 0; 0 0 3 0 0 0 0 0 0 0; 0 0 0 3 0 0 0 0 0 0; 0 0 0 0 3 0

18、0 0 0 0; 0 0 0 0 0 3 0 0 0 0; 0 0 0 0 0 0 3 0 0 0; 0 0 0 0 0 0 0 3 0 0; 0 0 0 0 0 0 0 0 3 0; 0 0 0 0 0 0 0 0 0 3;BJ=inv(D)*(L+U); %计算雅可比迭代矩阵FJ=inv(D)*b; %计算矩阵N=1000;ep=1e-10; % 误差 k=0;while k=iter_max %k从0开始到iter_max循环x1=BJ*x0+fJ; %计算if norm(x1-x0),inf)ep %与误差做比较break; %满足，跳出循环endx0=x1; %不满足，将赋给k=k+

19、1; %k加1endx_star=x1, iter=k 2、SOR（超松弛）迭代法：（1）算法步骤： 根据矩阵，算出矩阵、和，即，。 计算超松弛迭代矩阵：. 其中，为松弛因子，。 计算矩阵：. 其中，为松弛因子，。 代入给出的初值，得到超松弛迭代公式的矩阵形式：. 由式（3.2）就可以得到向量序列：.（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 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

20、 -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 -1/2;0 0 0 0 0 0 0 -1/4 -1/2 3;x0=0 0 0 0 0 0 0 0 0 0; b=5 5 5 5 5 5 5 5 5 5;L= 0 0 0 0 0 0 0 0 0 0;1/2 0 0 0 0 0 0 0 0 0; 1/4 1/2 0 0 0 0 0 0 0 0;0 1/4 1/2 0 0

21、0 0 0 0 0; 0 0 1/4 1/2 0 0 0 0 0 0;0 0 0 1/4 1/2 0 0 0 0 0; 0 0 0 0 1/4 1/2 0 0 0 0;0 0 0 0 0 1/4 1/2 0 0 0; 0 0 0 0 0 0 1/4 1/2 0 0;0 0 0 0 0 0 0 1/4 1/2 0;U= 0 1/2 1/4 0 0 0 0 0 0 0; 0 0 1/2 1/4 0 0 0 0 0 0; 0 0 0 1/2 1/4 0 0 0 0 0; 0 0 0 0 1/2 1/4 0 0 0 0; 0 0 0 0 0 1/2 1/4 0 0 0; 0 0 0 0 0 0 1/2

22、 1/4 0 0; 0 0 0 0 0 0 0 1/2 1/4 0; 0 0 0 0 0 0 0 0 1/2 1/4; 0 0 0 0 0 0 0 0 0 1/2; 0 0 0 0 0 0 0 0 0 0;D=3 0 0 0 0 0 0 0 0 0; 0 3 0 0 0 0 0 0 0 0; 0 0 3 0 0 0 0 0 0 0; 0 0 0 3 0 0 0 0 0 0; 0 0 0 0 3 0 0 0 0 0; 0 0 0 0 0 3 0 0 0 0; 0 0 0 0 0 0 3 0 0 0; 0 0 0 0 0 0 0 3 0 0; 0 0 0 0 0 0 0 0 3 0; 0 0 0

23、0 0 0 0 0 0 3;w=1.3; %超松弛因子Bw=(inv(D-w*L)*(1-w)*D+w*U); %计算超松弛迭代矩阵Fw=w*(inv(D-w*L)*b; %计算矩阵iter_max=1000;ep=1e-10; %误差k=0;while k=iter_maxx1=Bw*x0+fw; %if norm(x1-x0),inf)ep %与误差做比较break;endx0=x1; k=k+1;endx_star=x1, iter=k3、运行结果：（1）雅可比迭代法：迭代次数：iter =31（2）超松弛迭代法：迭代次数：iter =254、结果分析：雅可比迭代法和超松弛迭代法从迭代次

24、数上比较，可以看出超松弛迭代比雅可比迭代的次数少，表明超松弛迭代法比雅可比迭代法的收敛速度快。四、插值逼近题目: ，将10等分，作Lagrange插值，将插值函数的图形与的图形比较，并给出结论。1、算法步骤： 已知， 求出各个节点的拉格朗日插值基函数多项式：， 就称为以为节点的拉格朗日插值基多项式。 利用基函数的线性组合得到插值多项式：。2、程序代码：function y=lagrange(x0,y0,x)n=length(x0);m=length(x);for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j=k p=p*(z-x0(j

25、)/(x0(k)-x0(j); %求出，p即为 end end s=p*y0(k)+s; %计算，s即为 end y(i)=s;endx=-5:1:5;y=1./(1+x.2); %插值函数x0=-5:0.1:5;y0=lagrange(x,y,x0);y1=1./(1+x0.2); %绘制图形plot(x0,y0,-r) %插值曲线hold onplot(x0,y1,-b) %原曲线 3、 运行结果：插值曲线为红色，原曲线为蓝色：4、结果分析：当很大时，在内能很好的逼近，而在这个区间外，误差反而很大。由图像我们可以看出，在内部，与偏差小；而在端点附近，偏差大。这就是龙格现象。五、复化求积题目

26、: 分别用复化梯形公式、复化辛卜生公式计算,其中. (用区间逐步分半递推算法) 1、复化梯形公式：（1）算法步骤： 将区间分成等分，步长为，分点为.其中，我们又题意已知：。 建立与的递推公式：。 给定误差，当时，停止计算，取为所给积分的近似值。（2）程序代码：clc a=1;b=2;m=1;h=0.5;ep=0.5e-7;f(a)=exp(1);f(b)=2*exp(2);x0=h*(f(a)+f(b);iter_max=100;i=0;while i=iter_max k=1; F=0; while k=2(m-1) F=F+(a+(2*k-1)*h)*exp(a+(2*k-1)*h); k

27、=k+1; end x1=0.5*x0+h*F; if abs(x1-x0)ep break; end m=m+1; h=h/2; x0=x1; i=i+1;endx1i2、复化辛卜生公式：（1）算法步骤： 将区间分成等分，为偶数，步长为，分点为.其中，我们又题意已知：。 在每个小区间上用辛卜生公式，即。再将各个小区间的积分相加，即。 给定误差，当时，停止计算，取为所给积分的近似值.（2）程序代码：clca=1;b=2;n=2;h=(b-a)/n;f1=exp(1);f2=2*exp(2);f= 6.7225;S1=(b-a)/6*f1+4*f+f2;i=0;ep=0.5e-7;while i

28、100; P=(a+h)*exp(a+h); Q=0; j=2; while jn P=P+(a+(j+1)*h)*exp(a+(j+1)*h); Q=Q+(a+j*h)*exp(a+j*h); j=j+2; end S=h/3*f1+4*P+2*Q+f2; if(abs(S1-S)ep) break; end h=h/2; n=2*n; S1=S; i=i+1;endSi3、 运行结果：（1）复化梯形法： x1 = 7.3891；i= 13（2）复化辛卜生法：S =7.3891；i=64、结果分析：从结果上我们可以看出，为了达到相同的精度，用复化梯形公式时，需将1,2分成；用复化辛卜生公式时，需将1,2分成等分。说明了复化辛卜生公式的计算量比复化梯形公式的计算量要小。