1、1数理学院 2014 级信息与计算科学课 程 设 计姓名： 刘金玉 学号： 3141301240 班级： 1402 成绩： 2实验要求1 应用自己熟悉的算法语言编写程序，使之尽可能具有通用性。2 上机前充分准备，复习有关算法，写出计算步骤，反复检查，调试程序。 （注：在练习本上写，不上交）3 完成计算后写出实验报告，内容包括:算法步骤叙述，变量说明，程序清单，输出计算结果，结构分析和小结等。 （注：具体题目具体分析，并不是所有的题目的实验报告都包含上述内容！）4 独立完成，如有雷同，一律判为零分！5 上机期间不允许做其他任何与课程设计无关的事情，否则被发现一次扣 10 分，被发现三次判为不及格
2、！非特殊情况，不能请假。旷课 3 个半天及以上者，直接判为不及格。3目 录一、基本技能训练 41、误差分析 .42、求解非线性方程 .63、插值 .124、数值积分 .12二、提高技能训练 161、 .162、 .18三、本课程设计的心得体会（500 字左右） .214一、基本技能训练1、误差分析实验 1.3 求一元二次方程的根实验目的：研究误差传播的原因与解决对策。问题提出：求解一元二次方程 20axbc实验内容：一元二次方程的求根公式为 21,24bcx用求根公式求解下面两个方程： 210()3实验要求：（1） 考察单精度计算结果（与真解对比） ；（2） 若计算结果与真解相差很大，分析其原
3、因，提出新的算法（如先求 再1x根据根与系数关系求 ）以改进计算结果。2x实验步骤：方程（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*a5运行结果：x1 = 0.561552812808830x2 = -3.561552812808830然后由符号解求的该方程的真值程序为：syms x;y=x2+3*x-2;s=solve(y,x);vpa(s)运行结果为：X1= 0.56155281280883027491070492798704X2= -3.561
4、552812808830274910704927987方程（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） ，使用求根公式求得的结果等于精确值，故求根公式法可
5、靠。而对于方程（2） ，计算值与真值相差很大，故算法不可靠。改进算法，对于方程（2）:先用迭代法求 x1，再用 ，求 x21axc程序为：syms k6a= ;for i=1:100a(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） ，两种算法结果有相当大的偏差，求根公式所求的一个根直接为零，求根公式的算
6、法是不精确的。原因：方程（2）用求根公式计算时， 公式中，b 是大数，出现了24bac大数吃掉小数的误差，也出现了两个相近的数相减的误差，所以出现 x2=0 这样大的误差。改进的结果会比较准确。2、求解非线性方程实验 2.1 Gauss 消去法的数值稳定性实验实验目的：观察和理解高斯消元过程中出现小主元即 很小时，引起方程组()ka解的数值不稳定性实验内容：求解线性方程组 其中(1,2)iiAxb（1） ，1510.39.43260.1A159.746827（2） 2107013.962515A 85.901b实验要求：（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的（2）用高斯列主元消
7、去法求得 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 -
8、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;8b1=59.17;46.78;1;2;n=4;for k=1:
9、n-1a=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)=0A1(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);elsebreakendendL1=tril(A1,0);for i=1:nL1(i,i)=1;endL=L1U=triu(A1,0)for j=1:n-1b1(j)=b1(j)/L(j,j);b
10、1(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:2b1(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运行结果：9180026.7904.751.32.49L 1.29520431.8.0U18;.;3865x方程组（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
11、=1:n-1a=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)=0A2(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);elsebreakendendL1=tril(A2,0);for i=1:nL1(i,i)=1;endL=L1U=triu(A2,0)for j=1:n-1b2(j)=b2(j)/L(j,j
12、);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:210b2(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运行结果：12100.53.4.3L10712.5.630.U120.;0;x（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
13、.17;46.78;1;2;n=4;for k=1:n-1A1(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:nL1(i,i)=1;endL=L1U=triu(A1,0)for j=1:n-1b1(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:211b1(j)=b1(j)/U(j,j);b1(1:j
14、-1)=b1(1:j-1)-b1(j)*U(1:j-1,j);endb1(1)=b1(1)/U(1,1);x1=b1运行结果:151507.63 2.681.9L 151815150.3.4302.97.60 8.U 123.684;1.05;x方程组（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-1A2(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
15、)-A2(k+1:n,k)*A2(k,k+1:n);endL1=tril(A2,0);for i=1:nL1(i,i)=1;endL=L1U=triu(A2,0)12for j=1:n-1b2(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:2b2(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运行结果：12100.35.498.4L 1212121070.6
16、.3 .9875904.67U 520.361;.0;.x（4） 分析小元对计算结果的影响通过观察计算结果，分析可知，小元对计算结果的影响很大，小元的存在会使得到的计算结果有很大的误差。3、插值4、数值积分实验 2.4：复化求和公式计算定积分实验内容：计算下列各式右端定积分的近似值 312 201(1) ln () =4dxdx131 20 12(3) (4) elnx xdd实验要求：（1） 分别用复化梯形公式、复化 Simpson 公式计算，要求绝对误差限为；7.51（2） 将计算结果与精确解作比较，并比较各种算法的计算量。实验步骤：（1） 复化梯形公式和复化 Simpson 程序：式（1
17、）：a=2;b=3;syms x;f=1/(x2-1);n=8;h=(b-a)/n;s=0;w=0;for i=1:n-1x=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-1if rem(i,2)=0x1=a+i*h;f1=f1+subs(f,x1);else rem(i,2)=0x2=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 14f3 = 0.40
18、55式(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-1x=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-1if rem(i,2)=0x1=a+i*h;f1=f1+subs(f,x1);else rem(i,2)=0x2=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 =
19、3.1256式（3）：a=0;b=1;syms x;f=3x;n=100;h=(b-a)/n;15s=0;w=0;for i=1:n-1x=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-1if rem(i,2)=0x1=a+i*h;f1=f1+subs(f,x1);else rem(i,2)=0x2=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.92
20、71式（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-1x=a+h*i;s=s+subs(f,x)*2;endw=(s+subs(f,2)+subs(f,3)*h/2;f4=w16f1=0;f2=0;for i=1:n-1if rem(i,2)=0x1=a+i*h;f1=f1+subs(f,x1);else rem(i,2)=0x2=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.
21、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 方程的计算17在人造卫星轨道理论中对应的二体问题是可积系统，其中有 6 个积分常数为 ，其中 表示轨道半长径， 表示轨道偏心率， 表示
22、轨(,)TaeiMaei道倾角， 是升交点赤经， 是近地点辐角， 是平近点角。其中第六个M积分常数，通常可以用其他的一些量替换，如过近地点时刻 ，真近点角度 f pt与偏近点角度 。这几个是相互等价的关系，我们这里要讲的就是平近点角 M E与偏近点角 的一个关系，有如下方程 sinEe这就是著名的 Kelper 方程，是在人造卫星运动理论或者行星运动理论中最基本的方程之一。 因为如果我们已经知道卫星轨道量去计算卫星的星历时， 第 6 个根数通常是使用平近点角。这时候需要把平近点角化为偏近点角，也就是需要计算上述的 Kepler 方程。具体理论这里不再阐述，而 Kepler 方程本身是我们这里感
23、兴趣的。可采用 Newton 迭代法计算 Kepler 方程，在实际工程中，一般也是这样做的，因为 Newton 迭代法收敛速度快而且精度比较高。也可用其它算法计算。实验目的与要求（1）了解 Kepler 方程；（2）能够用牛顿迭代方法计算 Kepler 方程；（3）通过调节轨道偏心率 e 查看运算收敛情况。实验内容及数据来源编写解 kepler 方程的方法，给定偏心率为 0.01、平近点角为 32 度时计算出偏近点角 。E实验步骤：（1）对 Kepler 方程的了解：Kepler 方程又称作开普勒方程，是二体问题运动方程的一个积分。它反映天体在其轨道上的位置与时间 t 的函数关系。对于椭圆轨
24、道，开普勒方程可以表示为 EesinE=M，式中 E 为偏近点角，M 为平近点角，都是从椭圆轨道的近地点开始起算，沿逆时针方向为正，E 和 M 都是确定天体在椭圆轨道上的运动和位置的基本量。(2) 牛顿迭代法计算开普勒方程程序：N =100000;x0=0;y=32-x0+0.01*sin(x0);18E=1.0e-6;k=1;while(norm(y-x0)=Ef=diff(y);f1=diff(f);if(f1=0)y=x0-f/f1;endk=k+1;endfprintf(迭代的次数为 k=%dn,k);fprintf(n);fprintf(方程的根为 x*=%.7fn,y);运行结果：
25、迭代的次数为 k=2方程的根为 x*=32.0000000（3）调节偏心率查看收敛情况：调节不用的偏心率的运行结果都趋于 32，由此判断出该方程为收敛方程。实验结论：调节不同的偏心率程序运行出来的结果没有很大的偏差，所以该方程是收敛的。2（ 1） 、一维插值问题的应用（1） （机翼加工问题）书籍机翼轮廓上的数据如下表所示 x/m 0 3 5 7 9 11 12 13 14 15y/m 0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6加工时需要 x 每改变 0.1m 时 y 值，并画出相应的轮廓曲线。试验程序：x0=0 3 5 7 9 11 12 13 14 15;y0
26、=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); %三次样条插值19subplot(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)运行结果：0 5 10 150123 piecewise linear0 5 10 150123 spline（2） 、二维插值
27、问题的应用(山区地貌图)在某山区（平面区域 内，单位:m）测得一些,28,4地点的高度（单位:m）如下表所示。2400 1430 1450 1470 1320 1280 1200 1080 9402000 1450 1480 1500 1550 1510 1430 1300 12001600 1460 1500 1550 1600 1550 1600 1600 16001200 1370 1500 1200 1100 1550 1600 1550 1380800 1270 1500 1200 1100 1350 1450 1200 1150400 1230 1390 1500 1500 140
28、0 900 1100 10600 1180 1320 1450 1420 1400 1300 700 900y/x 0 400 800 1200 1600 2000 2400 2800试作出该山区的地貌图和等值线图。试验程序：20x=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
29、;1460 1500 1550 1600 1550 1600 1600 1600;1450 1480 1500 1550 1510 1430 1300 1200;1430 1450 1470 1320 1280 1200 1080 940;xi=0:50:2800;xi=xi;yi=0:50:2400;z=interp2(x,y,h,xi,yi,cubic);figure(1)mesh(xi,yi,z)figure(2)contour(x,y,h)运行结果：地貌图及等值线图如下：0100020003000010002000300060080010001200140016001800210 50
30、0 1000 1500 2000 25000500100015002000三、本课程设计的心得体会（500 字左右）通过本次一周计算方法课程设计的训练，我进一步学习和掌握对程序的设计和编写，从中体会到了程序设计的方便和巧妙，了解到在进行编写一个程序之前，要有明确的目标和整体的设计思想，另外某些具体的细节内容也是相当的重要，这些宝贵的编程思想和从中摸索到的经验都是在编程的过程中获得的宝贵财富，这些经验对我以后的编程会有很大的帮助的，我要好好利用。虽然这次课程设计是在参考程序的基础上进行的，但是我觉得对自己是一个挑战和锻炼。我很欣慰自己能在程序中加入自己的想法和有关程序内容，也就是对它的程序改进
31、了 一 番 ， 并 有 创 新 。 但 是 我 感 觉 自 己 的 创 新 不 够 典 型 ， 总 之 还 不是 很 满 意 。 另 外 由 于 时 间 的 紧 迫 和 对 知 识 的 了 解 不 够 广 泛 ， 造 成 了 系 统 中 还 存 在许 多 不 足 ， 功能上还不够完善。以后我会继续努力，大胆创新。这次课程设计让我充分认识到了自己的不足，认识到了动手能力的重要性。我会在以后的学习中更加努力自己，提高自己，让自己写出更好更完善的程序，为以后的编程打好基础。同时，一周的课程设计周，使我对计算方法这门课程有了更多的了解和收获，在不断的思考题目和编程中，又温习了一遍老师所教22授的知识，对我大有裨益。通过这次课程设计实习，我从中学到了很多，也真正领悟到了“态度决定一切”这句话的真正含义。首先这次课程设计使我懂得了理论与实际相结合的重要性。只有理论知识是远远不够的，只有把所学的理论知识与实践相结合起来，把理论付诸实际行动，从实践中得出结论，才能更加深刻的理解书本所学知识，从而提高自己的实际动手能力和独立思考的能力。
