数值分析上机报告.docx
《数值分析上机报告.docx》由会员分享，可在线阅读，更多相关《数值分析上机报告.docx（22页珍藏版）》请在冰豆网上搜索。
数值分析上机报告
数值分析上机报告
指导教师：
姓名：
学号：
专业：
联系电话：
电子邮箱：
目录
序言1
一、插值法3
二、松弛因子对SOR法收敛速度的影响9
三、Runge-Kutta4阶算法11
总结14
附录15
习题一程序：
15
习题二程序18
习题三程序：
19
序言
随着电子计算机技术的发展，数值方法在工程技术领域中的应用越来越广泛，且已成为数学与计算机之间的桥梁。
解决工程问题，往往需要处理很多数学模型，这就要花费大量的人力和时间，但是还有不少数学模型无法用解析法得到解。
使用数值方法并利用计算机，就可以克服这些困难。
事实上，科学计算已经与理论分析、科学实验成为平行的研究和解决科技问题的科学手段，经常被科技工作者所采用。
作为科学计算的核心内容——数值分析（数值计算方法），已逐渐成为广大科技工作者必备的基本知识并越来越被人重视。
应用数值计算方法解决科学研究和工程中的实际问题，首先要建立描述具体问题的适当数学模型；其次，要选择一定的计算方法并制定相应的计算方案；并编制、设计合理的计算机程序；最后由计算机计算出数值结果。
其中，计算方案的设计和计算方法的选择是上述求解过程中极为重要的环节，是程序设计和分析数值计算结果准确性的基础。
本实验所使用的是MATLAB程序设计语言。
MATLAB语言有不同于其他高级语言的特点，它被称为第四代计算机语言。
正如第三代计算机语言如FORTRAN语言与C语言等使人们摆脱了对计算机硬件的操作一样，MATLAB语言使人们从繁琐的程序代码中解放出来。
它的丰富的函数使开发者无需重复编程，只要简单地调用和使用即可。
MATLAB语言最大的特点是简单和直接。
MATLAB语言的主要特点有：
1．编程效率高；2．用户使用方便；3．扩充能力强，交互性好；4．移植性好，开放性好；5．语句简单，内涵丰富；6．高效方便的矩阵和数组运算；7．方便的绘图功能。
结合此次上机实习题目和Matlab工具的内置函数丰富，图形功能强大，运算符丰富等特点，选取Matlab语言来完成此次数值分析上机作业。
此次上机作业分为三个大题：
（1）对数据进行插值计算，了解不同插值方法的注意点，并比较几种插值法的优点和缺点；
（2）理解松弛因子对SOR法收敛速度的影响；
（3）使用Runge-Kutta4阶算法对初值问题按不同步长求解，观察稳定区间的作用。
一、插值法
1.习题题目：
已知：
a=-5,b=5,以下是某函数f（x）的一些点（xk,yk）,其中xk=a+0.1（k-1）,k=1,..,101，xk=a+0.1k,
请用插值类方法给出函数f（x）的一个解决方案和具体结果。
并通过实验考虑下列问题
（1）Ln（x）的次数n越高，逼近f（x）的程度越好？
（2）高次插值收敛性如何？
（3）如何选择等距插值多项式次数？
（4）若要精度增高，你有什么想法？
比如一定用插值吗？
（5）逼近某个函数不用插值方式，有何变通之举？
（6）函数之间的误差如何度量，逼近的标准又是什么？
（7）如何比较好的使用插值多项式呢？
xk=
Columns1through7
-5.0000-4.9000-4.8000-4.7000-4.6000-4.5000-4.4000
Columns8through14
-4.3000-4.2000-4.1000-4.0000-3.9000-3.8000-3.7000
Columns15through21
-3.6000-3.5000-3.4000-3.3000-3.2000-3.1000-3.0000
Columns22through28
-2.9000-2.8000-2.7000-2.6000-2.5000-2.4000-2.3000
Columns29through35
-2.2000-2.1000-2.0000-1.9000-1.8000-1.7000-1.6000
Columns36through42
-1.5000-1.4000-1.3000-1.2000-1.1000-1.0000-0.9000
Columns43through49
-0.8000-0.7000-0.6000-0.5000-0.4000-0.3000-0.2000
Columns50through56
-0.100000.10000.20000.30000.40000.5000
Columns57through63
0.60000.70000.80000.90001.00001.10001.2000
Columns64through70
1.30001.40001.50001.60001.70001.80001.9000
Columns71through77
2.00002.10002.20002.30002.40002.50002.6000
Columns78through84
2.70002.80002.90003.00003.10003.20003.3000
Columns85through91
3.40003.50003.60003.70003.80003.90004.0000
Columns92through98
4.10004.20004.30004.40004.50004.60004.7000
Columns99through101
4.80004.90005.0000
y（xk）=yk=
Columns1through7
25.000024.010023.040022.090021.160020.250019.3600
Columns8through14
18.490017.640016.810016.000015.210014.440013.6900
Columns15through21
12.960012.250011.560010.889910.23979.60938.9991
Columns22through28
8.40927.84057.29416.77056.26935.78665.3144
Columns29through35
4.84034.35223.84633.34022.88322.55542.4475
Columns36through42
2.61543.02193.49203.71493.32322.0435-0.1277
Columns43through49
-2.8066-5.2470-6.5469-5.9893-3.38620.73655.2312
Columns50through56
8.698510.00008.69855.23120.7365-3.3862-5.9893
Columns57through63
-6.5469-5.2470-2.8066-0.12772.04353.32323.7149
Columns64through70
3.49203.02192.61542.44752.55542.88323.3402
Columns71through77
3.84634.35224.84035.31445.78666.26936.7705
Columns78through84
7.29417.84058.40928.99919.609310.239710.8899
Columns85through91
11.560012.250012.960013.690014.440015.210016.0000
Columns92through98
16.810017.640018.490019.360020.250021.160022.0900
Columns99through101
23.040024.010025.0000
2.计算及结果：
n次Lagrange插值多项式为：
首先，利用MATLAB描绘出所有点的散点图，如图1.1。
图1.1101个插值点的散点图
分别在[-5,5]区间上依次按等间距选取51个点、21个点、11个点及5个点，进行Lagrange的50次、20次、10次和4次插值，并分别绘制不同次数插值函数的函数曲线，对不同次数的插值结果进行比较：
图1.250次插值多项式曲线图1.320次插值多项式曲线
图1.410次插值多项式曲线图1.54次插值多项式曲线
从这几个图可以看出，高次插值容易出现不收敛现象，即龙格现象，而低次插值很显然又不能很好地反映被插值函数。
基于这种情况，可以考虑采用分段插值。
从给定数据点的散点图看出，可以将[-5,5]上的离散点分成7段分别进行拉格朗日二次插值，这7个区间分别为[-5,-1.4]，[-1.4,-0.8]，[-0.8,-0.3]，[-0.3,0.3]，[0.3,0.8]，[0.8,1.4]，[1.4,5]，每个插值区间选取3个点进行拉格朗日二次插值。
分段拉格朗日插值结果如图：
1.6拉格朗日2次分段插值与数据散点图
插值多项式为：
从计算结果可以看出：
（1）对于拉格朗日插值法来说，并不是插值次数n越高，f（x）的逼近程度越好。
当插值次数n较大时，会出现龙格（Runge）现象。
（2）对于高次插值多项式来说，会在某一区间内发散，而同时也会在某一区间内收敛性很高。
拉格朗日插值法得到插值多项式与原函数的误差在区间内分布不均匀，而且增加节点并不能保证两节点之间的n次插值函数能很好地逼近原函数f（x），而且高次插值的误差积累会更加严重，很难达到一致收敛。
（3）可以通过截断误差来估计等距插值多项式的次数。
（4）除了插值法以外，还可以用曲线拟合方法来逼近和最小二乘法。
（5）逼近函数的方法有好多种，可以用最佳平方逼近，最佳一致逼近。
插值多项式的值在节点上是与函数值相同的，而拟合可以在多个节点条件下作较低次的近似多项式，节点上的误差符合最小二乘法原理。
最佳平方逼近的误差在欧氏范数的意义下尽可能小；最佳一致逼近的误差在无穷范数的意义下尽可能小。
（6）为了衡量函数g（x）对f（x）的逼近程度，引入以下几种度量：
a.一致逼近度：
定义在[a,b]上的连续函数f（x）,g（x）之间的接近程度规定如下：
b.平方逼近度：
定义在[a,b]上的平方可积函数f（x）,g（x）之间的接近程度规定如下：
c.P次幂逼近度：
定义在[a,b]上的P次幂可积函数f（x）,g（x）之间的接近程度规定如下：
（7）对于本题来说，可以采用分段二次插值方法，因为本题的节点有101个，如果直接采用拉格朗日插值法会因为节点太多次数抬高而产生龙格现象，分段插值则会因为插值点越多，插值效果越好。
二、松弛因子对SOR法收敛速度的影响
1.习题题目：
用SOR法求解方程组Ax=b,其中
要求程序中不存系数矩阵A,分别对不同的阶数取w=1.1,1.2,...,1.9进行迭代，记录近似解x（k）达到||x（k）-x（k-1）||<10-6时所用的迭代次数k，观察松弛因子对收敛速度的影响，并观察当w0或w2会有什么影响？
2．计算及结果：
超松弛法SOR的基本思想就是在Gauss-Seidel迭代格式修正项上乘一个因子来加快收敛的速度。
Gauss-Seidel方法迭代的矩阵形式为：
为了研究收敛性，可以将Gauss-Seidel迭代格式改写为：
那么添加松弛因子
后的迭代格式为：
最后解出SOR超松弛法的迭代格式：
根据SOR法的迭代格式可以建立函数[x,k]=SOR（n,w）（见附录），分析题目很容易得到该方程组的解为
，因此只要比较方程组阶数n和松弛因子
对收敛次数k的影响。
分别取n=5、10、15、20，再分别对每个n取不同松弛因子1.1，1.2，…，1.9，计算收敛次数k，其结果见表2-1.
ω
n
1.1
1.2
1.3
1.4
1.5
1.6
1.7
1.8
1.9
5
8
11
14
18
24
31
43
67
142
10
12
15
18
20
25
32
46
70
146
15
12
15
19
24
28
34
45
70
148
20
12
15
19
24
31
38
47
77
150
表2-1n和
对收敛速度的影响
由表2-1可得，n越大即方程组阶数越高，松弛因子
在[1.1,1.9]区间内取值越大，迭代次数增大，收敛速度降低，迭代运算量增加。
当n=5和n=10时，分别取
=-1和
=3，利用MATLAB计算发现迭代次数超过500，可能不收敛。
所以当松弛因子
或
时，SOR迭代不收敛。
三、Runge-Kutta4阶算法
1.习题题目：
用Runge-Kutta4阶算法对初值问题y/=-20*y，y（0）=1按不同步长求解，用于观察稳定区间的作用，推荐两种步长h=0.1,0.2。
注：
此方程的精确解为：
y=e-20x
Runge-Kutta4阶算法的公式为：
考虑该算法的绝对稳定区间，引入试验方程
，当步长为h时，将Runge-Kutta4阶算法应用于上式，并假设产生误差
，若此误差以后逐步衰减，就称该算法相对于λh“绝对稳定”，λh的全体构成“绝对稳定区间”。
Runge-Kutta4阶算法是收敛的，绝对稳定区间为（-2.78，0）。
（1）步长为0.1时，y（xk）表示真实值，yk表示为Runge-Kutta法得到的近似值，εk表示为误差。
计算结果如表3-1和图3.1所示：
表3-1h=0.1时的计算结果
xk
近似解yk
精确解y（xk）
误差εk
x（0）=0
1.0000
1.0000
0
x
（1）=0.1
0.3333
0.1353
0.1980
x
（2）=0.2
0.1111
0.0183
0.0928
x（3）=0.3
0.0370
0.0025
0.0346
x（4）=0.4
0.0123
0.0003
0.0120
x（5）=0.5
0.0041
0.0000
0.0041
x（6）=0.6
0.0014
0.0000
0.0014
x（7）=0.7
0.0005
0.0000
0.0005
x（8）=0.8
0.0002
0.0000
0.0002
x（9）=0.9
0.0001
0.0000
0.0001
x（10）=1.0
0.0000
0.0000
0.0000
图3.1h=0.1近似解与精确解的比较
由上述结果可以看到误差|εk|是逐渐缩小的，说明当h=0.1时Runge-Kutta法是收敛的。
（2）步长h=0.2时，计算结果如表3-2和图3.2所示：
表3-2h=0.2时的计算结果
xk
近似解yk
精确解y（xk）
误差εk
x（0）=0
1
1.0000
0.0000
x
（1）=0.2
5
0.0183
4.9817
x
（2）=0.4
25
0.0003
24.9997
x（3）=0.6
125
0.0000
125
x（4）=0.8
625
0.0000
625
x（5）=1.0
3125
0.0000
3125
x（6）=1.2
15625
0.0000
15625
x（7）=1.4
78125
0.0000
78125
x（8）=1.6
390625
0.0000
390625
x（9）=1.8
1953125
0.0000
1953125
x（10）=2.0
9765625
0.0000
9765625
图3.2h=0.2近似解与精确解的比较
由计算结果可以看出，当h=0.2时，随着x的增大，近似解与精确解之间的误差εk越来越大，说明此时Runge-Kutta法不收敛。
结论分析：
四阶Runge-Kutta法的绝对稳定区间是（-2.78，0），即当
时，四阶Runge-Kutta法才收敛，本题中，h=0.1时
，满足收敛条件。
当h=0.2时，
，不满足收敛条件。
所以只有当步长满足一定条件时，Runge-Kutta法才稳定收敛，并不是步长越小越好。
总结
通过完成这一段时间的上机报告，我对数值分析有了更深的认识,虽然以前也会用到Matlab编程,但通过学习数值分析,我了解到了数值分析的原理,以及进行数值分析的原因,认识到数值分析后面的数学原理.
这次上机报告让我对这一学期学到的数值分析知识进行了一次汇总,了解到自己还有什么理解不透彻的地方,并加以解决,之前一直在课堂上听讲,这次自己编写程序,加深了对数值分析的理解.
在这学期的数值分析课程中,谢谢赵老师的精心授课,帮助学生们较深刻地理解了数值分析这门课.
附录
习题一程序：
（1）n次拉格朗日插值子函数：
functionf=lagrange_n（n）%n次拉格朗日插值子函数，n为拉格朗日插值次数x=[-5:
0.1:
5];
y=[2524.0123.0422.0921.1620.2519.3618.4917.6416.81...
1615.2114.4413.6912.9612.2511.5610.889910.23979.6093...
8.99918.40927.84057.29416.77056.26935.78665.31444.84034.3522...
3.84633.34022.88322.55542.44752.61543.02193.49203.71493.3232...
2.0435-0.1277-2.8066-5.247-6.5469-5.9893-3.38620.73655.23128.6985...
108.69855.23120.7365-3.3862-5.9893-6.5469-5.247-2.8066-0.1277...
2.04353.32323.71493.4923.02192.61542.44752.55542.88323.3402...
3.84634.35224.84035.31445.78666.26936.77057.29417.84058.4092...
8.99919.609310.239710.889911.5612.2512.9613.6914.4415.21...
1616.8117.6418.4919.3620.2521.1622.0923.0424.0125];
h=100/n;
x1=x（1:
h:
101）;
y1=y（1:
h:
101）;
s=0;
fori=1:
n+1
z=ones（1,101）;
forj=1:
n+1
ifj~=i
z=z.*（x-x1（j））/（x1（i）-x1（j））;%拉格朗日插值基函数
end
end
s=s+z*y1（i）;
end
f=s;
plot（x,y,'r.',x,f）;%绘制原数据散点图和插值函数曲线
legend（'y','lan'）
（2）分段拉格朗日2次插值子程序：
functionf=sublagrange（x,y,x0）%拉格朗日插值子函数
m=length（x）;
n=length（y）;
ifm~=n
error（'向量x和y的长度不一致'）
end
f=0;
fori=1:
n
z=ones（1,length（x0））;
forj=1:
n
ifj~=i
z=z.*（x0-x（j））/（x（i）-x（j））;%拉格朗日插值基函数
end
end
f=f+z*y（i）;
end
%以下是利用拉格朗日插值子函数进行分段2次插值的m文件程序
x=[-5:
0.1:
5];
y=[2524.0123.0422.0921.1620.2519.3618.4917.6416.81...
1615.2114.4413.6912.9612.2511.5610.889910.23979.6093...
8.99918.40927.84057.29416.77056.26935.78665.31444.84034.3522...
3.84633.34022.88322.55542.44752.61543.02193.49203.71493.3232...
2.0435-0.1277-2.8066-5.247-6.5469-5.9893-3.38620.73655.23128.6985...
108.69855.23120.7365-3.3862-5.9893-6.5469-5.247-2.8066-0.1277...
2.04353.32323.71493.4923.02192.61542.44752.55542.88323.3402...
3.84634.35224.84035.31445.78666.26936.77057.29417.84058.4092...
8.99919.609310.239710.889911.5612.2512.9613.6914.4415.21...
1616.8117.6418.4919.3620.2521.1622.0923.0424.0125];
x1=[-5,-1.6,-1.4];%第一个插值区间
y1=[25,2.4475,3.0219];
x0=[-5:
0.1:
-1.4];
f1=sublagrange（x1,y1,x0）;%利用预先编写的子函数进行拉格朗日插值
po1=polyfit（x1,y1,2）;%求二次多项式系数
px1=poly2str（po1,'x'）%写出二次多项式表达式
plot（x1,y1,'ro',x0,f1,'r',x,y,'b+'）;%绘制第一个插值区间的插值多项式曲线
holdon;
x2=[-1.4,-1.2,-0.8];%第二个插值区间
y2=[3.0219,3.7149,-2.8066];
x0=[
升级会员 