计算方法实验报告doc.docx

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

2、迭代公式得. 如果，则迭代结束，取为解的近似值；否则，用代替，重复步骤和步骤. （2）程序代码 在区间内， 代码 clc x0-3.5; 初值 iter_max100; 迭代的最大次数 ep1e-6; 允许精度 k0; while kiter_max k从0开始到iter_max循环 x1sinx012*x0-1.1/3; 代入，算出的值 if absx1-x0ep 与允许精度作比较 break; 条件成立，跳出循环 end x0x1; 条件不成立，用代替 kk1; k加1 end x_starx1, iterk 为解的近似值，iter为迭代次数 运行结果x_star -3.4101 ；ite

3、r 14 在区间内， 代码 clc x00.5; 初值 iter_max100; 迭代的最大次数 ep1e-6; 允许精度 k0; while kiter_max k从0开始到iter_max循环 x11/12*x0.3-sinx01; 代入，算出的值 if absx1-x0ep 与允许精度作比较 break; 条件成立，跳出循环 end x0x1; 条件不成立，用代替 kk1; k加1 end x_starx1, iterk 为解的近似值，iter为迭代次数 结果x_star 0.07696；iter 6 在区间内， 代码 clc x03.5; 初值 iter_max100; 迭代的最大次数

4、 ep1e-6; 允许精度 k0; while kiter_max k从0开始到iter_max循环 x1sinx012*x0-1.1/3; 代入，算出的值 if absx1-x0ep 与允许精度作比较 break; 条件成立，跳出循环 end x0x1; 条件不成立，用代替 kk1; k加1 end x_starx1, iterk 为解的近似值，iter为迭代次数 运行结果x_star 3.4101 ；iter 10 2、 牛顿迭代法 （1）算法步骤 选定初值，计算，. 按公式迭代，得新的近似值，并计算，. 对于给定的允许精度，如果，则终止迭代，取；否则，在转回步骤计算. （2）程序代码 在

5、区间内， clc x1-3.5; 初值 k0; while k100 k从0开始到100循环 x0x1; 将初值赋给 f0x0.3-sinx0-12*x01; 计算 f13*x0.2-cosx0-12; 计算 x1x0-f0/f1; 计算得到新的近似值 if absx1-x0 1.0e-6 与允许精度作比较 break; 条件成立，跳出循环 end kk1; 条件不成立，k加1 end x_starx1, iterk 为解的近似值，iter为迭代次数 运行结果x_star -3.4911；iter 2 在区间内， clc x10.5; 初值 k0; while k100 k从0开始到100循环

6、 x0x1; 将初值赋给 f0x0.3-sinx0-12*x01; 计算 f13*x0.2-cosx0-12; 计算 x1x0-f0/f1; 计算得到新的近似值 if absx1-x0 1.0e-6 与允许精度作比较 break; 条件成立，跳出循环 end kk1; 条件不成立，k加1 end x_starx1, iterk 为解的近似值，iter为迭代次数 运行结果x_star 0.07696 ；iter 3 在区间内， clc x13.5; 初值 k0; while k100 k从0开始到100循环 x0x1; 将初值赋给 f0x0.3-sinx0-12*x01; 计算 f13*x0.2

7、-cosx0-12; 计算 x1x0-f0/f1; 计算得到新的近似值 if absx1-x0 1.0e-6 与允许精度作比较 break; 条件成立，跳出循环 end kk1; 条件不成立，k加1 end x_starx1, iterk 为解的近似值，iter为迭代次数 运行结果x_star 3.4101；iter 3 3、运行结果 x_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 3 4、结果分析 从这题的结果可以看出，牛顿迭代法的迭代速度比

8、一般迭代法的速度要快，牛顿法的迭代次数比较少。 例如，求在的根时，一般迭代法迭代了14次（iter 14），而牛顿法只用了2次（iter 2）。 总而言之，一般迭代法是一种逐次逼近的方法，具有原理简单、编制程序方便等优点，但存在是否收敛和收敛速度的快慢等问题。 牛顿迭代法是一种特殊的迭代法，用于求方程的单根时具有2阶收敛。 因此是一种求非线性方程解的好方法，还可以用于求重根和复根，而且可以推广到求非线性方程组的解. 二、解方程组直接法 1、已知 对矩阵A做LU分解。 2、用追赶法解下述方程组，并给出n10的结果，其中 ， 1、杜利特尔分解法（直接三角分解法） 设方程组的系数矩阵的各阶主子式，则

9、可知该方程组存在唯一的杜利特尔分解，其中 ，. 我们记矩阵的第行第列元素为，则矩阵的第一行和的第一列可由公式（1.1）求出 （2.1） 再由公式（2.2）求出的第行、的第列元素 （2.2） （1）算法步骤 由于系数矩阵的各阶主子式，则可知该方程组存在唯一的杜利特尔分解。 利用公式（2.1）和（2.2）求出和的元素和。 求解单位下三角形方程组，即按公式 （2.3） 求出。 求解上三角形方程组，即按公式 （2.4） 可求得方程组的解。 （2）程序代码 function L,ULUA A n,nsizeA; 定义为阶矩阵 Lzerosn,n; 矩阵初始化 Uzerosn,n; 矩阵初始化 for i

10、1n Li,i1; 矩阵的对角线元素为1 end for k1n for jkn Uk,jAk,j-sumLk,1k-1.*U1k-1,j ; 求出元素 end for ik1n Li,kAi,k-sumLi,1k-1.*U1k-1,k /Uk,k; 求出的元素 end end A4 2 1 5;8 7 2 10;4 8 3 6;12 6 11 20; 矩阵 L,ULUA 2、追赶法 设所求方程组为 ， （2.5） （1）算法步骤 由方程组的第1个方程，将未知元用表示为 记，得 。 以此类推，用表示， 用表示。 我们可以得到一般式（其中记） （2.6） （2.7） 将代入方程组（2.5）的第n

11、个方程，并解出得 ， 上式右端恰好是式（2.6）的。 于是，可由式（2.7）求出三对角方程组的解，即 ， （2.8） （2）程序代码 clc a0 1 1 1 1 1 1 1 1 1; b2 2 2 2 2 2 2 2 2 2; c1 1 1 1 1 1 1 1 1 0; r-7 -5 -5 -5 -5 -5 -5 -5 -5 -5; u0 0 0 0 0 0 0 0 0 0; v0 0 0 0 0 0 0 0 0 0; x0 0 0 0 0 0 0 0 0 0; u1r1/b1; 计算 v1c1/b1; 计算 for k210 ukrk-uk-1*ak/bk-vk-1*ak; 计算 vkck

12、/bk-vk-1*ak; 计算 end x10u10; for k19 xkuk-vk*xk1; 由，解出 end x 运行结果 x-3.5000 -1.0000 -3.0000 -1.6000 -2.8333 -1.8571 -2.7500 -2.0000 -0.8182 -2.0909 3、运行结果 （1）杜利特尔分解法 （2）追赶法 。 4、结果分析 （1）杜利特尔分解法在实现分解后，解具有相同系数矩阵的方程组，非常方便，每解一个方程只需用式（2.3）和式（2.4）解两个三角形方程组即可，大大减少了运算工作量。 另外，矩阵和的元素可采用紧凑格式存放在系数矩阵的相应元素位置上，节省了存储单

13、元。 （2）追赶法追赶法仅适用于三对角方程组的求解。 我们由式（2.6）和式（2.8）可以算出，追赶法的乘、除运算次数仅为，加、减运算次数仅为。 故追赶法具有运算量小和存储量小的优势。 三、解方程组迭代法 1、用迭代法解Axb，其中b5，5，5T， 给定误差，用Jacobi和SOR两种迭代法计算，并给出n10的结果。 对于阶方程组，假定系数矩阵的对角元时，可得雅可比迭代格式为 ， （3.1） 记为方程组的系数矩阵的对角元组成的对角阵，和分别为由的元素组成的严格下三角阵和上三角阵，则系数矩阵 可分解为 . 1、Jacobi（雅可比）迭代法 （1）算法步骤 根据矩阵，算出矩阵、和，即 ， 。 计算

14、雅可比迭代矩阵 . 计算矩阵 . 代入给出的初值，得到雅可比迭代公式（3.1）的矩阵形式 . 3.2 由式（3.2）就可以得到向量序列 . （2）程序代码 clc A3 -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

15、;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; x00 0 0 0 0 0 0 0 0 0 ; 初值 b5 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

16、 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

17、;0 0 0 0 0 0 0 0 0 0; D3 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; BJinvD*LU; 计算雅可比迭代矩阵 FJinvD*b; 计算矩阵 N1000; ep1e-10; 误差 k0; while ki

18、ter_max k从0开始到iter_max循环 x1BJ*x0fJ; 计算 if normx1-x0, inf ep 与误差做比较 break; 满足，跳出循环 end x0x1; 不满足，将赋给 kk1; k加1 end x_starx1, iterk 2、SOR（超松弛）迭代法 （1）算法步骤 根据矩阵，算出矩阵、和，即 ， 。 计算超松弛迭代矩阵 . 其中，为松弛因子，。 计算矩阵 . 其中，为松弛因子，。 代入给出的初值，得到超松弛迭代公式的矩阵形式 . 由式（3.2）就可以得到向量序列 . （2）程序代码 clc A3 -1/2 -1/4 0 0 0 0 0 0 0;-1/2 3

19、-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; x00 0 0 0 0 0 0 0 0 0 ; b5

20、 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

21、 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; D3 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

22、 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; w1.3; 超松弛因子 BwinvD-w*L*1-w*Dw*U; 计算超松弛迭代矩阵 Fww*invD-w*L*b; 计算矩阵 iter_max1000; ep1e-10; 误差 k0; while kiter_max x1Bw*x0fw; if normx1-x0, inf ep 与误差做比较 break; end x0x1; kk1; end x_star

23、x1, iterk 3、运行结果 （1）雅可比迭代法 迭代次数iter 31 （2）超松弛迭代法 迭代次数iter 25 4、结果分析 雅可比迭代法和超松弛迭代法从迭代次数上比较，可以看出超松弛迭代比雅可比迭代的次数少，表明超松弛迭代法比雅可比迭代法的收敛速度快。 四、插值逼近 题目 ，将10等分，作Lagrange插值，将插值函数的图形与的图形比较，并给出结论。 1、算法步骤 已知， 求出各个节点的拉格朗日插值基函数多项式 ， 就称为以为节点的拉格朗日插值基多项式。 利用基函数的线性组合得到插值多项式 。 2、程序代码 function ylagrangex0,y0,x nlengthx0;

24、 mlengthx; for i1m zxi; s0.0; for k1n p1.0; for j1n if jk pp*z-x0j/x0k-x0j; 求出，p即为 end end sp*y0ks; 计算，s即为 end yis; end x-515; y1./1x.2; 插值函数 x0-50.15; y0lagrangex,y,x0; y11./1x0.2; 绘制图形 plotx0,y0, -r 插值曲线 hold on plotx0,y1, -b 原曲线 3、 运行结果 插值曲线为红色，原曲线为蓝色 4、结果分析 当很大时，在内能很好的逼近，而在这个区间外，误差反而很大。 由图像我们可以看

25、出，在内部，与偏差小；而在端点附近，偏差大。 这就是龙格现象。 五、复化求积 题目 分别用复化梯形公式、复化辛卜生公式计算 , 其中. 用区间逐步分半递推算法 1、复化梯形公式 （1）算法步骤 将区间分成等分，步长为 ， 分点为 . 其中，我们又题意已知。 建立与的递推公式 。 给定误差，当时，停止计算，取为所给积分的近似值。 （2）程序代码 clc a1; b2; m1; h0.5; ep0.5e-7; faexp1; fb2*exp2; x0h*fafb; iter_max100; i0; while iiter_max k1; F0; while k2m-1 FFa2*k-1*h*exp

26、a2*k-1*h; kk1; end x10.5*x0h*F; if absx1-x0ep break; end mm1; hh/2; x0x1; ii1; end x1 i 2、复化辛卜生公式 （1）算法步骤 将区间分成等分，为偶数，步长为 ， 分点为 . 其中，我们又题意已知。 在每个小区间上用辛卜生公式，即 。 再将各个小区间的积分相加，即 。 给定误差，当时，停止计算，取为所给积分的近似值. （2）程序代码 clc a1; b2; n2; hb-a/n; f1exp1; f22*exp2; f 6.7225; S1b-a/6*f14*ff2; i0; ep0.5e-7; while i100; Pah*expah; Q0; j2; while jn PPaj1*h*expaj1*h; QQaj*h*expaj*h; jj2; end Sh/3*f14*P2*Qf2; ifabsS1-Sep break; end hh/2; n2*n; S1S; ii1; end S i 3、 运行结果 （1）复化梯形法 x1 7.3891；i 13 （2）复化辛卜生法S 7.3891；i6 4、结果分析 从结果上我们可以看出，为了达到相同的精度，用复化梯形公式时，需将1,2分成；用复化辛卜生公式时，需将1,2分成等分。 说明了复化辛卜生公式的计算量比复化梯形公式的计算量要小。