数值分析课程上机作业计算报告
班 级：
学 号： 姓 名：
专 业：大地测量学及测绘工程
指导老师： 联系电话：
西南交通大学-数值分析--报告
序 言
通过数值分析的理论知识的学习，此次实验将我们学过的理论知识运用于实践之中。本次实验，我选用的计算机语言为MATLAB，其主要有一下几个特点。 1．编程效率高
MATLAB是一种面向科学与工程计算的高级语言，允许使用数学形式的语言编写程序，且比BASIC、FORTRAN和C等语言更加接近我们书写计算公式的思维方式，用MATLAB编写程序犹如在演算纸上排列出公式与求解问题。因此，MATLAB语言也可通俗地称为演算纸式科学算法语言。由于它编写简单，所以编程效率高，易学易懂 2. 用户使用方便
MATLAB语言与其他语言相比，较好的解决了上述问题，把编辑、编译、链接和执行融为一体。它能在同一画面上进行灵活操作，快速排除输入程序中的书写错误、语法错误以至语义错误，从而加快了用户编写、修改和调试程序的速度，可以说在编程和调试过程中它是一种比VB还要简单的语言。 3. 方便的绘图功能
MATLAB的绘图是十分方便的，它有一系列绘图函数（命令），例如线性坐标、对数坐标、半对数坐标及极坐标，均只需调用不同的绘图函数（命令），在图上标出图题、XY轴标注，格（栅）绘制也只需调用相应的命令，简单易行。另外，在调用绘图函数时调整自变量可绘出不变颜色的点、线、复线或多重线。这种为科学研究着想的设计是通用的编程语言所不能及的。
2 / 21
目 录
1.实验一 ·························································································· 1
1.1题目 ······················································································ 1 1.2计算思路 ················································································ 1 1.3计算结果 ················································································ 1 1.4总结 ······················································································ 6 2.第二题 ·························································································· 7
2.1题目 ······················································································ 7 2.2 松弛思想分析 ········································································· 7 2.3问题的求解 ············································································· 7 2.4总结 ····················································································· 10 3.第三题 ························································································· 11
3.1题目 ····················································································· 11 3.2 Runge-Kutta法的基本思想 ························································· 11 3.3 问题的求解 ··········································································· 11 3.4问题的总结 ············································································ 14 总结 ······························································································· 15 附件 ······························································································· 16
实验一程序设计 ·········································································· 16 实验二程序设计 ·········································································· 16 实验三程序设计 ·········································································· 17
实验一：插值问题
1.1题目
已知：a=-5,b=5, 以下是某函数 f(x)的一些点（xk,yk）, 其中xk=a+0.1(k-1) ,k=1,..,101；（数据略）。
请用插值类方法给出函数f(x)的一个解决方案和具体结果。并通过实验考虑下列问题：
（1）Ln(x)的次数n越高，逼近f(x)的程度越好？ （2）高次插值收敛性如何？
（3）如何选择等距插值多项式次数 ？
（4）若要精度增高，你有什么想法？ 比如一定用插值吗？ （5）逼近某个函数不用插值方式，有何变通之举？ （6）函数之间的误差如何度量，逼近的标准又是什么？ (7) 如何比较好的使用插值多项式呢？
1.2计算思路
本题我选用拉格朗日插值函数，拉格朗日插值函数的构建算法是从101组数据中等距选取n组数据来构造n阶的拉格朗日插值函数，然后在[-5,5]区间选取m个插值点带入插值多项式计算出m组（x,y），本题中，我分别取了n=6,8,10，20四种不同的阶数，然后利用matlab绘图命令查看插值函数图像与原函数图像的逼近效果。
1
1.3 计算结果
当阶数为5阶时，插值图像与原函数图象如下所示（红色曲线为插值函数曲线，蓝色为原函数曲线）
表1（迭代数据）
i 1 26 51 71 101 Xi -5.0000 -2.5000 0 1.0000 5.0000 Yi 25.0000 6.2693 10.0000 6.2693 25.0000
由图象可以看出，除了在插值点附近的拟合效果还可以，其他的点都差距比较大，总的来说效果不是很好。
1
西南交通大学-数值分析--报告
当阶数为7阶时图象如下：
表2（迭代数据）
i 1 17 33 49 65 81 97 Xi -5.0000 -3.4000 -1.8000 -0.2000 1.4000 3.0000 4.6000 Yi 25.0000 11.5600 2.8832 5.2312 3.0219 8.9991 21.1600
由此图象可以看出，利用7阶插值来进行计算，在0值附近的拟合不是很好，但总得来说比5阶函数较好。
3 / 21
西南交通大学-数值分析--报告
当阶数为10阶时，图象如下：
表3（迭代数据）
i 1 12 23 34 ... 89 100 Xi -5.0000 -3.9000 -2.8000 -1.7000 … 3.8000 4.9000 Yi 25.0000 15.2100 7.8405 2.5554 … 14.4400 24.0100
利用10阶函数进行插值计算时，在函数的两端会出现相比于原函数偏离比较大的点，这就是runge（龙格）现象，但不是很明显。
4 / 21
西南交通大学-数值分析--报告
当阶数为20阶时，图象如下：
表4（迭代数据）
i 1 6 11 16 ... 91 96
Xi -5.0000 -4.5000 -4.0000 -3.5000 … 4.0000 4.5000 Yi 25.0000 20.2500 16.0000 12.2500 … 16.0000 20.2500
由此图像可以看出，runge（龙格）现象已经非常明显了，由于在两端的龙格现象太明显，所以此图是在去除了两端偏离比较大的点。
5 / 21
西南交通大学-数值分析--报告
1.4 总结
利用拉格朗日插值函数构造插值多项式是比较常用的插值计算的方法之一，其解决了不需要求解方程组的问题，它也有一定的缺点：
（1）从实验中可以看出，并不是拉格朗日的插值阶数越高，其插值效果越好，虽然在一些局部拟合的效果较好，但是当阶数超过7阶时就容易出现runge现象，由上图中的阶数n=20时就可以看出；
（2）但是当插值阶数太少时，插值的效果又不是很好，由上图中阶数n=5时可以看出，插值效果非常不好；
（3）由上面的图象可以分析，利用拉格朗日进行插值计算时，所选取的插值点并不是在区间里平均选取最好，根据不同的函数曲线，应当用不同的方法选点。
6 / 21
西南交通大学-数值分析--报告
实验二
2.1 题目
松弛因子对SOR法收敛速度的影响。
用SOR法求解方程组Ax=b,其中
??41??-3?????1?41???-2????-2?...??, B??? ???.?...????1?41???-2???-3?1?4?????要求程序中不存系数矩阵A,分别对不同的阶数取w=1.1, 1.2, ...,1.9进行迭代，记录近似解x(k)达到||x(k)-x(k-1)||<10-6时所用的迭代次数k，观察松弛因子对收敛速度的影响，并观察当w?0或w?2会有什么影响？
2.2 松弛思想分析
SOR法（松弛法），是由高斯-塞德尔迭代方法演变来的。高斯-塞德尔迭代方法的迭代格式为xi?k?1??xi????x，而松弛法的迭代格式为xi?kk?1??xi???w?x，通
k过大量的实验发现，ω的取值很关键，如果取了好的ω，迭代方程的收敛速度会加快，根据理论基础可以知道，当0<ω<2时松弛法才能够收敛。
2.3问题的求解
分别取不同阶数系数矩阵和不同松弛因子进行计算，结果如下图所示：
当取阶数为5时，ω=1. 1时：
7 / 21
西南交通大学-数值分析--报告
迭代次数：
迭代过程：
x = 0.8250 0.7769 x = 0.9561 0.9453 x = 0.9893 0.9868 x = 0.9974 1.0025 x = 1.0010 1.0001 x = 0.9999 1.0000 x = 1.0000 1.0000 x = 1.0000 1.0000
当阶数n=6，ω=1. 3时：
迭代次数：
0.7636 0.7600 0.9426 1.0176 1.0069 1.0005 1.0001 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
8 / 21
1.0340 1.0014 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 西南交通大学-数值分析--报告
迭代过程：
x = 0.9750 0.9669 0.9642 0.9634 0.9631 1.2880 x = 0.9967 0.9973 0.9979 0.9983 1.1041 0.9474 x = 1.0001 1.0002 1.0001 1.0344 0.9629 1.0037 x = 1.0000 1.0000 1.0111 0.9812 1.0062 1.0009 x = 1.0000 1.0036 0.9917 1.0050 1.0000 0.9997 x = 1.0012 0.9966 1.0030 0.9995 0.9997 1.0000 x = 0.9985 1.0015 0.9994 0.9999 1.0000 1.0000 x = 1.0009 0.9997 1.0000 1.0001 1.0000 1.0000 x = 0.9996 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0001 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
当阶数n=5时，ω=-1时： 请输入阶数n：5 请输入w的值：-1
num= 1, error=2.53e+000 num= 20, error=3.38e+006 num= 2, error=3.93e+000 num= 21, error=9.27e+006 num= 3, error=7.35e+000 num= 22, error=2.54e+007 num= 4, error=1.80e+001 num= 23, error=6.96e+007 num= 5, error=4.84e+001 num= 24, error=1.90e+008 num= 6, error=1.32e+002 num= 25, error=5.20e+008 num= 7, error=3.63e+002 num= 26, error=1.42e+009 num= 10; error=1.00e+003 num= 27, error=3.87e+009 num= 11, error=2.77e+003 num= 30, error=1.06e+010 num= 12, error=7.68e+003 num= 31, error=2.88e+010 num= 13, error=2.12e+004 num= 32, error=7.84e+010 num= 14, error=5.87e+004 num= 33, error=2.13e+011 num= 15, error=1.62e+005 num= 34, error=5.81e+011 num= 16, error=4.47e+005 num= 35, error=1.58e+012 num= 17, error=1.23e+006 num= 36, error=4.29e+012 x =
1.0e+012 *
-0.3019 0.8370 -1.5384 2.1285 -1.9856 由数据中可以分析得出，当ω=-1时,迭代是不收敛的。
9 / 21
西南交通大学-数值分析--报告
当阶数n=5时，ω=2.5时 请输入阶数n：5 请输入w的值：2.5
num= 1, error=1.38e+001 num= 21, error=6.34e+003 num= 2, error=2.11e+001 num= 22, error=9.29e+003 num= 3, error=2.58e+001 num= 23, error=1.25e+004 num= 4, error=2.80e+001 num= 24, error=1.83e+004 num= 5, error=4.63e+001 num= 25, error=2.82e+004 num= 6, error=6.61e+001 num= 26, error=3.74e+004 num= 7, error=8.74e+001 num= 27, error=6.20e+004 num= 10, error=1.38e+002 num= 30, error=1.18e+005 num= 11, error=1.72e+002 num= 31, error=1.96e+005 num= 12, error=2.63e+002 num= 32, error=2.94e+005 num= 13, error=4.26e+002 num= 33, error=4.09e+005 num= 14, error=8.10e+002 num= 34, error=5.63e+005 num= 15, error=1.51e+003 num= 35, error=7.87e+005 num= 16, error=2.52e+003 num= 36, error=1.24e+006 num= 17, error=3.58e+003 num= 20, error=4.46e+003 x =
1.0e+005 *
0.4613 1.5993 -1.4788 1.3274 -2.6027 由数据中可以分析得出，当ω=2.5时,迭代也是不收敛的。
2.4 总结
（1）在本题中，迭代次数和矩阵的未知数有以及ω的取值有一定的联系，不一定是未知数多迭代次数就多，也有可能在取相同ω时，未知数多迭代次数反而较少，这和SOR法的迭代原理有关。
（2）当ω≤0或ω≥2.0，SOR法不收敛，这是必要条件。
10 / 21
西南交通大学-数值分析--报告
实验三
3.1 题目
用Runge-Kutta 4阶算法对初值问题y/=-20*y，y(0)=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。 注：此方程的精确解为：y=e-20x
3.2 Runge-Kutta法的基本思想
Runge-Kutta法的基本思想是通过f (x, y)某些点函数值的适当线性组合替换Euler 法中的f(xk,yk)，可能使得方法的精确度更高。由于有了这种思想，因此标准的四阶Runge-Kutta的算法如下(h为求解步长）：
K1=hf(xk,yk)
如果yk某种方法第k步的近似值，y(xk)是其准确值，其绝对误差为δk，即有δk=y(xk)-yk。假定第k步之后的计算不再有舍入误差，只是由δk引起的扰动δm，都有︳δm︳<︳δk︳，则称此方法是绝对稳定的。数值解法的稳定性一般与步长h以及f(x,y)有关。
3.3 问题的求解
给定初值x为0，我选定迭代终点为1，,分别选取h=0.1、h=0.2进行计算，并对迭代过程和迭代图像进行分析。计算结果如下：
11 / 21
西南交通大学-数值分析--报告
当h=0.1时 请输入a的值：0 请输入b的值：1 请输入步长h：0.1
x[ 1]=0.000000 y[ x[ 2]=0.100000 y[ x[ 3]=0.200000 y[ x[ 4]=0.300000 y[ x[ 5]=0.400000 y[ x[ 6]=0.500000 y[ x[ 7]=0.600000 y[ x[ 8]=0.700000 y[ x[ 9]=0.800000 y[ 1]=1.000000 yi[ 2]=0.333333 yi[ 3]=0.111111 yi[ 4]=0.037037 yi[ 5]=0.012346 yi[ 6]=0.004115 yi[ 7]=0.001372 yi[ 8]=0.000457 yi[ 9]=0.000152 yi[ 1]= 1.000000 error[ 2]= 0.135335 error[ 3]= 0.018316 error[ 4]= 0.002479 error[ 5]= 0.000335 error[ 6]= 0.000045 error[ 7]= 0.000006 error[ 8]= 0.000001 error[ 9]= 0.000000 error[ 12 / 21
1]= 0.000000 2]= 0.197998 3]= 0.092795 4]= 0.034558 5]= 0.012010 6]= 0.004070 7]= 0.001366 8]= 0.000456 9]= 0.000152
西南交通大学-数值分析--报告
当h=0.2时
13 / 21
西南交通大学-数值分析--报告
由图象和数据中可以看出，当步长为0.2时，迭代是不收敛的。
3.4问题的总结
（1）用标准的四阶 Runge-Kutta 算法计算常微分方程时，不同的步长算法的稳定性有影响。不够稳定的步长下面的计算，误差会越来越大，结果失真严重。
（2）一般情况下，步长越小，标准的四阶 Runge-Kutta 算法的稳定性越高，精度也越高。
14 / 21
西南交通大学-数值分析--报告
总结
通过此次数值分析上机实习，认识到了计算机编程在数值分析课程中的重要性。也让我意识到了学习编程语言的重要性。要想要将数学应用于实际工程中，特别是对于工科的学生来讲，这是一门非常重要的实践课程。
在此次报告中，首次接触了Matlab这门软件，之所以选这门，是因为很多工程中对数据的处理需要使用，这对我本身的专业是非常重要的。
本学期学习的数值分析方法，让我在对数学的理解上又有了新的认识。在一连串看似杂乱无章的数据中，用数值分析进行处理，有了很多的收获。在学习的过程中，意识到了这门课对我专业的重要性，让我对自己的专业有了更好的发挥。对我写论文以及解决实际工程问题有很大的帮助。
最后，感谢老师给我这样的机会去接触这门语言，虽然只了解了皮毛，可是仍然收获颇多。由于初次接触这门软件，在报告中仍难免会有不完善甚至错误的地方，望谅解！
15 / 21
西南交通大学-数值分析--报告
附录
实验一程序
function yy=L(xi) f=fopen('xk.txt','r'); x=fscanf(f,'%f'); g=fopen('yk.txt','r'); y=fscanf(g,'%f');
n=input('请输入插值函数阶数：'); s=0;
xi=[-3.5:0.1:3.5];
for i=1:fix(101/(n-1)):101 t=ones(1,length(xi));
for j=1:fix(101/(n-1)):101 if j~=i
t=t.*(xi-x(j))/(x(i)-x(j)); end end
s=s+t*y(i); end yy=s;
plot(x,y,'b',xi,yy,'r')
16 / 21
西南交通大学-数值分析--报告
实验二程序
n=input('请输入阶数n：'); w=input('请输入w的值：'); for i=1:n for j=1:n if i==j
a(i,j)=-4;
elseif (i==j+1|i==j-1) a(i,j)=1; else
a(i,j)=0; end end end
for i=1:n
if (i==1|i==n) b(i)=-3; else
b(i)=-2; end end
x=zeros(1,n); for num=1:30 error=0; for i=1:n s=0; xb=x(i); for j=1:n
if i~=j,s=s+a(i,j)*x(j);end end
x(i)=w*(b(i)-s)/a(i,i)+(1-w)*x(i); error=error+abs(x(i)-xb); end
fprintf('num=%3.o, error=%7.2e\\n',num,error) if error/3<0.000001,break;end end x
17 / 21
西南交通大学-数值分析--报告
实验三程序
a=input('请输入a的值：'); b=input('请输入b的值：'); h=input('请输入步长h：'); x=[a:h:b]; y(1)=1;
yi(1)=exp(-20*x(1)); error(1)=0; n=(b-a)/h+1; for i=2:n
k1=-20*y(i-1)*h;
k2=-20*(y(i-1)+k1/2)*h; k3=-20*(y(i-1)+k2/2)*h; k4=-20*(y(i-1)+k3)*h;
y(i)=y(i-1)+(k1+2*k2+2*k3+k4)/6; yi(i)=exp(-20*x(i)); error(i)=y(i)-yi(i);
fprintf('x[%2.0f]=.6f y[%2.0f]=.6f error[%2.0f]=.6f\\n',i-1,x(i-1),i-1,y(i-1),i-1,yi(i-1),i-1,error(i-1)) end
plot(x,y,'g',x,yi,'r')
18 / 21
yi[%2.0f]=.6f
西南交通大学-数值分析--报告
