数值分析课程设计.docx

资源描述

数值分析课程设计.docx

《数值分析课程设计.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计.docx（105页珍藏版）》请在冰点文库上搜索。

数值分析课程设计.docx

数值分析课程设计

1.水手、猴子和椰子问题：

五个水手带了一只猴子来到南太平洋的一个荒岛上，发现那里有一大堆椰子。

由于旅途的颠簸，大家都很疲惫，很快就入睡了。

第一个水手醒来后，把椰子平分成五堆，将多余的一只给了猴子，他私藏了一堆后便又去睡了。

第二、第三、第四、第五个水手也陆续起来，和第一个水手一样，把椰子分成五堆，恰多一只猴子，私藏一堆，再去入睡，天亮以后，大家把余下的椰子重新等分成五堆，每人分一堆，正好余一只再给猴子，试问原先共有几只椰子？

试分析椰子数目的变化规律，利用逆向递推的方法求解这一问题。

解：

（1）算法分析：

求解这一问题可以用迭代递推算法。

首先分析椰子数目的变化规律，设最初的椰子数为

，即第一个水手所处理之前的椰子数，用

分别表示五个水手对椰子动了手脚以后剩余的椰子数目，则根据问题有

再用

表示最后每个水手平分得到的椰子数，于是有

所以

利用逆向递推的方法，有

有了逆向递推关系式，求解这一问题似乎很简单，但由于椰子数为一正整数，用任意的

作为初值递推出的

数据不一定是合适的。

这里用for循环语句结合break语句来寻找合适的

和

对任意的

递推计算出

，当计算结果为正整数时，结果正确,否则选取另外的

再次重新递推计算，直到计算出的结果

为正整数为止。

（2）MATLAB编程代码：

n=input（'输出n的值：

'）;

forx=1:

p=5*x+1;

fork=1:

p=5*p/4+1;

end

ifp==fix（p）

break

end

disp（[x,p]）

输出结果：

输出n的值：

1500

102315621

（3）结果分析：

在解决本题的过程中，运用了迭代的方法，每步还要判断X是否能被4整除，从而试探出结果。

2.当

时，选择稳定的算法计算积分

解：

用Romberg求解积分的程序：

新建m文件：

functionRomberg（a,b,n）

S=0;S1=0;S2=0;

fork=1:

n-1

S=S+f（a+k*（b-a）/n）;

end

forl=1:

2*n-1

S1=S1+f（a+1*（b-a）/（2*n））;

end

forj=1:

4*n-1

S2=S2+f（a+j*（b-a）/（4*n））;

end

formatlong

Tn=（b-a）/（2*n）*（f（a）+f（b）+2*S）;

T2n=（b-a）/（4*n）*（f（a）+f（b）+2*S1）;

T4n=（b-a）/（8*n）*（f（a）+f（b）+2*S2）;

Sn=（4*T2n-Tn）/（4-1）;

S2n=（4*T4n-T2n）/（4-1）;

T=（4^2*S2n-Sn）/（4^2-1）

新建m文件：

functionr=f（x）;

i=0：

100

r=exp（-i*x）/（exp（-x）+10）;

在CommandWindow窗口中输入：

>>Romberg（0,1,100）

输出结果为：

Columns1through4

.0954********

Columns5through8

-0.006725017522131-0.012846155883611-0.017013786280844-0.019979037056005

Columns9through12

.022*********

Columns13through16

-0.026957794976612-0.027621631857977-0.028163954943811-0.028609350346095

Columns17through20

.028*********

Columns21through24

.029*********

Columns25through28

-0.030279748540443-0.030322150088357-0.030348431688510-0.030360484153549

Columns29through32

-0.030359926829548-0.030348154524004-0.030326375025712-0.030295639345720

Columns33through36

-0.030256866275125-0.030210862467838-0.030158338971573-0.030099924918900

Columns37through40

-0.030036178931670-0.029967598672291-0.029894628883873-0.029817668190992

Columns41through44

.029*********

Columns45through48

.029*********

Columns49through52

.028*********

Columns53through56

.028*********

Columns57through60

.028*********

Columns61through64

.027*********

Columns65through68

.027*********

Columns69through72

-0.026789263359460-0.026674458567829-0.026559563207565-0.026444607685491

Columns73through76

-0.026329620537330-0.026214628557143-0.026099656916244-0.025984729272586

Columns77through80

.025*********

Columns81through84

.025*********

Columns85through88

.024*********

Columns89through92

.024*********

Columns93through96

.024*********

Columns97through100

.023*********

Column101

.023*********

3.绘制静态和动态的Koch分形曲线

问题描述：

从一条直线段开始，将线段中间的三分之一部分用一个等边三角形的另两条边代替，形成具有5个结点的新的图形；在新的图形中，又将图中每一直线段中间的三分之一部分都用一个等边三角形的另两条边代替，再次形成新的图形，这时，图形中共有17个结点。

这种迭代继续进行下去可以形成Koch分形曲线。

在迭代过程中，图形中的结点将越来越多，而曲线最终显示细节的多少取决于所进行的迭代次数和显示系统的分辨率。

Koch分形曲线的绘制与算法设计和计算机实现相关。

图1.1Koch曲线的形成过程

解：

在CommandWindow窗口中输入：

p=[00;100];n=2;

A=[cos（pi/3）-sin（pi/3）;sin（pi/3）cos（pi/3）];

fork=1:

d=diff（p）/3;m=4*n-3;

q=p（1:

n-1,:

）;p（5:

m,:

）=p（2:

n,:

）;

p（2:

m,:

）=q+d;

p（3:

m,:

）=q+d+d*A';

p（4:

m,:

）=q+2*d;

n=m;

end

plot（p（:

1）,p（:

2）,'k'）

axisequal

axisoff

（3）运行结果：

实验二

2.1小行星轨道问题：

一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在五个不同的对小行星作了五次观察，测得轨道上五个点的坐标数据（单位：

万公里）如下表所示：

X坐标

53605

58460

62859

66662

68894

Y坐标

6026

11179

16954

23492

68894

由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为：

现需要建立椭圆的方程以供研究。

（1）分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程组

以及方程组的系数矩阵和右端项b；

（2）用MARLAB求低阶方程的指令A＼b求出待定系数

;

（3）分别用直接法、Jacobi迭代法、Gauss-Seidel迭代法求出待定系数

解：

（1）

（2）A\b可得：

（3）直接法

新建m文件：

function[RA,RB,n,x]=liezy（A,b）

B=[Ab];n=length（b）;RA=rank（A）;

RB=rank（B）;zhica=RB-RA;

ifzhica>0,

disp（'因为RA~=RB，所以此方程无解'）

return

end

ifRA==RB

ifRA==n

disp（'RA=RB,有唯一解'）

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

b=B（1:

n,n+1）;A=B（1:

n,1:

n）;x（n）=b（n）/A（n,n）;

forq=n-1:

-1:

x（q）=（b（q）-sum（A（q,q+1:

n）*x（q+1:

n）））/A（q,q）;

end

else

disp（'因为RA=RB

end

在CommandWindow窗口中输入：

A=[28734960256460474603631267610721012052；341757160013070486801249004111692022358；3951253881213142297228743811612571833908；4443822244313204740855187406413332446984；474638323694927664724746383236137788137788];b=[-1;-1;-1;-1;-1];

2.2

（1）用Gauss列主元消去法、Gauss按比例列主元消去法、Cholesky分解求解下列线性方程组，并彼此互相验证。

（2）判断用Jacobi迭代法、Gauss-Seidel迭代法、SOR法（分别取

）解下列线性方程组的收敛性.若收敛，再用Jacobi迭代法、Gauss-Seidel迭代法、SOR法（分别取

）分别解线性方程组，并比较各种方法的收敛速度.

（3）用Cholesky分解求解下列线性方程组

（方程组的精确解是

）

解：

（1）Gauss列主元消去法

新建m文件：

function[RA,RB,n,x]=liezy（A,b）

B=[Ab];n=length（b）;RA=rank（A）;

RB=rank（B）;zhica=RB-RA;

ifzhica>0,

disp（'因为RA~=RB，所以此方程无解'）

return

end

ifRA==RB

ifRA==n

disp（'RA=RB,有唯一解'）

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

b=B（1:

n,n+1）;A=B（1:

n,1:

n）;x（n）=b（n）/A（n,n）;

forq=n-1:

-1:

x（q）=（b（q）-sum（A（q,q+1:

n）*x（q+1:

n）））/A（q,q）;

end

else

disp（'因为RA=RB

end

在CommandWindow窗口中输入：

>>A=[1,-1,2,1;-1,3,0,-3;2,0,9,-6;1,-3,-6,19];b=[1;3;5;7];[RA,RB,n,x]=liezy（A,b）

输出结果为：

RA=RB,有唯一解

RA=

RB=

-8.0000

0.3333

3.6667

2.0000

Gauss按比例列主元消去法

新建m文件：

function[RA,RB,n,X]=bilizy（A,b）

B=[Ab];n=length（b）;

RA=rank（A）;RB=rank（B）;

zhicha=RB-RA;

ifzhicha>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';

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:

m=B（k,p）/B（p,p）;

B（k,p:

n+1）=B（k,p:

n+1）-m*B（p,p:

n+1）;

end

b=B（1:

n,n+1）;A=B（1:

n,1:

n）;X（n）=b（n）/A（n,n）;

forq=n-1:

-1:

X（q）=（b（q）-sum（A（q,q+1:

n）*X（q+1:

n）））/A（q,q）;

end

else

disp（'请注意：

因为RA=RB

end

在CommandWindow窗口中输入：

>>A=[1,-1,2,1;-1,3,0,-3;2,0,9,-6;1,-3,-6,19];b=[1;3;5;7];[RA,RB,n,X]=bilizy（A,b）

输出结果为：

因为RA=RB=n，所以此方程组有唯一解

1-1211

-130-33

209-65

1-3-6197

1-1211

022-24

025-83

0-2-8186

1-1211

022-24

003-6-1

00-61610

RA=

RB=

-8.0000

0.3333

3.6667

2.0000

Cholesky分解

新建m文件：

functionCholesky（A,b）

n=size（A,1）

G=chol（A）;

（1）=b

（1）/G（1,1）;

fori=2:

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:

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）

在CommandWindow窗口中输入：

>>A=[1,-1,2,1;-1,3,0,-3;2,0,9,-6;1,-3,-6,19];b=[1;3;5;7];Cholesky（A,b）

输出结果为：

-8.00000.33333.66672.0000

（2）Jacobi迭代法

新建m文件：

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（'因为谱半径不小于1，所以Jacobi迭代序列发散，谱半径SRH和B的所有的特征值H如下：

'）

else

disp（'因为谱半径小于1，所以Jacobi迭代序列收敛，谱半径SRH和B的所有特征值H如下：

'）

end

SRH

在CommandWindow窗口中输入：

>>A=[1,-1,2,1;-1,3,0,-3;2,0,9,-6;1,-3,-6,19];H=JacobiC（A）

输出结果为：

因为谱半径小于1，所以Jacobi迭代序列收敛，谱半径SRH和B的所有特征值H如下：

SRH=

0.9669

0.5360

-0.9033

-0.5996

Gauss-Seidel迭代法

新建m文件：

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（'因为谱半径不小于1，所以Gauss-Seidel迭代序列发散，谱半径SRH和B的所有的特征值H如下：

'）

else

disp（'因为谱半径小于1所以Gauss-Seidel迭代序列收敛，谱半径SRH和B的所有特征值H如下：

'）

end

SRH

在CommandWindow窗口中输入：

>>A=[1,-1,2,1;-1,3,0,-3;2,0,9,-6;1,-3,-6,19];H=GSC（A）

输出结果为：

因为谱半径小于1所以Gauss-Seidel迭代序列收敛，谱半径SRH和B的所有特征值H如下：

SRH=

0.9349

0.2815

-0.0000

SOR法

新建m文件：

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（'因为谱半径不小于1，所以SOR迭代序列发散，谱半径SRH和B的所有特征值H如下：

'）

else

disp（'因为谱半径小于1，所以SOR迭代序列收敛，谱半径SRH和B的所有特征值H如下：

'）

end

SRH

（w=0.8）在CommandWindow窗口中输入：

>>A=[1,-1,2,1;-1,3,0,-3;2,0,9,-6;1,-3,-6,19];H=SORC（A,0.8）

输出结果为：

因为谱半径小于1，所以SOR迭代序列收敛，谱半径SRH和B的所有特征值H如下：

SRH=

0.9565

0.5033

0.0502

0.0662

（w=1.2）在CommandWindow窗口中输入：

>>A=[1,-1,2,1;-1,3,0,-3;2,0,9,-6;1,-3,-6,19];H=SORC（A,1.2）

输出结果为：

因为谱半径小于1，所以SOR迭代序列收敛，谱半径SRH和B的所有特征值H如下：

SRH=

0.9019

0.0392

0.0077+0.2125i

0.0077-0.2125i

（w=1.3）在CommandWindow窗口中输入：

>>A=[1,-1,2,1;-1,3,0,-3;2,0,9,-6;1,-3,-6,19];H=SORC（A,1.3）

输出结果为：

因为谱半径小于1，所以SOR迭代序列收敛，谱半径SRH和B的所有特征值H如下：

SRH=

0.8775

0.0946

-0.0538+0.3077i

-0.0538-0.3077i

（w=1.6）在CommandWindow窗口中输入：

>>A=[1,-1,2,1;-1,3,0,-3;2,0,9,-6;1,-3,-6,19];H=SORC（A,1.6）

输出结果为：

因为谱半径小于1，所以SOR迭代序列收敛，谱半径SRH和B的所有特征值H如下：

SRH=

0.6089

0.5901+0.0370i

0.5901-0.0370i

-0.2197+0.5679i

-0.2197-0.5679i

Jacobi迭代法

新建m文件：

functionJacobi（A,b,X0,P,error,max1）

[nn]=size（A）;

X=zeros（n,1）;

fork=1:

max1

forj=1:

X（j）=（b（j）-A（j,[1:

j-1,j+1:

n]）*X0（[1:

j-1,j+1:

n]））/A（j,j）;

end

errX=norm（X-X0,P）;

X0=X;X1=A\b;

if（errX

disp（'迭代次数k，精确解X1，近似解X分别是：

'）

return

end

if（errX>=error）

disp（'Jacobi迭代次数已超过最大次数'）

end

在CommandWindow窗口中输入：

>>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分别是：

173

X1=

-8.0000

0.3333

3.6667

2.0000

-7.9710

0.3397

3.6576

1.9965

Gauss-Seidel迭代法

新建m文件：

functionGS（A,b,X0,P,error,max1）

[nn]=size（A）

展开阅读全文