计算方法课程设计.doc

数理学院 2014级信息与计算科学课 程 设 计姓名： 刘金玉 学号： 3141301240 班级： 1402 成绩： 实验要求1 应用自己熟悉的算法语言编写程序，使之尽可能具有通用性。2 上机前充分准备，复习有关算法，写出计算步骤，反复检查，调试程序。（注：在练习本上写，不上交）3 完成计算后写出实验报告，内容包括:算法步骤叙述，变量说明，程序清单，输出计算结果，结构分析和小结等。（注：具体题目具体分析，并不是所有的题目的实验报告都包含上述内容！）4 独立完成，如有雷同，一律判为零分！5 上机期间不允许做其他任何与课程设计无关的事情，否则被发现一次扣10分，被发现三次判为不及格！非特殊情况，不能请假。旷课3个半天及以上者，直接判为不及格。目 录一、基本技能训练41、误差分析42、求解非线性方程63、插值124、数值积分12二、提高技能训练161、162、18三、本课程设计的心得体会（500字左右）21一、基本技能训练1、误差分析实验1.3 求一元二次方程的根实验目的：研究误差传播的原因与解决对策。问题提出：求解一元二次方程实验内容：一元二次方程的求根公式为用求根公式求解下面两个方程：实验要求：（1） 考察单精度计算结果（与真解对比）；（2） 若计算结果与真解相差很大，分析其原因，提出新的算法（如先求再根据根与系数关系求）以改进计算结果。实验步骤：方程（1）：根据求根公式，写出程序：format longa=1;b=3;c=-2;x1=(-1)*b+sqrt(b2-4*a*c)/2*ax2=(-1)*b-sqrt(b2-4*a*c)/2*a运行结果：x1 = 0.561552812808830x2 = -3.561552812808830然后由符号解求的该方程的真值程序为：syms x;y=x2+3*x-2;s=solve(y,x);vpa(s)运行结果为：X1= 0.56155281280883027491070492798704X2= -3.561552812808830274910704927987方程（2）：format longa=1;b=-1010;c=1;x1=(-1)*b+sqrt(b2-4*a*c)/2*ax2=(-1)*b-sqrt(b2-4*a*c)/2*a运行结果为：x1 = 1.000000000000000e+010x2 = 0然后由符号解求的该方程的真值程序为：syms x;y=x2-1010*x+1;s=solve(y,x);vpa(s)运行结果：X1= 10000000000.0X2= 0.000000000116415321827由此可知，对于方程（1），使用求根公式求得的结果等于精确值，故求根公式法可靠。而对于方程（2），计算值与真值相差很大，故算法不可靠。改进算法，对于方程（2）:先用迭代法求x1，再用，求x2程序为：syms ka= ;for i=1:100 a(1)=4; a(i+1)=(1010*a(i)-1)(1/2);endx1=a(100) x2=1/(x1)运行结果为：x1 = 1.000000000000000e+010x2 = 1.000000000000000e-010实验结论：对于方程（1），两种方法在精确到小数点后15位时相同，说明两种算法的结果都是精确的。对于方程（2），两种算法结果有相当大的偏差，求根公式所求的一个根直接为零，求根公式的算法是不精确的。原因：方程（2）用求根公式计算时，公式中，b是大数，出现了大数吃掉小数的误差，也出现了两个相近的数相减的误差，所以出现x2=0这样大的误差。改进的结果会比较准确。2、求解非线性方程实验2.1 Gauss消去法的数值稳定性实验实验目的：观察和理解高斯消元过程中出现小主元即 很小时，引起方程组解的数值不稳定性实验内容：求解线性方程组 其中（1），（2）实验要求：（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的（2）用高斯列主元消去法求得L和U及解向量（3）用不选主元的高斯消去法求得L和U及解向量（4）观察小主元并分析对计算结果的影响实验步骤：（1） 计算矩阵的条件数程序：矩阵A1：A1=0.3*10(-15) 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1;cond(A1,1)cond(A1,2)cond(A1,inf)运行结果：ans = 136.294470872045e+000ans = 68.4295577125341e+000ans = 84.3114602051800e+000由矩阵条件数判断出矩阵A1是病态矩阵。矩阵A2：A2=10 -7 0 1;-3 2.09999999999 6 2;5 -1 5 -1;0 1 0 2;cond(A1,1)cond(A1,2)cond(A1,inf)运行结果：ans = 19.2831683168042e+000ans = 8.99393809015365e+000ans = 18.3564356435280e+000 由矩阵条件数判断出矩阵A2是病态矩阵。（2） 高斯列主元消去法程序：方程组（1）：A1=0.3*10(-15) 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1;b1=59.17;46.78;1;2;n=4;for k=1:n-1 a=max(abs(A1(k:n,k); p,k=find(A1=a); B=A1(k,:);c=b1(k); A1(k,:)=A1(p,:);b1(k)=b1(p); A1(p,:)=B;b1(p)=c; if A1(k,k)=0 A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k); A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n); else break endendL1=tril(A1,0);for i=1:n L1(i,i)=1;endL=L1U=triu(A1,0)for j=1:n-1 b1(j)=b1(j)/L(j,j); b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j);endb1(n)=b1(n)/L(n,n);for j=n:-1:2 b1(j)=b1(j)/U(j,j); b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j);endb1(1)=b1(1)/U(1,1);x1=b1运行结果： 方程组（2）：A2=10 -7 0 1;-3 2.09999999999 6 2;5 -1 5 -1;0 1 0 2;b2=8;5.9000000000001;5;1;n=4;for k=1:n-1 a=max(abs(A2(k:n,k); p,k=find(A2=a); B=A2(k,:);c=b2(k); A2(k,:)=A2(p,:);b2(k)=b2(p); A2(p,:)=B;b2(p)=c; if A2(k,k)=0 A2(k+1:n,k)=A2(k+1:n,k)/A2(k,k); A2(k+1:n,k+1:n)=A2(k+1:n,k+1:n)-A2(k+1:n,k)*A2(k,k+1:n); else break endendL1=tril(A2,0);for i=1:n L1(i,i)=1;endL=L1U=triu(A2,0)for j=1:n-1 b2(j)=b2(j)/L(j,j); b2(j+1:n)=b2(j+1:n)-b2(j)*L(j+1:n,j);endb2(n)=b2(n)/L(n,n);for j=n:-1:2 b2(j)=b2(j)/U(j,j); b2(1:j-1)=b2(1:j-1)-b2(j)*U(1:j-1,j);endb2(1)=b2(1)/U(1,1);x2=b2运行结果： （3） 不选主元的高斯消去法程序： 方程组（1）：clearformat longA1=0.3e-15 59.14 3 1;5.291 -6.130 -1 2;11.2 9 5 2;1 2 1 1;b1=59.17;46.78;1;2;n=4;for k=1:n-1 A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k); A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n);endL1=tril(A1,0);for i=1:n L1(i,i)=1;endL=L1U=triu(A1,0)for j=1:n-1 b1(j)=b1(j)/L(j,j); b1(j+1:n)=b1(j+1:n)-b1(j)*L(j+1:n,j);endb1(n)=b1(n)/L(n,n);for j=n:-1:2 b1(j)=b1(j)/U(j,j); b1(1:j-1)=b1(1:j-1)-b1(j)*U(1:j-1,j);endb1(1)=b1(1)/U(1,1);x1=b1运行结果: 方程组（2）：clearformat longA2=10 -7 0 1;-3 2.09999999999 6 2;5 -1 5 -1;0 1 0 2;b2=8;5.9000000000001;5;1;n=4;for k=1:n-1 A2(k+1:n,k)=A2(k+1:n,k)/A2(k,k); A2(k+1:n,k+1:n)=A2(k+1:n,k+1:n)-A2(k+1:n,k)*A2(k,k+1:n);endL1=tril(A2,0);for i=1:n L1(i,i)=1;endL=L1U=triu(A2,0)for j=1:n-1 b2(j)=b2(j)/L(j,j); b2(j+1:n)=b2(j+1:n)-b2(j)*L(j+1:n,j);endb2(n)=b2(n)/L(n,n);for j=n:-1:2 b2(j)=b2(j)/U(j,j); b2(1:j-1)=b2(1:j-1)-b2(j)*U(1:j-1,j);endb2(1)=b2(1)/U(1,1);x2=b2运行结果： （4） 分析小元对计算结果的影响通过观察计算结果，分析可知，小元对计算结果的影响很大，小元的存在会使得到的计算结果有很大的误差。3、插值4、数值积分实验2.4：复化求和公式计算定积分实验内容：计算下列各式右端定积分的近似值实验要求：（1） 分别用复化梯形公式、复化Simpson公式计算，要求绝对误差限为；（2） 将计算结果与精确解作比较，并比较各种算法的计算量。实验步骤：（1） 复化梯形公式和复化Simpson程序：式（1）：a=2;b=3;syms x;f=1/(x2-1);n=8;h=(b-a)/n;s=0;w=0;for i=1:n-1 x=a+h*i; s=s+subs(f,x)*2;endw=(s+subs(f,2)+subs(f,3)*h/2;f4=2*wf1=0;f2=0;for i=1:n-1 if rem(i,2)=0 x1=a+i*h; f1=f1+subs(f,x1); else rem(i,2)=0 x2=a+i*h; f2=f2+subs(f,x2); endendf3=(h/3)*(subs(f,2)+subs(f,3)+4*f1+2*f2)*2运行结果：f4 = 0.4064 f3 = 0.4055式(2):a=0;b=1;syms x;f=1/(x2+1);n=100;h=(b-a)/n;s=0;w=0;for i=1:n-1 x=a+h*i; s=s+subs(f,x)*2;endw=(s+subs(f,2)+subs(f,3)*h/2;f4=4*wf1=0;f2=0;for i=1:n-1 if rem(i,2)=0 x1=a+i*h; f1=f1+subs(f,x1); else rem(i,2)=0 x2=a+i*h; f2=f2+subs(f,x2); endendf3=(h/3)*(subs(f,2)+subs(f,3)+4*f1+2*f2)*4运行结果：f4 = 3.1176f3 = 3.1256式（3）：a=0;b=1;syms x;f=3x;n=100;h=(b-a)/n;s=0;w=0;for i=1:n-1 x=a+h*i; s=s+subs(f,x)*2;endw=(s+subs(f,2)+subs(f,3)*h/2;f4=wf1=0;f2=0;for i=1:n-1 if rem(i,2)=0 x1=a+i*h; f1=f1+subs(f,x1); else rem(i,2)=0 x2=a+i*h; f2=f2+subs(f,x2); endendf3=(h/3)*(subs(f,2)+subs(f,3)+4*f1+2*f2)运行结果：f4 = 1.9805f3 = 1.9271式（4）：a=1;b=2;syms x;f=x*exp(x);n=100;h=(b-a)/n;s=0;w=0;for i=1:n-1 x=a+h*i; s=s+subs(f,x)*2;endw=(s+subs(f,2)+subs(f,3)*h/2;f4=wf1=0;f2=0;for i=1:n-1 if rem(i,2)=0 x1=a+i*h; f1=f1+subs(f,x1); else rem(i,2)=0 x2=a+i*h; f2=f2+subs(f,x2); endendf3=(h/3)*(subs(f,2)+subs(f,3)+4*f1+2*f2)运行结果：f4 = 7.6769f3 = 7.5809（2） 求精确解程序：y1=log(2)y2=piy3=2/log(3)y4=exp(2)运行结果：y1 = 0.6931y2 = 3.1416y3 = 1.8205y4 = 7.3891（3） 分析使用复化梯形公式和复化Simpson公式计算所得的结果与真实值比较有很大的误差。两种公式算法中，复化梯形公式计算量比复化Simpson公式的计算量少，但是两种公式算法精确度都不是很高，有待使用其他精确度高的算法替代。二、提高技能训练1、kepler方程的计算在人造卫星轨道理论中对应的二体问题是可积系统，其中有 6 个积分常数为，其中 表示轨道半长径， 表示轨道偏心率， 表示轨道倾角， 是升交点赤经， 是近地点辐角，是平近点角。其中第六个积分常数，通常可以用其他的一些量替换，如过近地点时刻 ，真近点角度 f 与偏近点角度。这几个是相互等价的关系，我们这里要讲的就是平近点角 M 与偏近点角 的一个关系，有如下方程这就是著名的 Kelper 方程，是在人造卫星运动理论或者行星运动理论中最基本的方程之一。 因为如果我们已经知道卫星轨道量去计算卫星的星历时， 第 6 个根数通常是使用平近点角。这时候需要把平近点角化为偏近点角，也就是需要计算上述的 Kepler 方程。具体理论这里不再阐述，而 Kepler 方程本身是我们这里感兴趣的。可采用 Newton 迭代法计算 Kepler 方程，在实际工程中，一般也是这样做的，因为 Newton 迭代法收敛速度快而且精度比较高。也可用其它算法计算。实验目的与要求（1）了解 Kepler 方程；（2）能够用牛顿迭代方法计算 Kepler 方程；（3）通过调节轨道偏心率 e 查看运算收敛情况。实验内容及数据来源编写解 kepler 方程的方法，给定偏心率为 0.01、平近点角为 32 度时计算出偏近点角。实验步骤：（1） 对Kepler 方程的了解： Kepler方程又称作开普勒方程，是二体问题运动方程的一个积分。它反映天体在其轨道上的位置与时间t的函数关系。对于椭圆轨道，开普勒方程可以表示为EesinE=M，式中E为偏近点角，M为平近点角，都是从椭圆轨道的近地点开始起算，沿逆时针方向为正，E和M都是确定天体在椭圆轨道上的运动和位置的基本量。(2) 牛顿迭代法计算开普勒方程程序：N =100000;x0=0;y=32-x0+0.01*sin(x0);E=1.0e-6;k=1;while(norm(y-x0)=E&kN) x0=y; f=diff(y); f1=diff(f); if(f1=0) y=x0-f/f1; end k=k+1;endfprintf(迭代的次数为k=%dn,k);fprintf(n);fprintf(方程的根为 x*=%.7fn,y);运行结果：迭代的次数为k=2方程的根为 x*=32.0000000（3） 调节偏心率查看收敛情况： 调节不用的偏心率的运行结果都趋于32，由此判断出该方程为收敛方程。实验结论：调节不同的偏心率程序运行出来的结果没有很大的偏差，所以该方程是收敛的。2（1）、一维插值问题的应用（1）（机翼加工问题）书籍机翼轮廓上的数据如下表所示 x/m035791112131415y/m01.21.72.02.12.01.81.21.01.6加工时需要x每改变0.1m时y值，并画出相应的轮廓曲线。试验程序：x0=0 3 5 7 9 11 12 13 14 15;y0=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6;x=0:0.1:15;y1=interp1(x0,y0,x); %分段线性插值y2=interp1(x0,y0,x,spline); %三次样条插值subplot(2,1,1) plot(x0,y0,k+,x,y1,r)gridtitle(piecewise linear) %图片命名subplot(2,1,2)plot(x0,y0,k+,x,y2,r)gridtitle(spline)运行结果：（2）、二维插值问题的应用 (山区地貌图)在某山区（平面区域内，单位:m）测得一些地点的高度（单位:m）如下表所示。240014301450147013201280120010809402000145014801500155015101430130012001600146015001550160015501600160016001200137015001200110015501600155013808001270150012001100135014501200115040012301390150015001400900110010600118013201450142014001300700900y/x040080012001600200024002800试作出该山区的地貌图和等值线图。试验程序：x=0:400:2800;y=0:400:2400;h=1180 1320 1450 1420 1400 1300 700 900 ; 1230 1390 1500 1500 1400 900 1100 1060; 1270 1500 1200 1100 1350 1450 1200 1150; 1370 1500 1200 1100 1550 1600 1550 1380; 1460 1500 1550 1600 1550 1