计算方法A上机实验报告.doc

计算方法A上机实验报告姓名：苏福 班级：硕4020 学号：3114161019一、上机练习目的1）复习和巩固数值计算方法的基本数学模型，全面掌握运用计算机进行数值计算的具体过程及相关问题。2）利用计算机语言独立编写、调试数值计算方法程序，培养学生利用计算机和所学理论知识分析解决实际问题的能力。二、上机练习任务1）利用计算机语言编写并调试一系列数值方法计算通用程序，并能正确计算给定题目，掌握调试技能。2） 掌握文件使用编程技能，如文件的各类操作，数据格式设计、通用程序运行过程中文件输入输出运行方式设计等。3）写出上机练习报告。三、上机题目1. 共轭梯度法求解线性方程组。 （第三章）2. 三次样条插值（第四章）3. 龙贝格积分（第六章）4. 四阶龙格库塔法求解常微分方程的初值问题四、上机报告题目1：共轭梯度法求解线性方程组1 算法原理共轭梯度法是把求解线性方程组的问题转化为求解一个与之等价的二次函数极小值的问题。从任意给定的初始点出发，沿一组关于矩阵共轭的方向进行线性搜索，在无舍入误差的假定下，最多迭代次（其中为矩阵的阶数），就可求得二次函数的极小值，也就求得了线性方程组的解。 定理：设是阶对称正定矩阵，则是方程组的解得充分必要条件是是二次函数的极小点，即共轭梯度法的计算公式：2. 程序框图3. MATLAB编程实现（1）编写共轭梯度法求解对称正定矩阵的线性方程组见附录（myge.m）：function x=myge(A,b)输入对称正定矩阵及对应的列向量，初始向量设为0，精度取为。函数的输出即为由共轭梯度发求解的近似解。（2）编写具体算例求解（example.m）：clcclear all%例题3.4.2A0=2,0,1;0,1,0;1,0,2;b0=3,1,3;myge(A0,b0);%习题3.2n=100; %矩阵阶数A=zeros(n,n);b=zeros(n,1);b(1)=-1;b(n)=-1;A(1,1)=-2;A(1,2)=1;A(n,n-1)=1;A(n,n)=-2;for i=2:n-1 A(i,i-1)=1; A(i,i)=-2; A(i,i+1)=1;endmyge(A,b);算例1（课本例题3.4.2）：算例2（课后习题计算实习3.2）：4. 算例结果算例1：x = 1.000000000000000 1.000000000000000 1.000000000000000迭代次数：k = 2算例2(n=100)：x = 0.999999999999999 1.000000000000005 0.999999999999982 1.000000000000038 0.999999999999976 1.000000000000013 1.000000000000000 1.000000000000008 0.999999999999984 1.000000000000017 0.999999999999977 1.000000000000052 0.999999999999942 1.000000000000058 0.999999999999974 1.000000000000011 0.999999999999992 1.000000000000028 0.999999999999983 1.000000000000008 1.000000000000021 0.999999999999977 1.000000000000027 0.999999999999987 1.000000000000031 0.999999999999966 1.000000000000043 0.999999999999985 1.000000000000019 0.999999999999993 1.000000000000022 0.999999999999990 1.000000000000019 1.000000000000006 1.000000000000004 1.000000000000011 0.999999999999996 1.000000000000022 0.999999999999991 1.000000000000026 0.999999999999996 1.000000000000018 1.000000000000004 1.000000000000014 1.000000000000012 1.000000000000008 1.000000000000009 1.000000000000015 1.000000000000010 1.000000000000012 1.000000000000012 1.000000000000010 1.000000000000015 1.000000000000007 1.000000000000015 1.000000000000007 1.000000000000012 1.000000000000007 1.000000000000019 1.000000000000001 1.000000000000023 0.999999999999992 1.000000000000021 0.999999999999997 1.000000000000021 0.999999999999984 1.000000000000029 0.999999999999985 1.000000000000030 0.999999999999986 1.000000000000026 0.999999999999974 1.000000000000054 0.999999999999951 1.000000000000062 0.999999999999943 1.000000000000078 0.999999999999934 1.000000000000060 0.999999999999968 1.000000000000021 0.999999999999990 1.000000000000019 0.999999999999999 1.000000000000012 0.999999999999978 1.000000000000042 0.999999999999963 1.000000000000057 0.999999999999933 1.000000000000063 0.999999999999960 1.000000000000034 0.999999999999974 1.000000000000025 0.999999999999969 1.000000000000054 0.999999999999955 1.000000000000022 0.999999999999988迭代次数：k =50对于n=200,400，编写的程序的计算结果显示可靠，满足精度要求。题目2: 三次样条插值1. 算法原理设在区间上给定个节点，在节点处的函数值为。若函数满足如下三条：（1） 在每个子区间上，是三次多项式；（2） ；（3） 在区间上，的二阶导数连续。则称为函数在区间上的三次样条插值函数。子区间上的的表达式为：关于参数的方程组（三弯矩方程组）：牛顿插值多项式：构造牛顿插值多项式首先列出差商表，进而由差商表写出牛顿插值多项式。阶差商为：零阶差商和一阶差商：2. 程序框图3. MATLAB编程实现给定函数，取等距节点。编写三次样条插值函数见附录（myspline.m）：function Sx=myspline(n)输入为区间个数，输出为各区间上的三次样条插值函数。编写牛顿插值函数见附录(myNewton.m)：function Nx=myNewton(n)输入为区间个数，输出为牛顿插值函数。编写main函数(main.m)：clcclear allx=-1:0.001:1;y=1./(1+25*x.2);plot(x,y,k-);hold on;Nx5=subs(myNewton(5),x);plot(x,Nx5,c-.);hold on;Nx10=subs(myNewton(10),x);plot(x,Nx10,r-);hold on;myspline(10);hold onlegend(f(x)=1/(1+25x2),f(x)=N5(x),f(x)=N10(x),f(x)=S10(x)4. 算例结果n = 5次牛顿插值多项式 Nx =0.1538*x + 1.058*(x + 1.0)*(x + 0.6) - 1.923*(x + 1.0)*(x + 0.6)*(x + 0.2) + 1.202*(x + 1.0)*(x + 0.6)*(x + 0.2)*(x - 0.2) + 0.1923 n = 10次牛顿插值多项式 Nx =0.1018*x + 0.2602*(x + 1.0)*(x + 0.8) + 0.7919*(x + 1.0)*(x + 0.6)*(x + 0.8) + 2.687*(x + 1.0)*(x + 0.6)*(x + 0.4)*(x + 0.8) - 6.363*(x + 1.0)*(x + 0.6)*(x + 0.2)*(x + 0.4)*(x + 0.8) - 17.68*x*(x + 1.0)*(x + 0.6)*(x + 0.2)*(x + 0.4)*(x + 0.8) + 84.84*x*(x + 1.0)*(x + 0.6)*(x + 0.2)*(x + 0.4)*(x + 0.8)*(x - 0.2) - 167.9*x*(x + 1.0)*(x + 0.6)*(x + 0.2)*(x + 0.4)*(x + 0.8)*(x - 0.2)*(x - 0.4) + 220.9*x*(x + 1.0)*(x + 0.6)*(x + 0.2)*(x + 0.4)*(x + 0.8)*(x - 0.2)*(x - 0.4)*(x - 0.6) - 220.9*x*(x + 1.0)*(x + 0.6)*(x + 0.2)*(x + 0.4)*(x + 0.8)*(x - 0.2)*(x - 0.4)*(x - 0.8)*(x - 0.6) + 0.1403n = 10次三次样条插值多项式 Sx = 0.08814*x + 0.3417*(x + 1.0)3 + 0.1266 0.1702*x - 0.3417*(x + 0.6)3 + 1.235*(x + 0.8)3 + 0.1922 0.4665*x + 2.071*(x + 0.6)3 - 1.235*(x + 0.4)3 + 0.37 0.9637*x - 2.071*(x + 0.2)3 + 15.48*(x + 0.4)3 + 0.5689 4.679*x - 38.99*(x + 0.2)3 - 15.48*x3 + 1.312 38.99*(x - 0.2)3 - 4.679*x + 15.48*x3 + 1.312 2.071*(x - 0.2)3 - 0.9637*x - 15.48*(x - 0.4)3 + 0.5689 1.235*(x - 0.4)3 - 0.4665*x - 2.071*(x - 0.6)3 + 0.37 0.3417*(x - 0.6)3 - 1.235*(x - 0.8)3 - 0.1702*x + 0.1922 0.1266 - 0.3417*(x - 1.0)3 - 0.08814*x 三种插值多项式与的曲线在同一坐标系:题目3：龙贝格积分1 算法原理通过逐步缩小步长的方法，以使积分的近似值满足精度要求。在实际中，通过将区间间逐次对半分来实现。在变步长积分的思想上，对于复化梯形积分：对于复化辛普森求积公式：同理可以证明：再令：即为龙贝格积分公式。龙贝格积分是将区间逐次分半，逐次递推计算，以得到较高精度的积分近似值。龙贝格的积分公式如下：若2. 程序框图3. MATLAB编程实现（1）编写自定义函数（f.m）：function fun=f(x)fun=log(1+x)/x; %改变函数表达式求不同的函数的积分取习题6.2的积分函数为例采用龙贝格积分方法求积分。（2）编写龙贝格积分程序见附录（myromberg.m）：4. 算例结果（1）R_k = 0.693147978563680（2）R_k = 0.272198867426936（3）R_k = 0.822467764886147（4）R_k = 1.370762774269822题目4：四阶龙格库塔法求解常微分方程的初值问题1. 算法原理对于高阶常微分方程的初值问题：将其转化为一阶微分方程组求解：标准的四级四阶龙格库塔法的向量形式：其分量形式：2. 程序框图3. MATLAB编程实现编写4级4阶龙格库塔积分法求解（步长0.05）：转化为一阶常微分方程组：编写MATLAB程序myRK.m（见附录），求出总长度为5的数值积分结果，并与解析解作比较。4. 算例结果龙格-库塔积分结果与实际函数值比较。 五、上机总结通过本次计算方法A课程的上机实验，复习和巩固数值计算方法的基本数学模型， 对共轭梯度法求解线性方程组，三次样条插值，牛顿插值，龙贝格积分法，标准四级四阶龙格库塔法求解高阶常微分方程的初值问题等知识有深刻的理解与认识，并通过MATLAB编程实现各种数值算法，对数值分析的应用有清晰的认识。附录1. 共轭梯度法myge.m%共轭梯度法function x=myge(A,b)length=size(b,1);x0=zeros(length,1); %迭代初始向量r0=b-A*x0;d0=r0;for k=1:length a=r0*r0/(d0*A*d0); x0=x0+a*d0; r1=b-A*x0; if norm(r1)1e-6 %精度设为1e-6 R_1=R_k; for i=1:2k Tsum=Tsum+f(a+(2*i-1)*(b-a)/2(k+1); end T_k=T_8/2+(b-a)/2(k+1)*Tsum; S_k=T_k+1/(4-1)*(T_k-T_8); C_k=S_k+1/(42-1)*(S_k-S_4); R_k=C_k+1/(43-1)*(C_k-C_2); T_8=T_k; S_4=S_k; C_2=C_k; Tsum=0; k=k+1;endkformat longR_k5. 四阶龙格库塔法myRK.m%myRK.mclcclear allh=0.05;n=5/h;y=zeros(3,n);y(:,1)=-1;3;2;K1=zeros(3,1);K2=zeros(3,1);K3=zeros(3,1);K4=zeros(3,1);for i=1:n x=i*h; for j=1:3 K1(1)=h*y(2,i);K1(2)=h*y(3,i);K1(3)=h*(y(3,i)+y(2,i)-y(1,i)+2*x-3); K2(1)=h*(y(2,i)+K1(2)/2); K2(2)=h*(y(3,i)+K1(3)/2); K2(3)=h*(y(3,i)+K1(3)/2+y(2,i)+K1(2)/2-(y(1,i)+K1(1)/2)+2*(x+h/2)-3); K3(1)=h*(y(2,i)+K2(2)/2); K3(2)=h*(y(3,i)+K2(3)/2); K3(3)=h*(y(3,i)+K2(3)/2+y(2,i)+K2(2)/2-(y(1,i)+K2(1)/2)+2*(x+h/2)-3); K4(1)=h*(y(2,i)+K3(2); K4(2)=h*(y(3,i)+K3(3); K4(3)=h*(y(3,i)+K3(3)+y(2,i)+K3(2)-(y(1,i)+K3(1)+2*(x+h)-3); end for k=1:3 y(k,i+1)=y(k,i)+(K1(k)+2*K2(k)+2*K3(k)+K4(k)/6; endendxh=h:h:(n+1)*h;plot(xh,y(1,:),r.);hold onxi=h:0.01:(n+1)*h;yi=xi.*exp(xi)+2*xi-1;plot(xi,yi,k);hold onlegend(4级R-K,y(x)=xex+2x-1)