计算方法上机实验报告

《计算方法》上机实验报告班级：XXXXXX小组成员：XXXXXXXXXXXXXXXXXXXXXXXXXXXX任课教师：XXX二〇一八年五月二十五日-2-前言通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton迭代法、Jacobi迭代法、Gauss-Seidel迭代法、Newton插值法、Lagrange插值法和Gauss求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB程序的编写。以下为本次上机实验报告，按照实验内容共分为六部分。实验一：一、实验名称及题目：Newton迭代法例2.7(P38):应用Newton迭代法求𝑥3−𝑥−1=0在𝑥=1附近的数值解𝑥𝑘，并使其满足|𝑥𝑘−𝑥𝑘−1|10−8.二、解题思路：设'x是0)(xf的根,选取0x作为'x初始近似值,过点)(,00xfx做曲线)(xfy的切线L,L的方程为))((')(000xxxfxfy,求出L与x轴交点的横坐标)(')(0001xfxfxx,称1x为'x的一次近似值,过点))(,(11xfx做曲线)(xfy的切线,求该切线与x轴的横坐标)(')(1112xfxfxx称2x为'x-3-的二次近似值,重复以上过程,得'x的近似值序列nx,把)(')(1nnnnxfxfxx称为'x的1n次近似值，这种求解方法就是牛顿迭代法。三、Matlab程序代码：functionnewton_iteration(x0,tol)symsz%定义自变量formatlong%定义精度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;whileabs(x1-x0)=tolx0=x1;y=subs(f,z,x0);y1=subs(f1,z,x0);x1=x0-y/y1;k=k+1;endx=double(x1)K四、运行结果：实验二：-4-一、实验名称及题目：Jacobi迭代法例3.7（P74）：试利用Jacobi迭代公式求解方程组[5−1−1−1−110−1−1−1−15−1−1−1−110][𝑥1𝑥2𝑥3𝑥4]=[−412834]要求数值解𝑋(𝑘)满足||X−𝑋(𝑘)||2≤10−4，其中𝑋=(1,2,3,4)𝑇为方程组的精确解.二、解题思路：首先将方程组中的系数矩阵A分解成三部分，即：UDLA，D为对角阵，L为下三角矩阵，U为上三角矩阵。之后确定迭代格式，fXBXkk)()1(*，（2,1,0k,k即迭代次数），B称为迭代矩阵。最后选取初始迭代向量)0(X，开始逐次迭代。最后验证精度。（迭代阵：bDUXDXkk1)(1)1(。）雅克比迭代法的优点明显，计算公式简单，每迭代一次只需计算一次矩阵和向量的乘法，且计算过程中原始矩阵A始终不变，比较容易并行计算。然而这种迭代方式收敛速度较慢，而且占据的存储空间较大。三、Matlab程序代码：functionjacobi(A,b,x0,eps,x1)D=diag(diag(A));%求A的对角矩阵L=-tril(A,-1);%求A的下三角矩阵-5-U=-triu(A,1);%求A的上三角矩阵B=D\(L+U);f=D\b;x=B*x0+f;n=1;%迭代次数whilenorm(x-x1)=epsx=B*x+f;n=n+1;endformatlongnxjingdu=norm(x-x1)四、运行结果：实验三：一、实验名称及题目：Gauss-Seidel迭代法例3.8（P75）:试利用Gauss-Seidel迭代公式求解方程组[5−1−1−1−110−1−1−1−15−1−1−1−110][𝑥1𝑥2𝑥3𝑥4]=[−412834]，并使其数值解𝑋(𝑘)满足精度要求||X−𝑋(𝑘)||2≤10−4，其中𝑋=(1,2,3,4)𝑇为方程组的精确解.-6-二、解题思路：Gauss-Seidel迭代法与Jacobi迭代法思路相近，首先将方程组中的系数矩阵A分解成三部分，即：UDLA，D为对角阵，L为下三角矩阵，U为上三角矩阵。之后确定迭代格式，fXBXkk)()1(*，（2,1,0k,k即迭代次数），B称为迭代矩阵。最后选取初始迭代向量0X，开始逐次迭代。最后验证精度。（迭代阵：bLDUXLDXkk1)(1)1()()(。）Gauss-Seidel迭代法与Jacobi迭代法相比速度更快，但不全如此。有例子表明：Gauss-Seidel迭代法收敛时，Jacobi迭代法可能不收敛；而Jacobi迭代法收敛时，Gauss-Seidel迭代法也可能不收敛。三、Matlab程序代码：functiongauss_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;%迭代次数whilenorm(x1-x)=epsx=B*x+f;n=n+1;endformatlongnxjingdu=norm(x1-x)-7-四、运行结果：实验四：一、实验名称及题目：Lagrange插值法例4.1（P88）：给定函数f(x)=x(1+cosx)及插值节点𝑥0=0,𝑥1+𝜋8,𝑥2=𝜋4,𝑥3=3𝜋8,𝑥4=𝜋2.试构造Lagrange插值多项式，给出其误差估计，并由此计算f(3𝜋16)及其误差.二、解题思路：一般来说，如果我们有n个点nnyxyx,,...,1,1，各ix互不相同。那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：njjjxlyxL0)()(，其中每个)(xlj为拉格朗日基本多项式（或称插值基函数），其表达式为：njiiijijxxxxxl,0)(。-8-三、Matlab程序代码：functiony=lagrange(x0,x)n=length(x0);%向量长度s=0;fork=1:n%k从1到n的循环p=1.0;forj=1:nifj~=k%“~=”不等于的意思p=p*(x-x0(j))/(x0(k)-x0(j));endendy0=x0(k)*(1+cos(x0(k)));s=p*y0+s;endformatlongswucha=abs(x*(1+cos(x))-s)四、运行结果：-9-五、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;fork=1:np=1.0;forj=1:nifj~=kp=p.*(x-x0(j))/(x0(k)-x0(j));endendy0=x0(k)*(1+cos(x0(k)));s=p*y0+s;endplot(x,s,'r');gridon;title('Lagrange²åֵͼÏñ')xlabel('X'),ylabel('Y');axisnormal;-10-实验五：一、实验名称及题目：Newton插值法例4.3（P96）：已知f(x)=√1+𝑐ℎ2𝑥,试取插值节点𝑥0=0.35,𝑥1=0.50,𝑥2=0.65,𝑥3=0.80,𝑥4=0.95,构造4次Newton插值多项式，由此计算f(0.7)的逼近值，并指出其绝对误差.二、解题思路：将拉格朗日插值公式njjjjnjnjiiijinyxlyxxxxxL000)()(中的)(xLn改写成：))...((...))(()()(10102010nnnxxxxaxxxxaxxaaxN,其中，naaa,,10为待定定系数。又)(00xfa,0001)()(,xxxxfxfxxfa，nnnnnxxxxxfxxxfxxxfa,...,,,...,,,...,,10100。将naaa,,10带入)(xNn可得：))...()(](,...,,[...)](,[)()(110100100nnxxxxxxxxxfxxxxfxfxf。三、Matlab程序代码：functionnewton_interpolation(x0,x)formatlongn=length(x0);symszf=sqrt(1+cosh(z)^2);a(1)=subs(f,z,x0(1));fork=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));%一阶差商end-11-forj=2:n-1fork=1:n-jd(k,j)=(d(k+1,j-1)-d(k,j-1))/(x0(k+j)-x0(k));%二阶差商及以上endenddouble(d)forj=2:na(j)=d(1,j-1);endb(1)=1;c(1)=a(1);forj=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插值图像绘制-12-实验六：一、实验名称及题目：Gauss求积公式例5.7（P140）：试构造Gauss型求积公式∫𝑓(𝑥)dx1−1≈𝐴0𝑓(𝑥0)+𝐴1𝑓(𝑥1)+𝐴2𝑓(𝑥2),并由此计算积分∫√𝑡(1+𝑡)2dt10.二、解题思路：设高斯-勒让德求积公式是：53053)(21011-fAfAfAdxxf,依次代入2,,1)(xxxf,解得88,95120AAA。利用换元公式badxabxabfabdttf11222)(变换原式的积分上下限，在套用高斯-勒让德求积公式求得积分。三、Matlab程序代码：functiongauss(a,b)symstf=sqrt(t)/(1+t)^2;P=[-sqrt(3/5)0sqrt(3/5)];A=[5/98/95/9];s=0;fori=1:3;x=P(i);y=subs(f,t,(b-a)*x/2+(a+b)/2);s=s+A(i)*y;endformatlongS=double((a-b)/2*s)-13-四、运行结果：结束语在本学期的《计算方法》课程学习中，我们感受到了巧妙的数学计算方法在解决实际问题中带来的便利与高效，借助计算机解决科学计算问题也是我们当代大学生应当掌握的必要技能。在本课程的上机实验过程中，我们亲自体验了课本中所学到的算法在计算机上的使用，使用计算机语言，如MATLAB和C语言等进行编程，加深了我们对各类算法及其理论的理解，并进一步激发了我们对《计算方法》这门课程进行持续学习的学习兴趣。