完整word版数值分析课程设计含代码.docx
《完整word版数值分析课程设计含代码.docx》由会员分享，可在线阅读，更多相关《完整word版数值分析课程设计含代码.docx（28页珍藏版）》请在冰点文库上搜索。
完整word版数值分析课程设计含代码
成绩评定表
学生姓名
班级学号
专业
信息与计算科学
课程设计题目
数值分析算法案例
评
语
组长签字：
成绩
日期
20年月日
课程设计任务书
学院
理学院
专业
信息与计算科学
学生姓名
班级学号
课程设计题目
数值分析算法案例
实践教学要求与任务:
要求：
格式以学校毕业论文格式要求为准，不准粘贴图片，尤其公式。
对每个试验，要求有：
实验基本原理，实验目的，实验内容及数据来源和实验结论。
以班级为单位统一装订封皮。
6月25日，十八周的周二交论文
每人至少四个实验，最少15页
任务（实验项目）：
线性方程组数值解法参考题目：
（1）列主元Gauss消去法
（2）LU分解法
插值法和数据拟合参考题目：
（1）Lagrange插值
（2）Newton插值（3）最小二乘拟合
数值积分参考题目：
（1）复化Simpon积分
（2）变步长的梯形积分公式（3）龙贝格求积公式
常微分方程数值解Runge-Kutta方法
数值方法实际应用用数值方法解决实际问题（自选）
工作计划与进度安排:
线性方程组数值解法（4学时）
插值法和数据拟合（4学时）
数值积分（4学时）
常微分方程数值解（4学时）
数值方法实际应用（4学时）
答辩（4学时）
指导教师：
201年月日
专业负责人：
201年月日
学院教学副院长：
201年月日
摘要
实验方法与理论方法是推动科学技术发展的两大基本方法，但有局限性。
许多研究对象，由于空间或时间的限制，既不可能用理论精确描述，也不能用实验手段实现。
数值模拟或称为科学计算突破了实验和理论科学的局限，在科技发展中起到越来越重要的作用。
可以认为，科学计算已于实验、理论一起成为科学方法上不可或缺的三个主要手段。
计算数学的研究是科学计算的主要组成部分，而数值分析则是计算数学的核心。
数值计算是研究使用计算机来解决各种数学问题的近似计算方法与理论，其任务是提供在计算机上可解的、理论可靠的、计算复杂性低的各种常用算法。
数值分析的主要内容：
1）、数值代数：
求解线性和非线性方程组的解，分直接方法和间接方法两大类；
2）、插值、曲线拟合和数值逼近；
3）、数值微分和数值积分；
4）、常微分和偏微分方程数值解法。
本文主要通过Matlab软件，对数值分析中的一些问题进行求解，如列主元Gauss消去法，Lagrange插值多项式，复化Simpson公式，Runge-Kutta方法以及数值分析在实际问题中的应用，并在求解的过程中更加熟识这门课程的主要内容，以及加强对课程知识的掌握。
在学习与设计计算方法时，从数学理论角度，学会分析方法的误差、收敛性和稳定性，保证计算方法的准确性；从实际应用的角度出发，掌握计算方法的结构与流程，能够把计算方法转换为可在计算机上直接处理的程序，保证算法的可用性。
关键词：
列主元Gauss消去法；Lagrange插值；复化Simpson公式；Runge-Kutta
实验一列主元Gauss消去法
1.1实验目的
1）理解列主元消去法的原理；
2）熟悉列主元消去法的计算步骤，能用代码编写；
3）解决实际问题。
1.2基本原理
在顺序Gauss消去法中，必须要求
；否则无法进行计算。
即使
，但其绝对值
很小，由于舍入误差的影响，也可能会引起很大的误差，从而使上述方法失效。
为了使消元过程中减小舍入误差和不至于中断，可以按照不同的自然顺序进行消元。
在第k步消元时，增广矩阵为
（1.1）
不一定选取
作为主元，而从同列
中选取绝对值最大的作为主元素，即
（1.2）
若
，此时矩阵不可逆，方程的解不确定，则停止计算；否则，当r>k时，则其增广矩阵中交换第k行和第r行，即
（1.3）
使
成为主元。
然后再按Gauss消去法进行消元运算。
于是就得到列主元Gauss消去法。
1.3实验内容
1.3.1程序来源
首先建立一个gaussMethod.m的文件，用来实现列主元的消去方法。
文件内容如下：
functionx=gaussMethod（A,b）
%高斯列主元消去法，要求系数矩阵非奇异的
n=size（A,1）;
ifabs（det（A））<=1e-8
error（'系数矩阵是奇异的'）;return;
end
fork=1:
n
ak=max（abs（A（k:
n,k）））;
index=find（A（:
k）==ak）;
iflength（index）==0
index=find（A（:
k）==-ak）;
end
%交换列主元
temp=A（index,:
）;
A（index,:
）=A（k,:
）;
A（k,:
）=temp;
temp=b（index）;b（index）=b（k）;b（k）=temp;%消元过程
fori=k+1:
n
m=A（i,k）/A（k,k）;%消除列元素
A（i,k+1:
n）=A（i,k+1:
n）-m*A（k,k+1:
n）;
b（i）=b（i）-m*b（k）;
end
end%回代过程
x（n）=b（n）/A（n,n）;
fork=n-1:
-1:
1;
x（k）=（b（k）-A（k,k+1:
n）*x（k+1:
n）'）/A（k,k）;
end;end
然后调用gaussMethod函数，来实现列主元的高斯消去法。
建立一个文件gauss，内容如下：
clear
disp（'**********************************************'）
x=gaussMethod（input（'请输入系数矩阵：
'）,input（'请输入常数列：
'））
disp（'**********************************************'）
1.3.2实例分析
例：
在Matlab上，利用列主元法求线性方程组的解:
解：
运行程序，按照提示输入方程的系数矩阵及常数列，如下所示：
**********************************************
请输入系数矩阵：
[1214;2043;4221;-3132]
请输入常数列：
[13;28;20;6]
x=
3-142
**********************************************
即该方程的解为：
1.4实验结论
把向量计算得到的解带入方程组，验证正确性，和其他的方法比较，列主元具有一定的简单性，比较容易实现。
避免使用其他方法的误差或不能进行性。
而列主元也有一定的限制，要求行列式的值不为0。
实验二拉格朗日插值多项式
2.1实验目的
1）熟悉简单的拉格朗日插值多项式的基本概念；
2）熟悉Lagrange公式及源代码，会利用它来计算基本函数；
3）能构造出正确的插值多项式；
2.2基本原理
设函数
在区间[a,b]上有定义，且已知在点
上的值
若存在一个次数不超过n的多项式
（2.1）
使其满足
（2.2）
则称
为
的n次插值多项式，称点
为插值节点，称条件（2.2）为插值条件。
包含插值节点的区间成为插值区间。
通过平面上不同的两点可以确定一条直线经过这两点，就是拉格朗日线性插值问题，对于不在同一直线的三点得到的插值多项式为抛物线。
拉格朗日是比较基础的方法，本身比较容易实现，容易理解。
给定n+1个不同节点，构造
的n次拉格朗日插值多项式：
，
（2.3）
2.3实验内容
2.3.1程序来源
首先建立一个Lagrange.m的文件，用来实现Lagrange插值。
文件内容如下：
%输入：
x是插值节点横坐标向量；y是插值节点对应纵坐标向量
%输出：
C是拉格朗日插值多项式的系数矩阵；L是插值函数系数矩阵
function[C,L]=Lagrange（x,y）
w=length（x）;
n=w-1;
L=zeros（w,w）;
fork=1:
n+1
V=1;
forj=1:
n+1
ifk~=j
V=conv（V,poly（x（j）））/（x（k）-x（j））;
end
end
L（k,:
）=V;
end
C=y*L
然后调用Lagrange函数，来实现Lagrange插值法。
建立一个文件Lg，内容如下：
clear
disp（'**********************************************'）
x=input（'请输入已知点的横坐标组：
'）;
y=input（'请输入已知点的纵坐标：
'）;
[C,L]=Lagrange（x,y）;
yi=polyval（C,input（'请输入需要计算得横坐标组：
'））
xx=1.5:
0.05:
6.5;
yy=polyval（C,xx）;
plot（xx,yy,x,y,'o'）
disp（'**********************************************'）
2.3.2实例分析
例1有4对数据（1.6,3.3），（2.7,4.22），（3.9,5.61），（5.6,2.94），写出这4个数据点的Lagrange插值公式，并计算出横坐标组xi=[2.101,4.234]时对应的纵坐标值。
解：
4个数据点的Lagrange插值公式为：
运行程序，按照提示输入已知点的横坐标组、纵坐标组及需要计算得横坐标组，如下所示：
**********************************************
请输入已知点的横坐标组：
[1.6,2.7,3.9,5.6]
请输入已知点的纵坐标：
[3.3,1.22,5.61,2.94]
C=
-1.053911.0551-34.493334.5053
请输入需要计算得横坐标组：
[2.101,4.234]
yi=
1.05966.6457
**********************************************
即
输出图形：
图2.1输出拟合曲线
例2将区间[-5，5]等分5份、10份，求函数
拉格朗日差值多项式，做出函数原图像，观察龙格现象。
解：
首先将区间五等分，取各端点坐标拟合曲线，输入：
clear
x=-5:
2:
5;
y=1./（1+x.^2）;
[C,L]=Lagrange（x,y）;
输出结果：
C=
0.00000.0019-0.0000-0.06920.00000.5673
即
然后输出图形，比较拟合效果。
在matlab中输入：
xx=-5:
0.1:
5;
yy=polyval（C,xx）;
holdon
plot（xx,yy,'--',x,y,'o'）
xp=-5:
0.01:
5;
z=1./（1+xp.^2）;
plot（xp,z,'r'）
legend（'拉格朗日插值曲线','插值点','原曲线'）
输出五等分插值龙格现象图形：
图2.2五等分插值龙格现象图形
十等分的过程与五等分基本相似，如下输入：
clear,clf
x=-5:
1:
5;
y=1./（1+x.^2）;
[C,L]=Lagrange（x,y）;
输出：
C=
-0.0000-0.00000.0013-0.0000-0.02440.00000.1974-0.0000-0.6742-0.00001.0000
即
然后输出图形，比较拟合效果。
在matlab中输入：
xx=-5:
0.1:
5;
yy=polyval（C,xx）;
holdon
plot（xx,yy,'--',x,y,'o'）
xp=-5:
0.01:
5;
z=1./（1+xp.^2）;
plot（xp,z,'r'）
legend（'拉格朗日插值曲线','插值点','原曲线'）
输出十等分插值图形：
图2.3十等分插值龙格现象图形
2.4实验结论
通过图像可以得出：
1）并不是插值节点越多，插值多项式逼近函数效果就越好；
2）误差较大的地方，是在插值区间两端点附近出现；
3）求插值多项式，不需要求解线性方程组，当数据比较多时，此公式才能
显示出优越性；
4）函数值可以用符号形式表示，数据点未确定的纵坐标可以用多项式表示；
实验三复化Simpson求积公式
3.1实验目的
1）了解单独的Simpson公式原理及几何意义；
2）学会复化Simpson求积公式的编程与应用；
3）掌握Matlab提供的计算积分的各种函数的使用方法；
3.2基本原理
复化Simpson公式是一种比较实用的积分方法，可以给出误差估计。
首先将区间[a,b]N等分，子区间的长度为
（3.1）
在每个子区间上采用Simpson公式，在用Simpson公式时，还要将子区间再二等分，因此有2N+1个分点。
即
（3.2）
经推导得到，
（3.3）
称为
为复化Simpson值，称（3.3）式为复化Simpson公式。
3.3实验内容
3.3.1程序来源
编写复化Simpson求积函数（函数名：
s_quad.m）
FunctionI=S_quad（x,y）
%复化求积公式
%x为被积函数自变量的等距节点；y为被积函数在节点处的函数值。
n=length（x）;
m=length（y）;%积分自变量的节点数应与它的函数值个数相同；
ifn~=m
error（'ThelengthofXandYmustbeequal'）;
return;
end
ifrem（n-1,2）~=0%如果n-1不能被2整除，则调用复化公式
error（'节点数不满足要求'）;
return;
end
N=（n-1）/2;
h=（x（n）-x
（1））/N;
a=zeros（1,n）;
fork=1:
N
a（2*k-1）=a（2*k-1）+1;
a（2*k）=a（2*k）+4;
a（2*k+1）=a（2*k+1）+1;
end
I=h/6*sum（a.*y）;
然后调用s_quad函数，来实现复化Simpson公式法。
建立一个文件SPS，内容如下：
clear
disp（'**********************************************'）
x=input（'请输入积分上下限及点间的间隔（例如-1:
0.1:
1）：
'）;
y=input（'请输入被积公式：
y='）;
I=S_quad（x,y）;
disp（'得出积分值I='）
disp（I）;
disp（'**********************************************'）
3.3.2实例分析：
例1用复化Simpson公式求积分
，在积分区间中点与点之间的间隔取为
0.1。
解：
运行程序，按照提示输入积分上下限、点间的间隔及被积公式，如下所示：
**********************************************
请输入积分上下限及点间的间隔（例如-1:
0.1:
1）：
-1:
0.1:
1
请输入被积公式：
y=exp（-x.^2）
得出积分值I=
1.4936
**********************************************
真值为：
1.4937
例2计算积分
，将区间8等分。
解：
运行程序，按照提示输入积分上下限、等分后的区间长度及被积公式，如下
所示：
**********************************************
请输入积分上下限及点间的间隔（例如-1:
0.1:
1）：
0:
0.125:
1
请输入被积公式：
y=x./（4+x.^2）
得出积分值I=
0.1116
**********************************************
真值为：
0.111572
3.4实验结论
复化Simpson计算所得的结果误差较小，精度较高，更适合科学计算与应用，且公式具有收敛性，稳定性良好。
实验四龙格-库塔（Runge-Kutta）方法
4.1实验目的
1）熟悉Runge-Kutta常微分方程初值问题的基本原理
2）了解Runge-Kutta常微分方程初值问题的计算
3）能编程实现Runge-Kutta常微分方程初值问题
4.2基本原理
龙格-库塔（Runge-Kutta）方法是一种在工程上应用广泛的高精度单步算法。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
该算法是构建在数学支持的基础之上的。
对于一阶精度的欧拉公式有：
（4.1）
当用点
处的斜率近似值
与右端点
处的斜率
的算术平均值作为平均斜率
的近似值，那么就会得到二阶精度的改进欧拉公式：
（4.2）
依次类推，如果在区间
内多预估几个点上的斜率值，并用他们的加权平均数作为平均斜率
的近似值，显然能构造出具有很高精度的高阶计算公式。
经数学推导、求解，可以得四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法。
四阶龙格库塔方法基本公式：
（4.3）
4.3实验内容
4.3.1程序来源
写出四阶Runge-Kutta方法源代码的
得.m文件：
functionR=rk4（f,a,b,ya,M）
%a,b分别为左右端点，ya为初始值，M为步长幅数。
h=（b-a）/M;
T=zeros（1,M+1）;
Y=zeros（1,M+1）;
T=a:
h:
b;
Y
（1）=ya;
forj=1:
M
k1=feval（f,T（j）,Y（j））;
k2=feval（f,T（j）+h/2,Y（j）+k1*h/2）;
k3=feval（f,T（j）+h/2,Y（j）+k2*h/2）;
k4=feval（f,T（j）+h,Y（j）+h*k3）;
Y（j+1）=Y（j）+（k1+2*k2+2*k3+k4）/6;
end
R=[T'Y'];
4.3.2实例分析：
例用四阶龙格库塔方法求解初值问题
取步长h=0.1。
并将a,b,ya,M分别取定0,1,0,10。
解：
首先写出执行函数的tt.m文件，其中tt表示f，即为求导函数：
functionf=tt（t,y）
f=t^2+t-y;
再输入调用四阶Runge-Kutta方法，执行语句如下：
>>R=rk4（'tt',0,1,0,10）
输出：
R=
00
0.100000000000000.05162708333333
0.200000000000000.16847829348958
0.300000000000000.30751722078089
0.400000000000000.46666197888861
0.500000000000000.64581185656207
0.600000000000000.84496198189452
0.700000000000001.06411211920748
0.800000000000001.30326225710000
0.900000000000001.56241239502055
1.000000000000001.84156253294245
4.4实验结论
龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算在一些点上的值，如四阶龙格-库塔方法每计算一步需要计算四次的值，这给实际计算带来一定的复杂性。
实验五数值方法实际应用
5.1实验目的
1）了解梯形公式原理及几何意义；
2）学会复化梯形求积公式的编程与应用；
3）掌握Matlab提供的计算积分的各种函数的使用方法；
5.2基本原理
在实际问题中，往往会遇到一些困难。
有些函数找不到用初等函数表示的原函数，例如，对于积分
（5.1）
而言，不存在用初等函数表示的原函数。
而有些函数虽然能找到原函数，但计算过于复杂，例如，椭圆型积分
（5.2）
而有些情况下，只能知道某些点处的函数值，并没有函数的具体表达式。
这些情况，使我们有必要研究积分的数值计算问题。
下面我们就以梯形公式为例做以说明。
所谓梯形求积公式就是用梯形面积来近似曲边梯形面积，利用梯形公式和连续增加[a,b]的区间数来逼近：
（5.3）
第j次循环在
个等距节点处对
采样。
5.3实验内容
卫星轨道是一个椭圆，椭圆周长计算公式是
，
这里a是椭圆半长轴，c是地球中心与轨道中心（椭圆中心）的距离，记h为近地点距离，H为远地点距离，R=637km为地球半径，则
我国第一颗人造卫星近地点距离h=439km，远地点距离H=2384，试求卫星轨道的周长。
解:
第一步：
先利用Matlab画出被积函数的图形。
输入程序如下：
clear
H=2384;
h=439;
R=6371;
a=（2*R+H+h）/2
c=（H-h）/2
x=0:
0.1:
pi/2;
y=sqrt（1-（c/a）^2*（sin（x））.^2）;
plot（x,y,'--'）
title（'梯形法则'）;
xlabel（'x'）;
ylabel（'y'）;
输出结果：
a=
7.782500000000000e+003
c=
9.725000000000000e+002
输出图形：
图5.1被积函数的图形
第二步：
应用数值积分梯形公式。
首先建立一个名为trapezg.m的M文件，程序如下：
functionI=trapezg（f_name3,a,b,n）
formatlong
%输出用15位数字表示
n=n;
h=（b-a）/n;
x=a+（0:
n）*h;
f=feval（f_name3,x）;
I=h/2*（f
（1）+f（n+1））;
ifn>1I=I+h*sum（f（2:
n））;
end
h1=（b-a）/100;
xc=a+（0:
100）*h1;
fc=feval（f_name3,xc）;
plot（xc,fc,'r'）;
holdon
xlabel（'x'）;
ylabel（'y'）;
plot（x,f）
title（'数值积分梯形效果图'）;
plot（x,zeros（size（x））,'.'）
fori=1:
n;plot（[x（i）,x（i）],[0,f（i）]）,
end
然后建立一个名为f_name3.m的M文件定义函数，Matlab命令如下：
functiony=f_name3（x）
y=sqrt（1-（9.725000e+002/7.782500e+003）^2*（sin（x））.^2）-0
升级会员 