计算方法实验报告

实验报告一、求方程 f(x)=x3-sinx-12x+1 的全部根, =1e-61、 用一般迭代法; 2、 用牛顿迭代法; 并比较两种迭代的收敛速度。一、首先，由题可求得： .12cos3)(2 xxf其次，分析得到其根所在的区间。 令 ，可得到 .0xf sin13 用一阶导数分析得到 和 两个函数的增减区间；再用二阶导数分析得2x到两个函数的拐点以及凹凸区间. 在直角坐标轴上描摹出 和 的图，在图上可以看到他们的交点，013sinx然后估计交点所在的区间，即是所要求的根的区间。经过估计，得到根所在的区间为， 和 .3,41,041、一般迭代法 （1）算法步骤：设 为给定的允许精度，迭代法的计算步骤为： 选定初值 .由 确定函数 ，得等价形式 .0xfxgxg 计算 .由迭代公式得 .g01 如果 ，则迭代结束，取 为解的近似值；否则，用 代替 ，重01x1x1x0复步骤和步骤.（2）程序代码： 在区间 内，3,4代码： clcx0=-3.5; %初值 0xiter_max=100; %迭代的最大次数ep=1e-6; %允许精度 k=0;while k=iter_max %k从0 开始到iter_max循环 x1=(sin(x0)+12*x0-1).(1/3); %代入 ，算出 的值0x1if abs(x1-x0)ep % 与允许精度作比较01xbreak; %条件 成立，跳出循环01xendx0=x1; %条件 不成立，用 代替011x0k=k+1; %k加1endx_star=x1, iter=k % 为解的近似值，iter为迭代次数1x运行结果：x_star = -3.4101 ；iter =14在区间 内， 1,0代码：clcx0=0.5; %初值 0xiter_max=100; %迭代的最大次数ep=1e-6; %允许精度 k=0;while k=iter_max %k从0 开始到iter_max循环x1=(1/12)*(x0.3-sin(x0)+1); %代入 ，算出 的值0x1if abs(x1-x0)ep % 与允许精度作比较01xbreak; %条件 成立，跳出循环01xendx0=x1; %条件 不成立，用 代替01x1x0k=k+1; %k加1endx_star=x1, iter=k % 为解的近似值，iter为迭代次数1x结果：x_star = 0.07696；iter =6在区间 内，4,3代码：clcx0=3.5; %初值 0xiter_max=100; %迭代的最大次数ep=1e-6; %允许精度 k=0;while k=iter_max %k从0 开始到iter_max循环x1=(sin(x0)+12*x0-1).(1/3); %代入 ，算出 的值0x1if abs(x1-x0)ep % 与允许精度作比较01xbreak; %条件 成立，跳出循环01xendx0=x1; %条件 不成立，用 代替 011x0k=k+1; %k加1endx_star=x1, iter=k % 为解的近似值，iter为迭代次数1x运行结果：x_star = 3.4101 ；iter =102、 牛顿迭代法（1）算法步骤： 选定初值 ，计算 ， .0x0xf0f 按公式 迭代，得新的近似值 ，并计算 ， .kkf1 1kx1kxf1kf 对于给定的允许精度 ，如果 ，则终止迭代，取 ；否则，kx1 1*k，在转回步骤计算.1k（2）程序代码：在区间 内，3,4clcx1=-3.5; %初值 1xk=0;while k=100 %k从 0开始到100循环x0=x1; %将初值 赋给 1xf0=x0.3-sin(x0)-12*x0+1; %计算 0xff1=3*x0.2-cos(x0)-12; %计算 x1=x0-f0/f1; %计算 得到新的近似值 1xif abs(x1-x0) 1.0e-6 % 与允许精度作比较0break; %条件 成立，跳出循环01xendk=k+1; %条件 不成立，k加101end x_star=x1, iter=k % 为解的近似值，iter为迭代次数1x运行结果：x_star = -3.4911；iter =2在区间 内，1,0clcx1=0.5; %初值 1xk=0;while k=100 %k从 0开始到100循环x0=x1; %将初值 赋给1xf0=x0.3-sin(x0)-12*x0+1; %计算 0xff1=3*x0.2-cos(x0)-12; %计算 x1=x0-f0/f1; %计算 得到新的近似值 1xif abs(x1-x0) 1.0e-6 % 与允许精度作比较0break; %条件 成立，跳出循环01xendk=k+1; %条件 不成立，k加101end x_star=x1, iter=k % 为解的近似值，iter为迭代次数1x运行结果：x_star =0.07696 ；iter =3在区间 内，4,3clcx1=3.5; %初值 1xk=0; while k=100 %k从 0开始到100循环x0=x1; %将初值 赋给1xf0=x0.3-sin(x0)-12*x0+1; %计算 0xff1=3*x0.2-cos(x0)-12; %计算 x1=x0-f0/f1; %计算 得到新的近似值 1xif abs(x1-x0) 1.0e-6 % 与允许精度作比较01xbreak; %条件 成立，跳出循环01xendk=k+1; %条件 不成立，k加101end x_star=x1, iter=k % 为解的近似值，iter为迭代次数1x运行结果：x_star =3.4101；iter =33、运行结果： 3,4x1,0x4,3xx_star iter x_star iter x_star iter一般迭代 3.4101 14 0.7696 6 3.4101 10牛顿法 3.4911 2 0.70696 3 3.4101 34、结果分析： 从这题的结果可以看出，牛顿迭代法的迭代速度比一般迭代法的速度要快，牛顿法的迭代次数比较少。例如，求在 的根时，一般迭代法迭代了 14 次（iter =14） ，而牛3,4顿法只用了 2 次（iter =2 ） 。总而言之，一般迭代法是一种逐次逼近的方法，具有原理简单、编制程序方便等优点，但存在是否收敛和收敛速度的快慢等问题。牛顿迭代法是一种特殊的迭代法，用于求方程的单根时具有 2 阶收敛。因此是一种求非线性方程解的好方法，还可以用于求重根和复根，而且可以推广到求非线性方程组的解.二、解方程组直接法1、已知421587036A对矩阵 A 做 LU 分解。2、用追赶法解下述方程组 ，并给出 n=10 的结果，其中xb2112A， 57b1、杜利特尔分解法（直接三角分解法）：设方程组 的系数矩阵 的各阶主子式 ，则可bAxAnkA,210det知该方程组存在唯一的杜利特尔分解： ，其中LU， .1121nnlllL nnuu 2211我们记矩阵 的第 行第 列元素为 ，则矩阵 的第一行和 的第一列可由AijijaUL公式（1.1）求出：（2.1）.,32,11niualjij 再由公式（2.2）求出 的第 行 、 的第 列 元素：Uk, Lk1,2n（2.2）.,2,1,1 nkiulaljkkpiiikpjkj （1）算法步骤：由于系数矩阵 的各阶主子式 ，则可知该方程组存AnkA,210det在唯一的杜利特尔分解： 。LU 利用公式（2.1）和（2.2）求出 和 的元素 和 。kiuil 求解单位下三角形方程组 ，即按公式by（2.3）1,32,kjjk nly求出 。n,21 求解上三角形方程组 ，即按公式yUx（2.4）1,2,1nkuxyxknjjkn可求得方程组的解： 。1,xn（2）程序代码：function L,U=LU(A)An,n=size(A); %定义 为 阶矩阵AnL=zeros(n,n); %矩阵 初始化LU=zeros(n,n); %矩阵 初始化Ufor i=1:nL(i,i)=1; %矩阵 的对角线元素为1endfor k=1:nfor j=k:nU(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j);%求出 元素 Ukiuendfor i=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)/U(k,k); %求出 的元素LiklendendA=4 2 1 5;8 7 2 10;4 8 3 6;12 6 11 20; %矩阵 AL,U=LU(A) % LUA2、追赶法设所求方程组为， （2.5）nnnrxbacxrxbac11121343322121（1）算法步骤： 由方程组的第 1 个方程，将未知元 用 表示为1x21bcr记 ， ，得1bru1cv。211xvu 以此类推， 用 表示， 用 表示。我们可以得到一般式（其中记2x3k）：0,1nca（2.6）1,2,11nkavbcurkkk （2.7）,xux 将 代入方程组（2.5）的第 n 个方程，并解出 得nnv11 nx，1nnvaburx上式右端恰好是式（2.6）的 。于是，可由式（2.7）求出三对角方程组的n解 ，即kx， （2.8）1,2,1nkxvuukkn（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 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); %计算 1bruv(1)=c(1)/b(1); %计算 1cvfor k=2:10u(k)=(r(k)-u(k-1)*a(k)/(b(k)-v(k-1)*a(k); %计算 kuv(k)=c(k)/(b(k)-v(k-1)*a(k); %计算 kvendx(10)=u(10); for k=1:9x(k)=u(k)-v(k)*x(k+1); %由 ，解出1kkxvux1kxend x运行结果：x=-3.5000 -1.0000 -3.0000 -1.6000 -2.8333 -1.8571 -2.7500 -2.0000 -0.8182 -2.09093、运行结果：（1）杜利特尔分解法： 12035414032UL（2）追赶法： 。09.281.75.83260.1.5319654321xxx4、结果分析：（1）杜利特尔分解法：在实现 分解后，解具有相同系数矩阵的方程组LUA，非常方便，每解一个方程只需用式（2.3）和式mibAx,21（2.4）解两个三角形方程组即可，大大减少了运算工作量。另外，矩阵和 的元素可采用紧凑格式存放在系数矩阵 的相应元素位置上，节省ULA了存储单元。（2）追赶法：追赶法仅适用于三对角方程组的求解。我们由式（2.6）和式（2.8）可以算出，追赶法的乘、除运算次数仅为 ，加、减运算次数45n仅为 。故追赶法具有运算量小和存储量小的优势。13n三、解方程组迭代法1、用迭代法解 Ax=b，其中 b=(5，5，5) T，132411342241134241342A 给定误差 ，用 Jacobi 和 SOR 两种迭代法计算，并给出10n=10 的结果。对于 阶方程组 ，假定系数矩阵 的对角元 时，nbAxAniai ,210可得雅可比迭代格式为：， （3.1） ,210;,21,0201mni xaxbaxijnijmjjmi T给 出 初 值记 为方程组 的系数矩阵 的对角元组成的对角阵，ndiagD,21bAA和 分别为由 的元素组成的严格下三角阵和上三角阵，则系数矩阵LUA A可分解为：.ULD1、Jacobi（雅可比）迭代法：（1）算法步骤： 根据矩阵 ，算出矩阵 、 和 ，即ALUD， ，00121321nnaaL 00122311nnaaU。nadigD,2 计算雅可比迭代矩阵 ：JB.ULJ1 计算矩阵 ：Jf.bDfJ1 代入给出的初值 ，得到雅可比迭代公式（3.1）的矩Tnxx0201,阵形式：. (3.2),1,1mfBJJm由式（3.2）就可以得到向量序列：. ,10xx（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 -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; %初值 Tnxx021,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 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 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 0 0 0 0 0 0 3;BJ=inv(D)*(L+U); %计算雅可比迭代矩阵 ULDBJ1FJ=inv(D)*b; %计算矩阵 bfJ1N=1000;ep=1e-10; % 误差 10k=0;while k=iter_max %k从0 开始到iter_max循环x1=BJ*x0+fJ; %计算 JmJfxB1if norm(x1-x0),inf)ep % 与误差 做比较01break; %满足 ，跳出循环01xendx0=x1; %不满足 ，将 赋给011xk=k+1; %k加1endx_star=x1, iter=k 2、SOR（超松弛）迭代法：（1）算法步骤： 根据矩阵 ，算出矩阵 、 和 ，即ALUD， ，00121321nnaaL 00122311nnaaU。nadigD,2 计算超松弛迭代矩阵 ：B.UDL1其中， 为松弛因子， 。, 计算矩阵 ：f.bLDfJ1其中， 为松弛因子， 。, 代入给出的初值 ，得到超松弛迭代公式的矩阵形式：Tnxx0201. ,1,mfBm由式（3.2）就可以得到向量序列：. ,10xx（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 -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 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 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 0 0 0 0 0 0 3;w=1.3; %超松弛因子 Bw=(inv(D-w*L)*(1-w)*D+w*U); %计算超松弛迭代矩阵 BFw=w*(inv(D-w*L)*b; %计算矩阵 fiter_max=1000;ep=1e-10; %误差 10k=0;while k=iter_maxx1=Bw*x0+fw; % fxBm1if norm(x1-x0),inf)ep % 与误差 做比较01xbreak;endx0=x1; k=k+1;endx_star=x1, iter=k3、运行结果：（1）雅可比迭代法： T2.4079) .8633.1.253.9 .2 3.58 .162 .8 2.4079(x迭代次数：iter =31（2）超松弛迭代法： T1 2.4079) .863.13.25.9 3.2 .58 3.162 .8 .4079(x迭代次数：iter =254、结果分析：雅可比迭代法和超松弛迭代法从迭代次数上比较，可以看出超松弛迭代比雅可比迭代的次数少，表明超松弛迭代法比雅可比迭代法的收敛速度快。四、插值逼近题目: ，将 10 等分，作 Lagrange 插值，5,15)(2xxf 5,将插值函数的图形与 的图形比较，并给出结论。)(fy1,()nnjkjkjxxy1、算法步骤：已知 ，xfy 求出各个节点的拉格朗日插值基函数多项式：， nikkiiiiiii ni xxxxxl 01110 就称为以 为节点的拉格朗日插值基多项式。li i 利用基函数 的线性组合得到插值多项式：xli。xlyxlyLnn02、程序代码：function y=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:nif j=kp=p*(z-x0(j)/(x0(k)-x0(j); %求出 ，p即为nikkixl0ilendends=p*y0(k)+s; %计算 ，s即为xlyxlyLnn0 xLnendy(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、结果分析：当 很大时，在 内 能很好的逼近 ，而在这个区间外，误n63.xxLnxf差反而很大。由图像我们可以看出，在 内部， 与 偏差小；而在5,L10f端点 附近，偏差大。这就是龙格现象。5x五、复化求积题目: 分别用复化梯形公式、复化辛卜生公式计算,21dxe其中 . (用区间逐步分半递推算法) 7105.| kx1、复化梯形公式：（1）算法步骤： 将区间 分成 等分，步长为ba,210