计算方法实验报告

实验报告一、求方程f(x)=x^3-sinx-12x+1的全部根, ε=1e -6 1、 用一般迭代法; 2、 用牛顿迭代法;并比较两种迭代的收敛速度。

一、首先，由题可求得：12cos 3)(2'--=x x x f . 其次，分析得到其根所在的区间。

① 令()0=x f ，可得到x x x sin 1123=+-.② 用一阶导数分析得到1123+-x x 和x sin 两个函数的增减区间；再用二阶导数分析得到两个函数的拐点以及凹凸区间.③ 在直角坐标轴上描摹出01123=+-x x 和0sin =x 的图，在图上可以看到他们的交点，然后估计交点所在的区间，即是所要求的根的区间。

经过估计，得到根所在的区间为[]3,4--，[]1,0和[]4,3.1、 一般迭代法 （1）算法步骤：设ε为给定的允许精度，迭代法的计算步骤为：① 选定初值0x .由()0=x f 确定函数()x g ，得等价形式()x g x =. ② 计算()0x g .由迭代公式得()01x g x =.③ 如果ε≤-01x x ，则迭代结束，取1x 为解的近似值；否则，用1x 代替0x ，重复步骤②和步骤③. （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); %代入0x ，算出1x 的值 if abs(x1-x0)<ep %01x x -与允许精度作比较 break ; %条件ε≤-01x x 成立，跳出循环 endx0=x1; %条件ε≤-01x x 不成立，用1x 代替0x k=k+1; %k 加1 endx_star=x1, iter=k %1x 为解的近似值，iter 为迭代次数 运行结果：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); %代入0x ，算出1x 的值if abs(x1-x0)<ep %01x x -与允许精度作比较 break ; %条件ε≤-01x x 成立，跳出循环 endx0=x1; %条件ε≤-01x x 不成立，用1x 代替0x k=k+1; %k 加1 endx_star=x1, iter=k %1x 为解的近似值，iter 为迭代次数 结果：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); %代入0x ，算出1x 的值 if abs(x1-x0)<ep %01x x -与允许精度作比较 break ; %条件ε≤-01x x 成立，跳出循环 endx0=x1; %条件ε≤-01x x 不成立，用1x 代替0x k=k+1; %k 加1 endx_star=x1, iter=k %1x 为解的近似值，iter 为迭代次数 运行结果：x_star = 3.4101 ；iter =102、 牛顿迭代法 （1）算法步骤：① 选定初值0x ，计算()0x f ，()0'x f .② 按公式()()k k k k x f x f x x '1-=+迭代，得新的近似值1+k x ，并计算()1+k x f ，()1'+k x f . ③ 对于给定的允许精度ε，如果ε≤-+k k x x 1，则终止迭代，取1*+≈k x x ；否则，1+=k k ，在转回步骤②计算.（2）程序代码： ①在区间[]3,4--内， clcx1=-3.5; %初值1x k=0;while k<=100 %k 从0开始到100循环 x0=x1; %将初值1x 赋给0xf0=x0.^3-sin(x0)-12*x0+1; %计算()0x f f1=3*x0.^2-cos(x0)-12; %计算()0'x f x1=x0-f0/f1; %计算得到新的近似值1xif abs(x1-x0)< 1.0e-6 %01x x -与允许精度作比较 break ; %条件ε≤-01x x 成立，跳出循环 endk=k+1; %条件ε≤-01x x 不成立，k 加1endx_star=x1, iter=k %1x 为解的近似值，iter 为迭代次数 运行结果：x_star = -3.4911；iter =2②在区间[]1,0内， clcx1=0.5; %初值1x k=0;while k<=100 %k 从0开始到100循环 x0=x1; %将初值1x 赋给0xf0=x0.^3-sin(x0)-12*x0+1; %计算()0x f f1=3*x0.^2-cos(x0)-12; %计算()0'x f x1=x0-f0/f1; %计算得到新的近似值1xif abs(x1-x0)< 1.0e-6 %01x x -与允许精度作比较 break ; %条件ε≤-01x x 成立，跳出循环 endk=k+1; %条件ε≤-01x x 不成立，k 加1 endx_star=x1, iter=k %1x 为解的近似值，iter 为迭代次数 运行结果：x_star =0.07696 ；iter =3③在区间[]4,3内， clcx1=3.5; %初值1x k=0;while k<=100 %k 从0开始到100循环 x0=x1; %将初值1x 赋给0xf0=x0.^3-sin(x0)-12*x0+1; %计算()0x f f1=3*x0.^2-cos(x0)-12; %计算()0'x f x1=x0-f0/f1; %计算得到新的近似值1xif abs(x1-x0)< 1.0e-6 %01x x -与允许精度作比较 break ; %条件ε≤-01x x 成立，跳出循环 endk=k+1; %条件ε≤-01x x 不成立，k 加1 endx_star=x1, iter=k %1x 为解的近似值，iter 为迭代次数 运行结果：x_star =3.4101；iter =33、运行结果：4、结果分析：从这题的结果可以看出，牛顿迭代法的迭代速度比一般迭代法的速度要快，牛顿法的迭代次数比较少。

例如，求在[]3,4--的根时，一般迭代法迭代了14次（iter =14），而牛顿法只用了2次（iter =2）。

总而言之，一般迭代法是一种逐次逼近的方法，具有原理简单、编制程序方便等优点，但存在是否收敛和收敛速度的快慢等问题。

牛顿迭代法是一种特殊的迭代法，用于求方程的单根时具有2阶收敛。

因此是一种求非线性方程解的好方法，还可以用于求重根和复根，而且可以推广到求非线性方程组的解.二、解方程组直接法1、已知42158721048361261120A ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦ 对矩阵A 做LU 分解。

2、用追赶法解下述方程组Ax b =，并给出n=10的结果，其中⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=2112112112 A ，⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡----=5557 b1、杜利特尔分解法（直接三角分解法）：设方程组b Ax =的系数矩阵A 的各阶主子式()()n k A k ,,2,10det =≠，则可知该方程组存在唯一的杜利特尔分解：LU A =，其中⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=-11112121nn n n l l l l L ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=nn n n u u u u u u U 22211211. 我们记矩阵A 的第i 行第j 列元素为ij a ，则矩阵U 的第一行和L 的第一列可由公式（1.1）求出：⎪⎩⎪⎨⎧====.,,3,2,,,,2,1,111111n i ua l n j a u i i j j （2.1） 再由公式（2.2）求出U 的第k 行()n k ,,2 =、L 的第k 列()1,,2-=n k 元素：⎪⎪⎪⎩⎪⎪⎪⎨⎧++=-=+=-=∑∑-=-=.,,2,1,,,,1,,1111n k k i u u l a l n k k j u l a u kk k p pkip ik ik k p pj kp kj kj （2.2）（1）算法步骤：由于系数矩阵A 的各阶主子式()()n k A k ,,2,10det =≠，则可知该方程组存在唯一的杜利特尔分解：LU A =。

① 利用公式（2.1）和（2.2）求出U 和L 的元素ki u 和ik l 。

② 求解单位下三角形方程组b Ly =，即按公式⎪⎩⎪⎨⎧=-==∑-=1111,,3,2,,k j jkj k k n k y l b y b y （2.3）求出n y y y ,,,21 。

④ 求解上三角形方程组y Ux =，即按公式⎪⎪⎪⎩⎪⎪⎪⎨⎧--=-==∑+=1,,2,1,,1n n k u x u y x u y x kk nk j jkj k k nnn n （2.4） 可求得方程组的解：11,,,x x x n n -。

（2）程序代码：function [L,U]=LU(A)A[n,n]=size(A); %定义A 为n n ⨯阶矩阵 L=zeros(n,n); %矩阵L 初始化 U=zeros(n,n); %矩阵U 初始化 for i=1:nL(i,i)=1; %矩阵L 的对角线元素为1endfor k=1:n for j=k:nU(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j)');%求出U 元素ki uendfor i=k+1:nL(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)'))/U(k,k);%求出L 的元素ik lend endA=[4 2 1 5;8 7 2 10;4 8 3 6;12 6 11 20]; %矩阵A [L,U]=LU(A) %LU A =2、追赶法设所求方程组为⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧=+=++=++=++=+-------n n n n n n n n n n n n r x b x a r x c x b x a r x c x b x a r x c x b x a r x c x b 11111213433323232221212111 ， （2.5）（1）算法步骤：① 由方程组的第1个方程，将未知元1x 用2x 表示为12111b x c r x -=记111b r u =，111b cv =，得 2111x v u x -=。