《计算方法》上机实验报告班级：XXXXXX小组成员：XXXXXXXXXXXXXXXXXXXXXXXXXXXX任课教师：XXX二〇一八年五月二十五日前言通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB 程序的编写。

以下为本次上机实验报告，按照实验内容共分为六部分。

实验一：一、实验名称及题目： Newton 迭代法例2.7(P38):应用Newton 迭代法求x 3−x −1=0在x =1附近的数值解x k ，并使其满足|x k −x k−1|<10−8. 二、解题思路：设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交点的横坐标)(')(0001x f x f x x -=,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标)(')(1112x f x f x x -=称2x 为'x的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把)(')(1n n n n x f x f x x -=+称为'x 的1+n 次近似值，这种求解方法就是牛顿迭代法。

三、Matlab 程序代码：function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1;f1=diff(f);%求导 y=subs(f,z,x0);y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1;while abs(x1-x0)>=tol x0=x1;y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; endx=double(x1) K四、运行结果：实验二：一、实验名称及题目：Jacobi 迭代法例3.7（P74）：试利用Jacobi 迭代公式求解方程组[5−1−1−1−110−1−1−1−15−1−1−1−110][x 1x 2x 3x 4]=[−412834]要求数值解X (k)满足||X −X (k )||2≤10−4，其中X =(1,2,3,4)T 为方程组的精确解. 二、解题思路：首先将方程组中的系数矩阵A 分解成三部分，即：U D L A ++=，D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。

之后确定迭代格式，f X B X k k +=+)()1(*，（ ⋅⋅⋅=2,1,0k , k 即迭代次数），B 称为迭代矩阵。

最后选取初始迭代向量)0(X ，开始逐次迭代。

最后验证精度。

（迭代阵：b D UXD Xk k1)(1)1(--++-=。

）雅克比迭代法的优点明显，计算公式简单，每迭代一次只需计算一次矩阵和向量的乘法，且计算过程中原始矩阵A 始终不变，比较容易并行计算。

然而这种迭代方式收敛速度较慢，而且占据的存储空间较大。

三、Matlab 程序代码：function jacobi(A,b,x0,eps,x1) D = diag(diag(A));%求A 的对角矩阵 L = -tril(A,-1);%求A 的下三角矩阵U = -triu(A,1);%求A的上三角矩阵B = D\(L+U);f = D\b;x = B*x0+f;n = 1;%迭代次数while norm(x-x1)>=epsx = B*x+f;n = n+1;endformat longnxjingdu=norm(x-x1)四、运行结果：实验三：一、实验名称及题目：Gauss-Seidel 迭代法例 3.8（P75）:试利用Gauss-Seidel迭代公式求解方程组[5−1−1−1−110−1−1−1−15−1−1−1−110][x1x2x3x4]=[−412834]，并使其数值解X(k)满足精度要求||X−X(k)||2≤10−4，其中X=(1,2,3,4)T为方程组的精确解.二、解题思路：Gauss-Seidel 迭代法与Jacobi 迭代法思路相近，首先将方程组中的系数矩阵A 分解成三部分，即：U D L A ++=，D 为对角阵，L 为下三角矩阵，U 为上三角矩阵。

之后确定迭代格式，f X B X k k +=+)()1(*，（ ⋅⋅⋅=2,1,0k , k 即迭代次数），B 称为迭代矩阵。

最后选取初始迭代向量0X ，开始逐次迭代。

最后验证精度。

（迭代阵：b L D UXL D Xk k1)(1)1()()(--++++-=。

）Gauss-Seidel 迭代法与Jacobi 迭代法相比速度更快，但不全如此。

有例子表明：Gauss-Seidel 迭代法收敛时，Jacobi 迭代法可能不收敛；而Jacobi 迭代法收敛时，Gauss-Seidel 迭代法也可能不收敛。

三、Matlab 程序代码：function gauss_seidel(A,b,x0,eps,x1) D = diag(diag(A));%求A 的对角矩阵 L = -tril(A,-1);%求A 的下三角矩阵 U = -triu(A,1);%求A 的上三角矩阵 B = (D-L)\U; f = (D-L)\b; x = B*x0+f;n = 1;%迭代次数while norm(x1-x)>=eps x = B*x+f; n = n+1; endformat long n xjingdu=norm(x1-x)四、运行结果：实验四：一、实验名称及题目： Lagrange 插值法例4.1（P88）：给定函数f (x )=x(1+cosx)及插值节点x 0=0,x 1+π8,x 2=π4,x 3=3π8,x 4=π2.试构造Lagrange 插值多项式，给出其误差估计，并由此计算f (3π16)及其误差.二、解题思路：一般来说，如果我们有n 个点()()n n y x y x ,,...,1,1，各i x 互不相同。

那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：∑==nj j j x l y x L 0)()(，其中每个)(x l j 为拉格朗日基本多项式（或称插值基函数），其表达式为：∏≠=--=nj i i ij ij x x x x x l ,0)(。

三、Matlab程序代码：function y=lagrange(x0,x)n=length(x0);%向量长度s=0;for k=1:n%k从1到n的循环p=1.0;for j=1:nif j~=k%“~=”不等于的意思p=p*(x-x0(j))/(x0(k)-x0(j));endendy0=x0(k)*(1+cos(x0(k)));s=p*y0+s;endformat longswucha=abs(x*(1+cos(x))-s)四、运行结果：五、Lagrange插值图像绘制%Lagrange插值图像算法x=linspace(0,1002,200);s=linspace(0,1000,200);x0=[0;pi/8;pi/4;3*pi/8;pi/2];n=length(x0);s=0;for k=1:np=1.0;for j=1:nif j~=kp=p.*(x-x0(j))/(x0(k)-x0(j));endendy0=x0(k)*(1+cos(x0(k)));s=p*y0+s;endplot(x,s,'r');grid on;title('Lagrange²åֵͼÏñ')xlabel('X'),ylabel('Y');axis normal;实验五：一、实验名称及题目：Newton 插值法例 4.3（P96）：已知f (x )=√1+cℎ2x ,试取插值节点x 0=0.35,x 1=0.50,x 2=0.65,x 3=0.80,x 4=0.95,构造4次Newton 插值多项式，由此计算f (0.7)的逼近值，并指出其绝对误差. 二、解题思路：将拉格朗日插值公式∑∑∏==≠==--=nj j j j nj nj i i i j in y x l y x x x x x L 000)()(中的)(x L n 改写成：))...((...))(()()(10102010---++--=-+=n n n x x x x a x x x x a x x a a x N ,其中，n a a a ⋅⋅⋅,,10为待定定系数。

又)(00x f a =,[]0001)()(,x x x x f x f x x f a ---==，[][][]nn n n n x x x x x f x x x f x x x f a --==-,...,,,...,,,...,,10100。

将n a a a ⋅⋅⋅,,10带入)(x N n 可得：))...()(](,...,,[...)](,[)()(110100100----++-+=n n x x x x x x x x x f x x x x f x f x f 。

三、Matlab 程序代码：function newton_interpolation(x0,x) format long n=length(x0); syms zf=sqrt(1+cosh(z)^2); a(1)=subs(f,z,x0(1)); for k=1:n-1y0=subs(f,z,x0(k)); y1=subs(f,z,x0(k+1));d(k,1)=(y1-y0)/(x0(k+1)-x0(k));%一阶差商 endfor j=2:n-1for k=1:n-jd(k,j)=(d(k+1,j-1)-d(k,j-1))/(x0(k+j)-x0(k));%二阶差商及以上endenddouble(d)for j=2:na(j)=d(1,j-1);endb(1)=1;c(1)=a(1);for j=2:nb(j)=(x-x0(j-1)).*b(j-1);c(j)=a(j).*b(j);endnp=double(sum(c))wucha=double(abs(np-subs(f,z,x)))四、运行结果：五、Newton插值图像绘制实验六：一、实验名称及题目：Gauss 求积公式例5.7（P140）：试构造Gauss 型求积公式∫f(x)dx 1−1≈A 0f (x 0)+A 1f (x 1)+A 2f (x 2), 并由此计算积分∫√t (1+t)2dt 10. 二、解题思路：设高斯-勒让德求积公式是：()⎪⎪⎭⎫ ⎝⎛++⎪⎪⎭⎫ ⎝⎛-≈⎰53053)(21011-f A f A f A dx x f ,依次代入2,,1)(x x x f =,解得88,95120===A A A 。