数值分析课程设计.docx
《数值分析课程设计.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计.docx（120页珍藏版）》请在冰点文库上搜索。
数值分析课程设计
课程设计
数值分析课程设计
201年月日
设计题目
数值分析课程设计
成绩
课
程
设
计
主
要
内
容
实验一.12，实验二.1234，实验三.123，实验四.1234，实验五.1234，实验六.24，实验七.567，实验八.126，共计25题。
题目通过matlib程序解答。
指
导
教
师
评
语
建议：
从学生的工作态度、工作量、设计（论文）的创造性、学术性、实用性及书面表达能力等方面给出评价。
签名：
200年月日
1.1水手、猴子和椰子问题：
五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。
由于旅途的颠簸，大家都很疲惫，很快就入睡了。
第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。
第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？
试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。
解：
算法分析：
根据题意可以采用递归算法进行计算。
设设椰子起初的数目为x0,第一至第五次猴子在夜里藏椰子后，椰子的数目分别为x0，x1，x2，x3，x4，再设最后每个人分得x个椰子，设计以下算法：
n=input（'n='）;
forx=1:
n
t=5*x+1;
fork=1:
5
t=5*t/4+1;
end
ift==fix（t）
break
end
end
disp（[x,t]）
结果：
输入n=10000，得到没人分1024个椰子，椰子总数为15621。
改变x的范围，得到204731246；307146871；.......
等解。
1.2当
时，选择稳定的算法计算积分
.
解：
题目分析：
观察积分中的函数，发现下列特性：
由此得到递归式：
并且由上式可知求
时，
的误差的影响被缩小了。
当n=100时
的近似值为0。
则得到matli程序：
fprintf（'稳定算法:
\n'）
y0=0;
n=100;
plot（n,y0,'r*'）;
holdon
fprintf（'y[100]=%10.6f',y0）;
while
（1）
y1=1/10*[（1-exp（-n））/n-y0];
fprintf（'y[%10.0f]=%10.6f',n-1,y1）;plot（n-1,y1,'r*'）
if（n<=1）break;end
y0=y1;n=n-1;
ifmod（n,3）==0,fprintf（'\n'）,end,end
得到图形：
由此看出这是稳定的算法。
2.1小行星轨道问题：
一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在五个不同的对小行星作了五次观察，测得轨道上五个点的坐标数据（单位：
万公里）如下表所示：
P1
P2
P3
P4
P5
X坐标
53605
58460
62859
66662
68894
Y坐标
6026
11179
16954
23492
68894
由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为：
现需要建立椭圆的方程以供研究。
（1）分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程组
以及方程组的系数矩阵和右端项b；
（2）用MARLAB求低阶方程的指令A＼b求出待定系数
;
（3）分别用直接法、Jacobi迭代法、Gauss-Seidel迭代法求出待定系数
.
解：
（1）分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程组
以及方程组的系数矩阵和右端项b；
线性方程组为：
=
先算出系数矩阵。
使用MATLIB计算。
X0=[5360558460628596666268894];
Y0=[606211179169542349268894];
A=zeros（5）;
fori=1:
5
A（i,1）=X0（i）^2;A（i,2）=2*X0（i）*Y0（i）;A（i,3）=Y0（i）^2;A（i,4）=2*X0（i）;A（i,5）=2*Y0（i）;
end;
formatlongg;A
得到结果
A=
28734960256499070203674784410721012124
341757160013070486801249700411169202235
395125388121314229722874381112571833908
4443822244313204740855187406413332446984
474638323694927664724746383236137788137788
b=[-1-1-1-1-1]
（2）用MARLAB求低阶方程的指令A＼b求出待定系数;
B=[-1-1-1-1-1]';formatlongg;x=A\B
x=
-8.06820280371841e-011
-7.63620099622306e-011
-3.0801152978055e-010
-8.89025159419867e-006
2.02829368401655e-005
（3）分别用直接法、Jacobi迭代法、Gauss-Seidel迭代法求出待定系数.
直接法可使用用Gauss列主元消元法解线性方程组Ax=b
%用Gauss列主元消元法解线性方程组Ax=b
%RA,RB分别表示系数矩阵A和增广矩阵B的秩
%N表示向量b的维数,X是解向量
function[RA,RB,N,X’]=Gauss（A,b）
B=[Ab];
N=length（b）;
RA=rank（A）;
RB=rank（B）;
Diff=RB-RA;
ifDiff>0
disp（'因为RA~=RB，所以此方程组无解.'）
return
end
ifRA==RB
ifRA==N
disp（'因为RA=RB=N，所以此方程组有唯一解.'）
X=zeros（N,1）;
C=zeros（1,N+1）;
fori=1:
N-1
[Yk]=max（abs（B（i:
N,i）））;
C=B（i,:
）;
B（i,:
）=B（k+i-1,:
）;
B（k+i-1,:
）=C;
forj=i+1:
N
m=B（j,i）/B（i,i）;
B（j,i:
N+1）=B（j,i:
N+1）-m*B（i,i:
N+1）;
end
end
b=B（1:
N,N+1）;
A=B（1:
N,1:
N）;
X（N）=b（N）/A（N,N）;
fori=N-1:
-1:
1
X（i）=（b（i）-sum（A（i,i+1:
N）*X（i+1:
N）））/A（i,i）;
end
else
disp（'因为RA=RBend
end
X=X’;
End
因为RA=RB=N，所以此方程组有唯一解.
RA=
5
RB=
5
N=
5
X=
1.0e-004*
0.000006437215780
-0.000001957441606
0.000002799074569
-0.254776815358136
-0.001104956851916
Jacobi迭代法解线性方程组Ax=b
%用Jacobi迭代法解线性方程组Ax=b
%Norm：
范数的名称，Norm=1,2,inf;
%error：
近似解x的误差；
%Max：
迭代的最大次数;
functionX=Jacobi（A,b,X0,Norm,Error,Max）
[NN]=size（A）;
X=zeros（N,1）;
fori=1:
Max
forj=1:
N
X（j）=（b（j）-A（j,[1:
j-1,j+1:
N]）*X0（[1:
j-1,j+1:
N]））/A（j,j）;
end
X,
errX=norm（X-X0,Norm）;
X0=X;
iferrXX1=A\b;
disp（'迭代次数i,精确解X1和近似解X分别是:
'）
formatlong
i,X1,X,
return
end
end
iferrX>=Error
disp（'请注意：
Jacobi迭代次数已经超过最大迭代次数Max.'）
end
end
在MATLAB命令窗口输入：
A=[28734960256499070203674784410721012124;341757160013070486801249700411169202235;
3951253881213142297228743811612571833908;4443822244313204740855187406413332446984;
474638323694927664724746383236137788137788];
b=[-1-1-1-1-1]';
X0=[00000]';
Jacobi（A,b,X0,inf,0.001,100）;
结果
X=
1.0e-005*
-0.000034800813758
-0.000076508244513
-0.000347900972187
-0.750052503675257
-0.725752605451854
迭代次数i,精确解X1和近似解X分别是:
i=
1
X1=
1.0e-004*
0.000006437215780
-0.000001957441606
0.000002799074569
-0.254776815358137
-0.001104956851916
X=
1.0e-005*
-0.000034800813758
-0.000076508244513
-0.000347900972187
-0.750052503675257
-0.725752605451854
用Gauss-Seidel迭代法解线性方程组Ax=b
%用Gauss-Seidel迭代法解线性方程组Ax=b
%Norm:
范数的名称，Norm=1,2,inf;
%Error是近似解X的误差；Max是迭代的最大次数
functionX=G_S（A,b,X0,Norm,Error,Max）
[NN]=size（A）;
X=zeros（N,1）;
fori=1:
Max
forj=1:
N
X（j）=0;
ifj>1
X（j）=X（j）+A（j,1:
j-1）*X（1:
j-1）;
end
ifjX（j）=X（j）+A（j,j+1:
N）*X0（j+1:
N）;
end
X（j）=（b（j）-X（j））/A（j,j）;
end
X,
errX=norm（X-X0,Norm）;
X0=X;
iferrXX1=A\b;
disp（'迭代次数i,精确解X1和近似解X分别是:
'）
formatlong
i,X1,X,
return
end
end
iferrX>=Error
disp（'请注意：
Gauss-Seidel迭代次数已经超过最大迭代次数Max.'）
end
end
X=
1.0e-004*
-0.000003480081376
0.000001448627970
0.000002306743862
-0.002590235506932
-0.129368843047805
迭代次数i,精确解X1和近似解X分别是:
i=
1
X1=
1.0e-004*
0.000006437215780
-0.000001957441606
0.000002799074569
-0.254776815358137
-0.001104956851916
X=
1.0e-004*
-0.000003480081376
0.000001448627970
0.000002306743862
-0.002590235506932
-0.129368843047805
2.2
（1）用Gauss列主元消去法、Gauss按比例列主元消去法、Cholesky分解求解下列线性方程组，并彼此互相验证。
（2）判断用Jacobi迭代法、Gauss-Seidel迭代法、SOR法（分别取
）解下列线性方程组的收敛性.若收敛，再用Jacobi迭代法、Gauss-Seidel迭代法、SOR法（分别取
）分别解线性方程组，并比较各种方法的收敛速度.
（3）用Cholesky分解求解下列线性方程组
（方程组的精确解是
）
解：
（1）Gauss列主元消去法已在上题中给出。
1.用按比例主元消元法解线性方程组Ax=b
%用按比例主元消元法解线性方程组Ax=b
%RA,RB分别表示系数矩阵A和增广矩阵B的秩
%N表示向量b的维数,X是解向量
function[RA,RB,N,X]=Ratio（A,b）
B=[Ab];
N=length（b）;
RA=rank（A）;
RB=rank（B）;
RDiff=RB-RA;
ifRDiff>0
disp（'因为RA~=RB，所以此方程组无解.'）
return
end
ifRA==RB
ifRA==N
disp（'因为RA=RB=n，所以此方程组有唯一解.'）
X=zeros（N,1）;
C=zeros（1,N+1）;
T=0;
D=B';
s=max（abs（D（1:
N,:
）））;
S=s';
fori=1:
N-1
[Yk]=max（（abs（B（i:
N,i）））./S（i:
N））;
K=min（k）;
C=B（i,:
）;
B（i,:
）=B（K+i-1,:
）;
B（K+i-1,:
）=C;
T=S（i）;
S（i）=S（K+i-1）;
S（K+i-1）=T;
forj=i+1:
N
m=B（j,i）/B（i,i）;
B（j,i:
N+1）=B（k,i:
N+1）-m*B（i,i:
N+1）;
end
end
b=B（1:
N,N+1）;
A=B（1:
N,1:
N）;
X（N）=b（N）/A（N,N）;
fori=N-1:
-1:
1
X（i）=（b（i）-sum（A（i,i+1:
N）*X（i+1:
N）））/A（i,i）;
end
else
disp（'请注意：
因为RA=RBend
X=X’
End
2.用Cholesky分解（平方根法）解线性方程组
%用Cholesky分解解线性方程组
%A是方程组的系数矩阵，b是方程组的右边向量
functionCholesky（A,b）
[NN]=size（A）;
RA=rank（A）;
ifRA~=N
disp（'因为A的n阶行列式D等于零，所以A不能进行Cholesky分解.A的秩RA如下:
'）
RA,
return
end
ifA~=A'
disp（'因为A不是对称矩阵，所以A不能进行Cholesky分解.'）
return
end
ifRA==N
fori=1:
N
h（i）=det（A（1:
i,1:
i））;
end
D1=h（1:
N）;
fori=1:
N
ifh（1,i）<=0
disp（'因为A的i阶主子式不大于零，所以A不能进行Cholesky分解.A的各阶顺序主子式值Dl依次如下：
'）
D1,
return
end
end
ifh（1,i）>0
disp（'因为A是对称正定矩阵，所以A能进行Cholesky分解.A的下三角矩阵G和方程组的解X如下：
'）
forj=1:
N
U（1,j）=A（1,j）;
end
fori=2:
N
L（1,1）=1;
L（i,i）=1;
L（i,1）=A（i,1）/U（1,1）;
end
fork=2:
N
fori=2:
N
forj=2:
N
ifi>j
L（1,1）=1;
L（2,1）=A（2,1）/U（1,1）;
L（i,1）=A（i,1）/U（1,1）;
L（i,k）=（A（i,k）-L（i,1:
k-1）*U（1:
k-1,k））/U（k,k）;
else
U（k,j）=A（k,j）-L（k,1:
k-1）*U（1:
k-1,j）;
end
end
end
end
D1=eye（N）;
fori=1:
N
D1（i,i）=sqrt（U（i,i））;
end
G=L*D1;
Y
（1）=b
（1）/G（1,1）;
fori=2:
N
S1=0;
forj=1:
i-1
S1=S1+G（i,j）*Y（j）;
end
Y（i）=（b（i）-S1）/G（i,i）;
end
GT=G';
X（N）=Y（N）/GT（N,N）;
fori=N-1:
-1:
1
S2=0;
forj=1:
N-i
S2=S2+GT（i,i+j）*X（i+j）;
end
X（i）=（Y（i）-S2）/GT（i,i）;
end
G,X,
end
end
end
A=[1-121;-130-3;209-6;1-3-619];b=[1357]';
[RA,RB,n,x]=liezy（A,b）；
[RA,RB,n,x]=bilizy（A,b）；
cholesky（A,b）
结果：
列主元
因为RA=RB=n,所以此方程组有唯一解.
RA=
4
RB=
4
n=
4
x=
-8.000000000000005
0.333333333333332
3.666666666666668
2.000000000000000
比例主元
因为RA=RB=n,所以此方程组有唯一解.
RA=
4
RB=
4
n=
4
x=
-8.000000000000000
0.333333333333333
3.666666666666667
2.000000000000000
cholesky分解
S1=
0
S1=
0
S1=
0
x=
003.6666666666666732.000000000000003
x=
00.3333333333333293.6666666666666732.000000000000003
x=
-8.0000000000000200.3333333333333293.6666666666666732.000000000000003
-8.0000000000000200.3333333333333293.6666666666666732.000000000000003
对比可知，三种方法可相互验证。
（2）
先用谱半径判别Jacobi迭代法的收敛性
%用谱半径测试Jacobi迭代法的收敛性
%A是方程组的系数矩阵
functionJacobiTest（A）
[mn]=size（A）;
ifm~=n
disp（'系数矩阵必须为方阵'）
return
end
D=diag（diag（A））;
I=eye（n,n）;
B=I-inv（D）*A;
E=eig（B）;
SRH=norm（E,inf）;
ifSRH>=1
disp（'因为谱半径不小于1，所以Jacobi迭代序列发散,谱半径SRH和迭代矩阵的所有特征值如下：
'）
SRH,Eig=E’,
else
disp（'因为谱半径小于1，所以Jacobi迭代序列收敛,谱半径SRH和迭代矩阵的所有特征值如下：
'）
SRH,Eig=E’,
end
end
用谱半径测试Gauss-Seidel迭代法的收敛性
%A是方程组的系数矩阵
functionG_STest（A）
[mn]=size（A）;
ifm~=n
disp（'系数矩阵必须为方阵'）
return
end
D=diag（diag（A））;
U=-triu（A,1）;
L=-tril（A,-1）;
B=inv（D-L）*U;
E=eig（B）;
SRH=norm（E,inf）
ifSRH>=1
disp（'因为谱半径不小于1，所以Jacobi迭代序列发散,谱半径SRH和迭代矩阵的所有特征值如下：
'）
SRH,Eig=E’,
else
disp（'因为谱半径小于1，所以Jacobi迭代序列收敛,谱半径SRH和迭代矩阵的所有特征值如下：
'）
SRH,Eig=E’,
end
end
用谱半径测试SOR方法的收敛性
%A是方程组的系数矩阵，w松弛因子
functionSORTest（A,w）
[mn]=size（A）;
ifm~=n
d
升级会员 