数值分析课程设计Word下载.docx
《数值分析课程设计Word下载.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计Word下载.docx（15页珍藏版）》请在冰豆网上搜索。
2基本原理12
3实验内容及数据来源12
4实验结论14
心得15
参考文献16
实验一列主元Gauss消去法
1实验目的与要求
熟悉列主元Gauss消去法的基本原理；
会用matlab程序实现列主元Gauss消去法的全过程。
2基本原理
为了提高计算的数值稳定性,在消元过程中采用选择主元的方法.常采用的是列主元消去法。
给定线性方程组Ax=b,记A
（1）=A,b
（1）=b,列主元Gauss消去法的具体过程如下:
首先在增广矩阵B
（1）=（A
（1）,b
（1））的第一列元素中,取
为主元素，
。
然后进行第一步消元得增广矩阵B
（2）=（A
（2）,b
（2））。
再在矩阵B
（2）=（A
（2）,b
（2））的第二列元素中,取
然后进行第二步消元得增广矩阵B（3）=（A（3），b（3））。
按此方法继续进行下去，经过n-1步选主元和消元运算，得到增广矩阵B（n）=（A（n）,b（n））。
则方程组A（n）x=b（n）是与原方程组等价的上三角形方程组，可进行回代求解。
易证,只要|A|
0，列主元Gauss消去法就可顺利进行。
可见，列主元Gauss消去法是在每一步消元前,在主元所在的一列选取绝对值最大的元素作为主元素.而全主元Gauss消去法是在每一步消元前,在所有元素中选取绝对值最大的元素作为主元素.但由于运算量大增,实际应用中并不经常使用.
3实验内容及数据来源
下面先用Matlab编写列主元Gauss消去的函数文件如下：
%列选主元的高斯消去法
functionX=lufact_my（A,B）
%InpiutA是系数矩阵，B是右端项
%Outputx是解
[N,N]=size（A）;
X=zeros（N,1）;
Y=zeros（N,1）;
C=zeros（1,N）;
R=1:
N;
k=1;
whilek<
=N-1
%求列中最大的值赋给max
[max1,j]=max（abs（A（k:
N,k）））;
%交换行
C=A（k,:
）;
%C为A的第k列的值
A（k,:
）=A（j+k-1,:
%将A的第K列赋为最大
A（j+k-1,:
）=C;
d=R（k）;
R（k）=R（j+k-1）;
R（j+k-1）=d;
%主元为0的情况
ifA（k,k）==0
'
Aissingular.nouniquesolution'
break
end
%化为上三角
form=k+1:
N
mult=A（m,k）/A（k,k）;
A（m,k）=mult;
A（m,k+1:
N）=A（m,k+1:
N）-mult*A（k,k+1:
N）;
m=m+1;
k=k+1;
end
%对右端项做处理，但要保证行的交换相同,要注意R（k）的作用
Y
（1）=B（R
（1））;
fork=2:
Y（k）=B（R（k））-A（k,1:
k-1）*Y（1:
k-1）;
X（N）=Y（N）/A（N,N）;
fork=N-1:
-1:
1
X（k）=（Y（k）-A（k,k+1:
N）*X（k+1:
N））/A（k,k）;
例题：
用列主元Gauss消去法求解方程
的解。
输入：
>
A=[223;
476;
789];
B=[4;
7;
9];
lufact（A,B）
输出：
ans=
-2.3333
-0.3333
3.1111
4实验结论
上述matlab程序可以比较好的实现列主元Gauss消去法，且经过实际例子检验出这种方法是正确的。
列主元Gauss消去法是在每一步消元前,在主元所在的一列选取绝对值最大的元素作为主元素.而全主元Gauss消去法是在每一步消元前,在所有元素中选取绝对值最大的元素作为主元素.但由于运算量大增,实际应用中并不经常使用.
实验二Lagrange插值
熟悉Lagrange插值法的基本原理；
会用matlab程序实现Lagrange插值法；
会用Lagrange插值法求解实际问题。
给定（n+1）个插值节点
和对应的函数值
，利用n次拉格朗日插值多项式公式：
，其中
，可以得到插值区间任意x的函数值y为
从公式中可以看出，生成的多项式与用来插值的数据密切相关，数据变化则函数就要重新计算，所以当插值数据特别多的时候，计算量就会比较大。
由于Matlab中并没有现成的拉格朗日插值命令，下面先编写拉格朗日的函数文件：
functions=Lagrange（x,y,x0）
nx=length（x）;
ny=length（y）;
ifnx~=ny
warning（'
矢量x与y的长度应该相同'
）
return
m=length（x0）;
fori=1:
m
t=0.0;
forj=1:
nx
u=1.0;
fork=1:
nx;
ifk~=j
u=u*（x0（i）-x（k））/（x（j）-x（k））;
t=t+u*y（j）;
s（m）=t;
Return
编好M文件后，就可以用Lagrange插值函数进行插值计算了，下面根据x和y的数据求0.9出的值。
其中，
x=[0.70.50.40.8]
y=[-0.916291-0.693147-0.356675-0.223144];
x=[0.70.50.40.8];
y=[-0.916291-0.693147-0.356675-0.223144];
Lagrange（x,y,0.9）
1.3930
下面重新给定数据：
求解3.6,4.9,7.3三个点上的值。
x=[578910]
y=[107641];
x=[578910];
y=[107641];
Lagrange（x,y,[3.44.97.3]）
006.7729
从上图中可以看出，拉格朗日插值的一个特点：
拟合出的多项式通过每一个测量数据点。
由于Lagrange插值的n次插值基函数
lk（x）（k=0,1,2,…n）都依赖于全部插值结点，利用公式很容易得到插值多项式，但在增加或减少结点时，插值基函数
lk（x）（k=0,1,2,…n）也随之变化，必须全部重新计算。
实验三龙贝格求积公式
熟悉龙贝格求积公式的基本原理；
会用matlab程序实现龙贝格求积公式；
会用龙贝格求积公式求解定积分。
由于变步长梯形求积法算法简单，但其精度较差，收敛速度较慢，故利用变步长梯形法算法简单的优点，重新形成一个新算法，这就是龙贝格求积公式。
它是对近似值进行修正以后得到的更近似的公式，能自动改变积分步长，以使其相邻两次值的绝对误差或相对误差小于预先设定的允许误差。
龙贝格求积公式为：
根据上述实验原理，可以编写龙贝格求积公式的M文件romber.m：
function[R,quad,err,h]=romber（fun,a,b,n,delta）
M=1;
h=b-a;
err=1;
J=0;
R=zeros（4,4）;
R（1,1）=h*（fun（a）+fun（b））/2;
while（（err>
delta）&
&
（J<
n））||（J<
4）
J=J+1;
h=h/2;
s=0;
forp=1:
M
x=a+h*（2*p-1）;
c=fun（x）;
s=s+c;
R（J+1,1）=R（J,1）/2+h*s;
M=2*M;
forK=1:
J
R（J+1,K+1）=R（J+1,K）+（R（J+1,K）-R（J,K））/（4.^K-1）;
err=abs（R（J,J）-R（J+1,K+1））;
end
quad=R（J+1,J+1）
用龙贝格算法计算定积分
在matlab的命令栏输入：
fun=@（x）sqrt（4*x-x.^3）;
romber（fun,0,1,5,10.^（-10））
输出结果：
quad=
1.2580
0.866000000
1.11771.20150000
1.20831.23851.2410000
1.24081.25161.25251.252600
1.25241.25621.25651.25661.25660
1.25651.25791.25801.25801.25801.2580
龙贝格求积公式算法简单，高速有效，易于编制程序，适合于计算机上操作。
采用了提高阶数与减少步长这两种提高精度的措施。
能比较准确的求出定积分的值。
但它有一个明显的缺点是：
每当把区间对分后，就要对被积函数f（x）计算它在新分点处的值，而这些函数值的个数是成倍增加的。
实验四Runge-Kutta方法
熟悉Runge-Kutta方法的基本原理；
会用matlab程序实现Runge-Kutta方法；
会用Runge-Kutta方法求解微分方程。
三阶库塔方法（它的截断误差为
）:
用方程
中函数
在前一节点
上取值
的线性组合构造一个表示
的近似公式，从而避免了求
时用
的高阶导数。
该方法有多样推导方法，这里用数值积分法推导。
为此，先将微分方程
略加变形得出：
，对该式两边在相邻两节点
和
之间求积分，移项得出：
，采用不同的方法计算定积分
，便可得出数值积分的不同近似结果。
还可以构造更高阶的龙格—库塔公式，它的一般形式是：
，其中，
根据上述实验原理，可以编写Runge-Kutta方法的M文件：
functionR=rk3（f,a,b,ya,N）
h=（b-a）/N;
T=zeros（1,N+1）;
Y=zeros（1,N+1）;
T=a:
h:
b;
Y
（1）=ya;
forj=1:
k1=feval（f,T（j）,Y（j））;
k2=feval（f,T（j）+h/3,Y（j）+h*k1/3）;
k3=feval（f,T（j）+2*h/3,Y（j）+2*h*k2/3）;
Y（j+1）=Y（j）+h*（k1+3*k3）/4;
R=[T'
Y'
];
根据给出的条件求下列微分方程的解：
建立本题数据的M文件：
functionz=f（x,y）
z=y-x/y;
以f.m文件命名。
然后，回到Matlab的命令窗口，
rk3（'
f'
0,1,1,10）
回车后得到结果：
01.0000
0.10001.1003
0.20001.2025
0.30001.3081
0.40001.4187
0.50001.5359
0.60001.6613
0.70001.7965
0.80001.9433
0.90002.1035
1.00002.2791
龙格—库塔方法可以解决大部分微分方程的求解问题，并且可以较好较快的解出微分方程的解。
并且，龙格—库塔方法是由Euler法改进的，故其拥有Euler法的所有优点。
实验五Gauss-Seidel迭代法实例
熟悉Gauss-Seidel迭代法的基本原理；
会用matlab程序实现Gauss-Seidel迭代法；
会用Gauss-Seidel迭代法求解微分方程。
设线性方程组为AX=B，则Gauss-Seidel迭代法的迭代公式：
对应的矩阵表达式为：
称为Gauss-Seidel迭代矩阵，
上式称为Gauss-Seidel迭代法,简称为G-S迭代法
下面编写Gauss-Seidel迭代法的M文件gseid.m：
functionX=gseid（A,B,P,delta,max1）
N=length（B）;
fork=1:
max1
ifj==1
X
（1）=（B
（1）-A（1,2:
N）*P（2:
N））/A（1,1）;
elseifj==N
X（N）=（B（N）-A（N,1:
N-1）*（X（1:
N-1））'
）/A（N,N）;
else
X（j）=（B（j）-A（j,1:
j-1）*X（1:
j-1）-A（j,j+1:
N）*P（j+1:
N））/A（j,j）;
err=abs（norm（X'
-P））;
relerr=err/（norm（X）+eps）;
P=X'
;
if（err<
delta）|（relerr<
delta）
X
用Gauss-Seidel迭代法求解线性方程组Ax=B,其中，
A=[4-11;
4-81;
-215];
B=[7-2115]'
A=[4-11;
B=[7-2115]'
P=[000]'
delta=0.0000001;
max1=80;
X=gseid（A,B,P,delta,max1）
X=
2.00004.00003.0000
线性方程组的解法很多，如雅克比迭代法，高斯迭代法以及超松弛迭代法，较雅克比迭代法来看，高斯迭代法收敛较快.同样取精确到小数点后同一位数的近似解,高斯迭代法所需要的迭代次数次比雅克比迭代法需要的迭代次数少。
心得
短短几天的数值分析课程设计很快结束了，通过此次的训练我了解到一些简单的matlab编程方法，丰富了编程知识及常识，了解到了基础知识的重要性，以前认为的matlab都只是用来作图形的，但是课设完了我知道能利matlab解决比较复杂的数值计算问题。
通过我的努力以及遇到不会的问题主动向老师提问，我顺利的完成了这次语言课程设计，感到要想对一门课程了解深刻就必须亲自身入其中的去体验，这样才能发现服，也才能真正的学到知识，也明白了那句“纸上谈来终觉浅，深知此是要躬行自己的很多缺点，才能一一的去克”的真正意义。
纵然把结论书本背得滚瓜烂熟而不运用与实际，也不过是纸上谈兵罢了，终究在未来需要动手能力和实践人才的社会竞争中被淘汰，所以我们须要课程实习，需要工程训练，这样才能成为真正的二十一世纪合格的大学生！
参考文献
[1]石辛民.郝整清.基于Matlab的实用数值计算.清华大学出版社.北京交通大学出版社,2006
[2]苏金朋.阮沈勇.Matlab实用教程.电子工程出版社,2005
[3]张明辉.Matlab6.1..中国水利水电出版社,2008
[4]张池平.计算方法.哈尔滨工业大学出版社,2009
[5]杨万利.数值分析教程.国防工业出版社,2005
升级会员 