数值分析课程设计大作业

本肥工学大孝

课程设计

设计题目《数值分析》课程设计

学生姓名____****

学号______####______________

专业班级_________

指导教师—

2013年07月20日

设计

《数值分析》课程设计成绩

题目

实验一：二题全做；实验一：四题全做；实验二：3.13.2已

课

程做；实验四4.14.24.34.4已做；实验五：5.15.25.35.4已

设做；实验六：6.26.36.4已做；实验七：7.37.57.67.8

计

已做；实验八：8.28.8已做。共计26道题。

主

要部分题通过matlab中notebook实现程序,notebook使用非

内常方便。由于篇幅有限部分题调用的函数源代码没有给出，

容

本作业所使用的m文件均在附录中

建议：从学生的工作态度、工作量、设计（论文）的创造性、学

术性、实用性及书面表达能力等方面给出评价。

指

导

教

师

评

语

签名：20年月日

1.1水手、猴子和椰子问题

算法分析:设椰子起初的数目为外，第一至第五次猴子在夜里藏椰子后，椰子的数目分别

为Po,Pl,p2,p3,P4,再设最后每个人分得X个椰子，由题意得：

41一

Z4+1左=°1，2,3,4.X=W(P5—1)所以p5=5x+l

利用逆向递推方法求解：

n=input(*n=*);

forx=l:n

p=5*x+l;

fork=l:5

p=5*p/4+1;

end

ifp==fix(p)

break

end

end

disp([x,p])

执行代码后得：n=102315621(输入n=1000000)

即最后每个人分得1023个椰子，椰子总数为15621

1.2当〃=。,1,2,‘100时，选择稳定的算法计算积分小〕；―也

由

rig-'+10

/.+10/=[----------dx=l,

'°0Joe-v+lO

,-(nl)x^

le++10e,11

:n

------------ax=e-^dx^-(l-e-)

°e~x+100n

得

，。4OF

/"4jWi=100,99，.'L

由上式可知求/“时，/"I的误差的影响被缩小了。n=100时4oo的近似值为。。

matlab代码为

fprintf('稳定算法:\n*)

yO=O;

n=100;

plot(n,yOJr*');

holdon

fprintf('y[100]=%10.6f\y0);

while(1)

yl=l/10*[(1-exp(-n))/n-yO];

1111

fprintf(y[%10.Of]=%10.6f,n-1zyl);plot(n-1,ylzr*)

if(n<=l)break;end

yO=yl;n=n-l;

ifmod(n,3)==0,fprintf(*\n1),end,end

(具体值已省略)

编程实现得下图。

由图可知，该算法是稳定的。

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

Koch曲线程序koch.m

functionkoch(a1,b1,a2,b2,n)

%koch(0,0,9,0,3)

%al,bl,a2,b2为初始线段两端点坐标,n为迭代次数

%例如a1=0;b1=0;a2=9;b2=0;n=3;

%第1-1次迭代时由各条线段产生的新四条线段的五点横、纵坐标存储在数组A、B中

[A,B]=sub_kochl(al,bl,a2,b2);

fori=l:n

forj=1:length(A)/5;

w=sub_koch2(A(l+5*(j-l):5*j),B(l+5*(j-l):5*j));

fork=l:4

[AA(5*4*(H)+5*(k-D+l：5*4*(j-D+5*(k-l)+5),BB(5*4*(j-l)+5*(k-D+l：5*4*(j-l)+5*(k-l)+5)]

=sub_koch1(w(k,1),w(k,2),w(k,3),w(k,4));

end

end

A=AA;

B=BB;

end

plot(A,B)

holdon

axisequal

%由以(ax,ay),(bx,by)为端点的线段生成新的中间三点坐标并把这五点横、纵坐标依次分别

存％储在数组A,B中

function[A,B]-sub_kochl(ax,ay,bx,by)

cx=ax+(bx-ax)/3;

cy=ay+(by-ay)/3;

ex=bx-(bx-ax)/3;

ey=by-(by-ay)/3;

L=sqrt((ex-cx).A2+(ey-cy).A2);

alpha=atan((ey-cy)./(ex-cx));

if(ex-cx)<0

alpha=alpha+pi;

end

dx=cx+cos(alpha+pi/3)*L;

dy二cy+sin(alpha+pi/3)*L;

A=[ax,cx,dx,ex,bx];

B=[ay,cy,dy,ey,by];

%把由函数sub_kochl生成的五点横、纵坐标A,B顺次划分为四组，分别对应四条折线段中

%每条线段两端点的坐标，并依次分别存储在4*4阶矩阵k中,k中第i(i=l,2,3,4)行数字代表

第％1条线段两端点的坐标

functionw=sub_koch2(A,B)

all=A(l);bll=B(l);

al2=A(2);bl2=B(2);

a21=A(2);b21=B(2);

a22=A(3);b22=B(3);

a31=A(3);b31=B(3);

a32=A(4);b32=B(4);

a41=A(4);b41=B(4);

a42=A(5);b42=B(5);

w=[all,bll,al2,b12;a21,b21,a22,b22;a31,b31,a32,b32;a41,b41,a42,b42];

%调用函数得到下图

n=5;

i=0;whilei<n

figure(i+l);koch(0,0,3,0,i)

i=i+l;pause(1)

end

0.51.52.5

将pause(l)去掉可得静态图

2.1小行星轨道问题

为了确定方程中的5个待定系数，需要将上述5个点的坐标代入上面的方程

2

々if+2a2xy+a3y+2a^x+2a5y=-1，得：

2

%遂：+2a2xiy1+a3yj+2a4Xj+2a.5y1=-1

2

+2a2x2y2+a3y2+2a4x2+2a5y2=-1

22

a1X3+2a2x3y3+a3y3+2a4x3+2a5y3=-1

22

a；x4+2a2x4y4+a3y4+2a4x4+2a5y4=-1

+2a?X5y5+a3y5?+2a4x5+2a5y5=-1

将这一包含5个未知数的线性方程组，写成矩阵的形式

,2

X12x*

2

X22x2y2

22xAX=b

X33y3

2

X42x4y4

2

X52X5丫5

(1)求解这一线性方程组，即可得到曲线方程的系数

X0=[5360558460628596666268894];

Y0=[606211179169542349268894];

A=zeros(5);X0(1);

fori=l:5

A(i,l)=X0(i)A2;A(i,2)=2*X0(i)*Y0(i);A(i,3)=Y0(i)A2;A(i,4)=2*X0(i

);A(i,5)=2*Y0(i);

end;

formatlongg;A

A=

28734960256499070203674784410721012124

341757160013070486801249700411169202235

3951253881213142297228743811612571833908

4443822244313204740855187406413332446984

474638323694927664724746383236137788137788

B=[-l-1-1-1-1]1/formatlongg;x=A\B

x

-8.06820280371841e-011

-7.63620099622306㊀-Oil

-3.0801152978055e-010

-8.89025159419867e-006

2.02829368401655㊀-005

(2)用Lu分解法解可得

formatlongg;

A=[28734960256499070203674784410721012124

3417571600130704868012497004111692022358;

3951253881213142297228743811612571833908;

4443822244313204740855187406413332446984;

474638323694927664724746383236137788137788];

B=[-l-1-1-1-1]1;

[L,U,flag]=LU_Decom(A),formatlongg;x=U\(L\B)

田L<5x5double>

12345

110000

21.18931000

31.37512.3175100

41.54653.98253.528810

51.651815.763970,0894224.19471

丑U<5x5double>

12345

12.8735e+096499070203674784410721012124

205.3409e+088.1264e+07-1.0589e+...7.9384e+03

3004.8576e+072.8381e+03-l.1608e+...

4000-317.9922716.0885

50000・8,6565e+...

flag=

OK

x=

-8.06820280370254e-011

-7.63620099622621e-011

-3.08011529780421e-010

-8.89025159420231e-006

2.02829368401615e-005

jacobi迭代法：

jacobic(A)

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

SRH=

4.41963931714337

ans

-4.41963931714337

1.5453292801696

0.757138763648732

1.11838895650533

0.9987823168197

GSC(A)

因为谱半径不小于1,所以Gauss-Seidel迭代序列发散，谱半径SRH和B的所有特征值H

如下：

SRH=

1.12218280703645

ans=

0

0.0914047045360394

1.10703104191813+0.1837839074507621

1.10703104191813-0.1837839074507621

0.99748223605848

2.2

⑴用Gauss列主元消去法、Gauss按比例列主元消去法、Cholesky分解求解下

列线性方程组，并彼此互相验证。

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

。=0.8,1.2,1.3,1.6）解下列线性方程组的收敛性.若收敛，再用Jacobi迭代法、

Gauss-Seidel迭代法、SOR法（分别取0=0.8,1.2,1.3,1.6）分别解线性方程组，并

比较各种方法的收敛速度.

再一冗2+2%3+%4=L

一再+3%2-3九4=3,

<

2石+9X3-6X4=5,

x1-3X2-6X3+19X4=7.

⑶用Cholesky分解求解下列线性方程组

42-402400王0

22-1-21320九2-6

-4-1141-8-35620

0-216-1-4-3323

21-8-1224-10-39

43-3-44111-44-22

025-3-101142x1-15

0063-3-4219445

(1)

A=[l-121;-130-3;209-6;1-3-619];b=[l357]1;

[RA,RB,n^x]=liezy(A,b),[RA,RB,n,x]=bilizy(Azb),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.666666666666673

2.000000000000003

x

00.3333333333333293.666666666666673

2.000000000000003

x=

-8.0000000000000200.3333333333333293.666666666666673

2.000000000000003

-8.0000000000000200.3333333333333293.666666666666673

2.000000000000003

对比可知，三种方法可相互验证。

(2)

x0=ones(4,l);[]

A=[l-121;-130-3;209-6;1-3-619];jacobic(A),GSC(A)

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

SRH=

0.966881024532242

ans=

0.966881024532242

0.536020655954138

-0.903317605697383

-0.599584074788997

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

下：

SRH=

0.934888367800313

ans=

0

0.934888367800313

0.281485901205534

-2.38957160758552e-017

jacobi迭代法：

A=[l-121;-130-3;209-6;1-3-619];b=[l357]1;x0=ones(4,1);

[x,n]=jacobi(A,b,xO)

x=

-7.99997357113847

0.333339142428815

3.66665839089021

1.99999680708368

n=

376

gauseidelfe

[x,n]=gauseidel(A,b,xO)

x

-7.99998673303378

0.333336112109258

3.66666262275452

1.99999846346783

n=

197

SOR法：

wl=0.8;w2=l.2;w3=l.3;w4=l.6;

SORC(A,wl)ASORC(A,w2),SORC(A,w3),SORC(A,w4)

w=

0.8

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

SRH=

0.956498743812875

ans=

0.956498743812875

0.503328665301619

0.050190208900126

0.0662163001140356

w=

1.2

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

SRH=

0.901948033909749

ans=

0.901948033909749

0.0392206356524567

0.00773145469258105+0.2125322052681921

0.00773145469258105-0.2125322052681921

w=

1.3

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

SRH=

0.877539492310148

ans=

0.877539492310148

0.0946172161251248

-0.0537947284866427+0.3076699917256641

-0.0537947284866427-0.3076699917256641

w=

1.6

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

SRH

0.608878775828679

ans=

0.590094939083101+0.0369506212620504i

0.590094939083101-0.0369506212620504i

-0.21966219054509+0.56787488560383i

-0.21966219054509-0.56787488560383i

w=0.8时[x,n]=SOR(Azb,xO,0.8,10人(-5),300)

x=

-7.99980126177984

0.333375649565234

3.66660553355675

1.99997665093436

n=

238

w=l.2时

[x,n]=SOR(A,b,xO,1.2,10A(-5),300)

x=

-7.99992226324999

0.333349197733832

3.66664330811682

1.99999119661784

n=

111

w=l.3时

A

[x,n]=SOR(A,bfxO,l.3r1O(-5))

x=

-7.99993630689315

0.333346074182271

3.66664773570359

1.99999290958378

n=

89

w=l.6时

[x,n]=SOR(A,b,xO,1.6,10A(-5))

x=

-7.99999410032529

0.333332406655426

3.66666433369261

1.99999911550516

n=

26

比较上述方法可得，w=1.6时的SOR法最为快速。

(3)Cholesky分解求解下列线性方程组

A=[42-402400;22-1-21320;-4-1141-8-356;

0-216-1-4-33;21-8-1224-10-3;43-3-44111-4;

025-3-101142;0063-3-4219];b=[0;-6;20;23;9;-22;-15;45];

Choleshy(A,b)

121.1481-140.112729.7515-60.152810.9120-26.7963

5.4259-2.0185

2.3设

20

A=1

2

(1)将矩阵A进行LU分解A=LU,其中U是上三角矩阵，L是主对角线

上的元素都是1的下三角矩阵。

(2)利用上述分解分别求解方程组

AX=4，AX=b2,AX

并由此求出逆矩阵A-，

(3)用LU分解求下列线性方程组的解

-42-3-121000O-王-5-

86-5-3650100x212

42-2-132-1031x33

0-215-13-1194%2

-426-167-3323X53

86-8571726-35%46

02-13-425301x713

1610-11-917342-122438

462-713920124X919

00-18-3-24-863-1_X10_-21

(方程组的精确解是x*=(1,-1,0,1,2,0,3,1,-1,2)T.)

(1)

A=[2023;181;2-315];[L,U,flag]=LU_Decom(A)

bl=[l00],;b2=[010],;b3=[001],;xl=A\bl,x2=A\b2,x3=A\b3

L

1.000000

0.05001.00000

0.1000-0.40511.0000

u=

20.00002.00003.0000

07.90000.8500

0015.0443

flag=

OK

X=[xlzx2,x3]

X=

0.0517-0.0164-0.0093

-0.00550.1237-0.0072

-0.00800.02690.0665

X即为A的逆。

(2)matlab中输入命令：[L,U,flag]=LU_Decom(A),formatlongg;x=U\(L\b)得

L为

田L<10x10double>

12345678910

11000000000

22100000000

31010000000

40-121000000

5-1210100000

621-32210000

701-21-1-0.40001000

8410-125.80003.0000100

9123-110.400019.0000-6.000010

1000-120-6.2000-12.00001.5000-0.52401

u为

田U<10x10double>

12345678910

142-3-1210000

2021-1230100

3001011-1031

40004-141232

5000031-21-12

60000051-120

70000000.40000.60002.80003

800000002.0000-13.0000-9.0000

900000000-125.0000-110.0000

10000000000-12.1400

flag=

OK

x=

1

-0.999999999999999

4.21884749357559e-015

0.999999999999999

1.99999999999999

1.68753899743024e-015

3

0.999999999999999

-1

2

2.4(i)用追赶法求解方程组xx=b

41

14

(a)A=

1

5-21

1

-25-21

0

1-25-21

0

(b)A=,b

1-25-21

0

1-25-2

0

1-25_20x20

(2)设计实验验证Hilbert矩阵的病态性，其中

c1]

1

n

11

H,,2n+1

iii

\nn+12〃一L

(1)

1

al=4*ones(30f1);a2=ones(29f1)*;a3=a2;b=ones(30zl)*;b(1f1)=2;b(30,1)=

2;A=diag(al)4-diag(a3,1)+diag(a2,-1);

x=chasing(A,b)

zhuiganfa(A,b)

ans=

0.5359

-0.1436

0.0385

-0.0103

0.0028

-0.0007

0.0002

-0.0001

0.0000

-0.0000

0.0000

-0.0000

0.0000

-0.0000

0.0000

0.0000

-0.0000

0.0000

-0.0000

0.0000

-0.0000

0.0000

-0.0001

0.0002

-0.0007

0.0028

-0.0103

0.0385

-0.1436

0.5359

(b)

1

al=5*ones(20f1)*;a2=-2*ones(19,1);a3=a2;a4=ones(18,1)';a5=a4;b=zeros

(20,1);b(l,l)=l;

A=diag(al)+diag(a3,1)+diag(a2,-1)+diag(a4,-2)+diag(a5,2);zhuigan

fa(A,b)

ans=

0.249999999999829

0.124999999999574

0.0624999999991047

0.0312499999981881

0.0156249999963656

0.00781249999272582

0.00390624998544897

0.00195312497089661

0.000976562441792561

0.000488281133584789

0.000244140392169412

0.00012206984683874

6.10342249274393—005

3.05157154798577㊀-005

1.5255063772205e-005

7.62194395065481e-006

3.79979610443202㊀-006

1.87754631042523e-006

8.94069671631063㊀-007

3.57627868652425—007

(2)

matlab代码为：

m=input(1inputm:=*);%输入矩阵的阶数N=[m];

N=[m];

fork=l:length(N)n=N(k);%矩阵的阶

H=hilb(n);%产生n阶Hilbert矩阵

disp(H)%输出n阶Hilbert矩阵

Hi=invhilb(n);%产生完全准确的n阶逆Hilbert矩阵

b=ones(n,l);%生成n阶全1向量

x_approx=H\b;%利用左除H求近似解

x__exact=Hi*b;%禾1」用准确逆Hilbert矩阵求准确解

ndb=norm(H*x_approx-b);nb=norm(b);ndx=norm(x_approx-

x__exact);

nx=norm(x__approx);

er_actual(k)=ndx/nx;%实际相对误差

K=cond(H);

er_approx(k)=K*eps;%最大可能的近似相对误差

er_max(k)=K*ndb/nb;%最大可能的相对误差

end

disp「Hilbert矩阵阶数，)，disp(N),formatshorte

disp(,实际误差er_actual*)zdisp(er__actual)zdisp(*')

disp(1近似的最大可能差er_approx*),disp(er_approx),disp(**)

disp(,最大可能误差er_max*),disp(er_max),disp(**)

执行代码，在matlab中输入4得

inputm:=

1.0000e+0005.0000e-0013.3333e-0012.5000e-001

5.0000e-0013.3333e-0012.5000e-0012.0000e-001

3.3333e-0012.5000e-0012.0000e-0011.6667e-001

2.5000e-0012.0000e-0011.6667e-0011.4286e-001

Hilb㊀rt矩阵阶数

4

实际误差er_actual

1.0284e-013

近似的最大可能差er_approx

3.4447e-012

最大可能误差er_max

4.7732e-011

多次执行代码，分别将m值设为81016100后得

\阶

481016100

实际误

1.0284e-0131.7310e-0071.9489e-0042.6458e+0029.1294e+123

差

近似的

最大可3.4447e-0123.3879e-0063.5583e-0031.3948e+0023.5120e+004

能误差

最大可

4.7732e-0113.8709e-0021.2703e+0036.2587e+0091.2236e+013

能误差

对于高阶系数矩阵来说，如果系数矩阵是Hilbert矩阵，则迭代结果误差较大，而且阶数越大,

误差就越大。

3.1

用Taylor级数的第，项多项式p(x)"1,2,,30),分别在区间[0,1/2]和

[0,l]上逼近正弦函数/(x)=sinx,%e[0,TT],并用计算机绘出上面31个函数的

图形。

symsx;

y=sin(x);

ezplot(y,[0,pi/2]);

gridon;

axis([0,pi/2,0,3]);

holdon;

form=l:2:30

p=taylor(y,x,m);

ezplot(p,[0,pi/2]);

axis([O,pi/2,O,1.5]);

end

fl=subs(sym(f),findsym(sym(f)),a);

f2=subs(sym(f),findsym(sym(f)),b);

if(fl==O)

root=a;

end

if(f2==0)

root=b;

end

if(fl*f2>0)

disp(*A^SIEpa0-Ey0p3E»yz60UO!*);

return;

else

tol=l;

fa=subs(sym(f),findsym(sym(f)),a);

fb=subs(sym(f),findsym(sym(f)),b);

fx=subs(sym(f),findsym(sym(f)),x);

dl=(fb-fa)/(b-a);

d2=(fx-fb)/(x-b);

d3=(d2-dl)/(x-a);

B=d2+d3*(x-b);

root=x-2*fx/(B+sign(B)*sqrt(BA2-4*fx*d3));

t=zeros(3);

t(1)=a；

t(2)=b;

t(3)=x;

while(tol>eps)

t(l)=t(2);

t(2)=t(3);

t(3)=root;

fl=subs(sym(f),findsym(sym(f)),t(1));

f2=subs(sym(f),findsym(sym(f))zt(2));

f3=subs(sym(f),findsym(sym(f)),t(3));

dl=(f2-fl)/(t(2)-t(l));

d2=(f3-f2)/(t(3)-t(2));

d3=(d2-dl)/(t(3)-t(l));

B=d2+d3*(t(3)-t(2));

root=t(3)-2*f3/(B+sign(B)*sqrt(BA2-4*f3*d3));

tol=abs(root-t(3));

end

end

f=inline('

xA10-55*xA9+1320*xA8-18150*xA7+157773*xA6-902055*xA5+3416930*xA4-8409

500*xA3+12753576*xA2-10628640*xAl+6328800');

Parabola(,

xA10-55*xA9+1320*xA8-18150*xA7+157773*xA6-902055*xA5+3416930*xA4-

8409500*xA3+12753576*xA2-10628640*x+63288001,10.6051,10.6051-

1.0128iz1.0e-8)

ans

2.41497968303632+2.7898128565398i

(3)

roots([1-5613209020553416930-840950012753576

-106286406328800])

ans=

21.7335121408968

7.35013402487903+7.79728509446717i

7.35013402487903-7.79728509446717i

4.35886185227404+3.3285467123814i

4.35886185227404-3.3285467123814i

5.1830729717618

2.4378000685152+2.79739872415432i

2.4378000685152-2.79739872415432i

0.394911498002447+1.01269497774644i

0.394911498002447-1.01269497774644i

原多项式方程的根变化之后多项式方程的根

10.6051+1.0127121.7335

10.6051-1.0127i7.3501+7.7973i

8.5850+2.7898i7.3501-7.79731

8.5850-2.789814.3589+3.3285i

5.5000+3.50581