数值分析课程设计.docx
《数值分析课程设计.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计.docx（108页珍藏版）》请在冰点文库上搜索。
数值分析课程设计
数值分析实验指导书
考核标准：
及格：
独立完成12—15题，其中八组实验中每组至少做1题；
中：
独立完成16—23题，其中八组实验中每组至少做1题；
良：
独立完成24—31题，其中八组实验中每组至少做2题；
优：
独立完成32—40题，其中八组实验中每组至少做3题。
结束课程时，抽查上机考核
实验一
1.1水手、猴子和椰子问题：
五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。
由于旅途的颠簸，大家都很疲惫，很快就入睡了。
第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。
第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？
试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。
分析：
E0（整数）是最初的总数，Ei（i=1,2,3,4,5）是第i次平分后每堆的椰子数，E6是天亮后最后一次平分的每堆的椰子数，4Ei=5E（i+1）+1，满足的E0最小是1023
解答：
MATLAB程序：
>>n=input（'inputn:
'）;
forx=1:
n
p=5*x+1;
fork=1:
5
p=5*p/4+1;
end
ifp==fix（p）
disp（[x,p]）
end
end
inputn:
3000
102315621
204731246
1.2当
时，选择稳定的算法计算积分
.
解答：
由
.，有
所以有递推关系
显然求
时，
的误差的影响被输小了，是稳定算法。
不过要事先估计初值
由于：
取
时，由上式得
取
由此可做递归。
MATLAB程序：
I（100）=0.0009368;
fori=100:
-1:
1
I（i-1）=[（1-exp（-i））/i-I（i）]/10;
disp（[i-1,I（i-1）]
End
99.000.00090632
98.000.00091947
97.000.00092846
96.000.00093808
95.000.00094786
94.000.00095785
93.000.00096805
92.000.00097846
91.000.00098911
90.000.00099999
1.3绘制静态和动态的Koch分形曲线
问题描述：
从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5个结点的新的图形；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替，再次形成新的图形，这时，图形中共有17个结点。
这种迭代继续进行下去可以形成Koch分形曲线。
在迭代过程中，图形中的结点将越来越多，而曲线最终显示细节的多少取决于所进行的迭代次数和显示系统的分辨率。
Koch分形曲线的绘制与算法设计和计算机实现相关。
图1.1Koch曲线的形成过程
解答：
%绘制Koch曲线的Matlab程序。
clear;closeall;clf;
y=[1,0];t=4;
fork=1:
t
x=y;
n=length（x）-1;
form=0:
n-1;
z=（x（m+2）-x（m+1））/3;
y（4*m+1）=x（m+1）;
y（4*m+2）=x（m+1）+z;
y（4*m+3）=y（4*m+2）+z*（（1-sqrt（3）*i）/2）;
y（4*m+4）=x（m+1）+2*z;
end
y（4*n+1）=x（n+1）;
plot（y）;
axisequal;
holdon;figure;
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迭代法求出待定系数
.
分析：
将每一组的x和y坐标带入可得：
解答：
>>A=[
28734960256460474603631267610721012052;
3417571600130704868012497004111692022358;
3951253881213142297228743811612571833908;
4443822244313204740855187406413332446984;
474638323694927664724746383236137788137788
];b=[-1-1-1-1-1]’;
A\b
x=
1.0e+123*
0.0000
0.0001
0.0007
2.0703
7.7487
LU分解法：
>>[L,U]=lu（A）;x=U\（L\b）
x=
1.0e+123*
0.0000
0.0001
0.0007
2.0703
7.7487
Jacobi迭代法：
function[x,n]=jacobi（A,b,x0,eps,varargin）
%A:
系数矩阵，b:
常熟向量，x0:
初值，eps:
精度控制，varargin:
迭代步数控制
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin<3
disp（'error'）;
return;
elseifnargin==5
M=varargin{1};
end
D=diag（diag（A））;%求A的对角矩阵
L=-tril（A,-1）;%求A得下三角阵
U=-triu（A,1）;%求A的上三角阵
B=D\（L+U）;
f=D\b;
x=B*x0+f;
n=1;%迭代次数
whilenorm（x-x0）>=eps
x0=x;
x=B*x0+f;
n=n+1;
if（n>=M）
disp（'迭代次数太多，可能不收敛!
'）;
return;
end
end
A=[
53605^22*53605*60266026^22*536052*6026;
58460^22*58460*1117911179^22*584602*11179;
62859^22*62859*1695416954^22*628592*16954;
66662^22*66662*2349223492^22*666622*23492;
68894^22*68894*6889468894^22*688942*68894
];
b=[-1;-1;-1;-1;-1];x0=[0;0;0;0;0];
[x,n]=jacobi（A,b,x0）
迭代次数太多，可能不收敛!
x=
1.0e+123*
0.0000
0.0001
0.0007
2.0703
7.7487
Gauss-Seidel迭代法：
function[x,n]=gauseseidel（A,b,x0,eps,M）
%A:
系数矩阵，b:
常熟向量，x0:
初值，eps:
精度控制，M:
迭代步数控制
ifnargin==3
eps=1.0e-6;
M=200;
elseifnargin<3
disp（'error'）;
return;
elseifnargin==4
M=200;
end
D=diag（diag（A））;%求A的对角矩阵
L=-tril（A,-1）;%求A得下三角阵
U=-triu（A,1）;%求A的上三角阵
G=（D-L）\U;
f=（D-L）\b;
x=G*x0+f;
n=1;%迭代次数
whilenorm（x-x0）>=eps
x0=x;
x=G*x0+f;
n=n+1;
if（n>=M）
disp（'迭代次数太多，可能不收敛!
'）;
return;
end
end
A=[
53605^22*53605*60266026^22*536052*6026;
58460^22*58460*1117911179^22*584602*11179;
62859^22*62859*1695416954^22*628592*16954;
66662^22*66662*2349223492^22*666622*23492;
68894^22*68894*6889468894^22*688942*68894
];
b=[-1;-1;-1;-1;-1];x0=[0;0;0;0;0];
[x,n]=gauseseidel（A,b,x0）
迭代次数太多，可能不收敛!
x=
1.0e+005*
0.0000
0.0000
-0.0000
-0.2442
-2.7716
n=
200
2.2
（1）用Gauss列主元消去法、Gauss按比例列主元消去法、Cholesky分解求解下列线性方程组，并彼此互相验证。
（2）判断用Jacobi迭代法、Gauss-Seidel迭代法、SOR法（分别取
）解下列线性方程组的收敛性.若收敛，再用Jacobi迭代法、Gauss-Seidel迭代法、SOR法（分别取
）分别解线性方程组，并比较各种方法的收敛速度.
（3）用Cholesky分解求解下列线性方程组
（方程组的精确解是
）
解答：
Gauss列主元消去法：
function[RA,RB,n,x]=GaussLZY（A,b）
B=[Ab];n=length（b）;RA=rank（A）;
RB=rank（B）;zhicha=RB-RA;
ifzhicha>0,
disp（'此方程组无解'）
return
end
ifRA==RB
ifRA==n
disp（'此方程组有唯一解'）
x=zeros（n,1）;C=zeros（1,n+1）;
forp=1:
n-1
[Y,j]=max（abs（B（p:
n,p）））;C=B（p,:
）;
B（p,:
）=B（j+p-1,:
）;B（j+p-1,:
）=C;
fork=p+1:
n
m=B（k,p）/B（p,p）;
B（k,p:
n+1）=B（k,p:
n+1）-m*B（p,p:
n+1）;
end
end
b=B（1:
n,n+1）;A=B（1:
n,1:
n）;x（n）=b（n）/A（n,n）;
forq=n-1:
-1:
1
x（q）=（b（q）-sum（A（q,q+1:
n）*x（q+1:
n）））/A（q,q）;
end
else
disp（'此方程组有无穷多解'）
end
End
Gauss按比例列主元消去法：
function[RA,RB,n,X]=GaussBLZY（A,b）
B=[Ab];n=length（b）;RA=rank（A）;
RB=rank（B）;zhicha=RB-RA;
ifzhicha>0,
disp（'此方程组无解'）
return
end
ifRA==RB
ifRA==n
disp（'此方程组有唯一解'）
X=zeros（n,1）;C=zeros（1,n+1）;%#ok<*NASGU>
T=0;D=B';s=max（abs（D（1:
n,:
）））;
S=s';
forp=1:
n-1
[Y,j]=max（（abs（B（p:
n,p）））./S（p:
n））;
J=min（j）;
C=B（p,:
）;B（p,:
）=B（J+p-1,:
）;B（J+p-1,:
）=C;
T=S（p）;S（p）=S（J+p-1）;S（J+p-1）=T;
fork=p+1:
n
m=B（k,p）./B（p,p）;
B（k,p:
n+1）=B（k,p:
n+1）-m.*B（p,p:
n+1）;
end
end
b=B（1:
n,n+1）;A=B（1:
n,1:
n）;X（n）=b（n）./A（n,n）;
forq=n-1:
-1:
1
X（q）=（b（q）-sum（A（q,q+1:
n）*X（q+1:
n）））/A（q,q）;
end
else
disp（'此方程组有无穷多解'）
end
end
Cholesky分解法：
functionCholesky（A,b）
n=size（A,1）;
G=chol（A）;
f
（1）=b
（1）/G（1,1）;
fori=2:
n;
S1=0;
fork=1:
i-1
S1=S1+G（k,i）*f（k）;
end
f（i）=（b（i）-S1）/G（i,i）;
end
x（n）=f（n）/G（n,n）;
fori=n-1:
-1:
1
S2=0;
fork=1:
n-i
S2=S2+G（i,i+k）*x（i+k）;
end
x（i）=（f（i）-S2）/G（i,i）;
end
disp（x）
命令窗口~
Gauss列主元消去法
>>A=[1-121;-130-3;209-6;1-3-619];
>>b=[1;3;5;7];
>>[RA,RB,n,x]=GaussLZY（A,b）
此方程组有唯一解
RA=
4
RB=
4
n=
4
x=
-8.00000000000001
0.333333333333332
3.66666666666667
2
Gauss按比例列主元消去法
>>A=[1-121;-130-3;209-6;1-3-619];
>>b=[1;3;5;7];
>>[RA,RB,n,X]=GaussBLZY（A,b）
此方程组有唯一解
RA=
4
RB=
4
n=
4
X=
-8
0.333333333333333
3.66666666666667
2
Cholesky分解
>>A=[1-121;-130-3;209-6;1-3-619];
>>b=[1;3;5;7];
>>Cholesky（A,b）
Columns1through3
-8.000000000000020.3333333333333293.66666666666667
Column4
2
Jacobi迭代法：
functionH=JacobiC（A）
[nn]=size（A）;
D=diag（diag（A））;I=eye（n,n）;E=inv（D）;B=I-E*A;
H=eig（B）;SRH=norm（H,inf）;
ifSRH>=1
disp（'发散，谱半径SRH和B的特征值H如下：
'）
else
disp（'收敛，谱半径SRH和B的特征值H如下：
'）
end
SRH
functionJacobi（A,b,X0,P,error,max1）
[nn]=size（A）;
X=zeros（n,1）;
fork=1:
max1
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,P）;
X0=X;X1=A\b;
if（errXdisp（'迭代次数k，精确解X1，近似解X分别是：
'）
k
X1
X
return
end
end
if（errX>=error）
disp（'Jacobi迭代次数已超过最大次数'）
end
Gauss-Seidel迭代法：
functionH=GSC（A）
D=diag（diag（A））;U=-triu（A,1）;L=-tril（A,-1）;E=inv（D-L）;B=E*U;
H=eig（B）;SRH=norm（H,inf）;
ifSRH>=1
disp（'发散，谱半径SRH和B的特征值H如下：
'）
else
disp（'收敛，谱半径SRH和B的特征值H如下：
'）
end
SRH
functionGS（A,b,X0,P,error,max1）
[nn]=size（A）;
X=zeros（n,1）;
fork=1:
max1
forj=1:
n
XX=0;
fori=1:
n
ifiXX=XX+A（j,i）*X（i）;
end
ifi>j
XX=XX+A（j,i）*X0（i）;
end
end
X（j）=（b（j）-XX）/A（j,j）;
end
X
errX=norm（X-X0,P）;
X0=X;X1=A\b;
if（errXdisp（'迭代次数k，精确解X1，近似解X分别是：
'）
k
X1
X
return
end
end
if（errX>=error）
disp（'Gauss-Seidel迭代次数已超过最大次数'）
end
SOS法：
functionH=SORC（A,w）
D=diag（diag（A））;U=-triu（A,1）;L=-tril（A,-1）;E=inv（D-w*L）;
B=E*（w*U+（1-w）*D）;
H=eig（B）;SRH=norm（H,inf）;
ifSRH>=1
disp（'发散，谱半径SRH和B的特征值H如下：
'）
else
disp（'收敛，谱半径SRH和B的特征值H如下：
'）
end
SRH
functionSOR（A,b,X0,w,P,error,max1）
[nn]=size（A）;
X=zeros（n,1）;
fork=1:
max1
forj=1:
n
XX=0;
fori=1:
n
ifiXX=XX+A（j,i）*X（i）;
end
ifi>j
XX=XX+A（j,i）*X0（i）;
end
end
X（j）=（1-w）*X0（j）+w*（b（j）-XX）/A（j,j）;
end
X
errX=norm（X-X0,P）;
X0=X;X1=A\b;
if（errXdisp（'迭代次数k，精确解X1，近似解X分别是：
'）
k
X1
X
return
end
end
if（errX>=error）
disp（'SOR迭代次数已超过最大次数'）
end
Jacobi迭代法
>>A=[1-121;-130-3;209-6;1-3-619];
>>H=JacobiC（A）
收敛，谱半径SRH和B的特征值H如下：
SRH=
0.966881024532242
H=
0.966881024532242
0.536020655954138
-0.903317605697383
-0.599584074788997
>>A=[1-121;-130-3;209-6;1-3-619];
>>b=[1;3;5;7];
>>X0=[0;0;0;0];
>>Jacobi（A,b,X0,inf,0.001,1000）
迭代次数k，精确解X1，近似解X分别是：
k=
173
X1=
-8.00000000000002
0.333333333333329
3.66666666666667
2
X=
-7.97095795774273
0.339716804643296
3.65757261775243
1.99649138206938
Gauss-Seidel迭代法
A=[1-121;-130-3;209-6;1-3-619];
>>H=GSC（A）
收敛，谱半径SRH和B的特征值H如下：
SRH=
0.934888367800313
H=
0
0.934888367800313
0.281485901205534
-2.38957160758552e-017
>>A=[1-121;-130-3;209-6;1-3-619];
>>b=[1;3;5;7];
>>X0=[0;0;0;0];
>>GS（A,b,X0,inf,0.001,1000）
迭代次数k，精确解X1，近似解X分别是：
k=
99
X1=
-8.00000000000002
0.333333333333329
3.66666666666667
2
X=
-7.9864386402221
0.336173*********
3.66253302023792
1.99842937222506
SOR法（
）
>>SORC（A,0.8）
收敛，谱半径SRH和B的特征值H如下：
SRH=
0.956498743812875
ans=
0.956498743812875
0.503328665301619
0.050190208900126
.0662********
>>A=[1-121;-130-3;209-6;1-3-619];
>>b=[1;3;5;7];
>>X0=[0;0;0;0];
>>SOR（A,b,X0,0.8,inf,0.001,1000）
迭代次数k，精确解X1，近似解X分别是：
k=
139
X1=
-8.00000000000002
0.333333333333329
3.66666666666667
2
X=
-7.97865551731761
0.337878096178542
3.66010097132651
1.99749231060481
>>A=[1-121;-130-3;209-6;1-3-619];
>>SORC（A,1.2）
收敛，谱半径SRH和B的特征值H如下：
SRH=
0.901948033909749
ans=
0.901948033909749
0.039
升级会员 