数值分析上机作业2.docx

数值分析上机作业2.docx

《数值分析上机作业2.docx》由会员分享，可在线阅读，更多相关《数值分析上机作业2.docx（28页珍藏版）》请在冰豆网上搜索。

数值分析上机作业2

实验一

一、实验名称：

矩阵的LU分解.

二、实验目的：

用不选主元的LU分解和列主元LU分解求解线性方程组Ax=b,并比较这两种方法.

三、实验内容与要求

（1）用所熟悉的计算机语言将不选主元和列主元LU分解编成通用的子程序，然后用编写的程序求解下面的84阶方程组

将计算结果与方程组的精确解进行比较，并就此谈谈你对Gauss消去法的看法。

（2）写出追赶法求解三对角方程组的过程，并编写程序求该实验中的方程组

实验步骤：

（1）程序：

不选主元：

functionx=GAuss（A,b）

n=length（b）;

x=b;

fori=1:

n-1

forj=i+1:

n

mi=A（j,i）/A（i,i）;

b（j）=b（j）-mi*b（i）;

fork=i:

n

A（j,k）=A（j,k）-mi*A（i,k）;

end

AB=[A,b];

end

end

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

fori=n-1:

-1:

1

s=0;

forj=i+1:

n

s=s+A（i,j）*x（j）;

end

x（i）=（b（i）-s）/A（i,i）;

end

选主元：

functionx=GAussZhu（A,b）

n=length（b）;

x=b;

fork=1:

n-1

a_max=0;

fori=k:

n

ifabs（A（i,k））>a_max

a_max=abs（A（i,k））;r=i;

end

end

ifr>k

forj=k:

n

z=A（k,j）;A（k,j）=A（r,j）;A（r,j）=z;

end

z=b（k）;b（k）=b（r）;b（r）=z;

end

fori=k+1:

n

m=A（i,k）/A（k,k）;

forj=k:

n

A（i,j）=A（i,j）-m*A（k,j）;

end

b（i）=b（i）-m*b（k）;

end

end

ifabs（A（n,n））==0

return;

end

AbZhu=[A,b];

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

fori=n-1:

-1:

1

forj=i+1:

n

b（i）=b（i）-A（i,j）*x（j）;

end

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

end

生成n阶以a,b,c为元素的三对角阵：

functionA=Sanduijiao（a,b,c,n）

A=diag（b*ones（1,n）,0）+diag（c*ones（1,n-1）,1）+diag（a*ones（1,n-1）,-1）;

运行：

A2=Sanduijiaozhen（8,6,1,84）;

b2

（1）=7;

fori=2:

83

b2（i）=15;

end

b2（84）=14;

zx21=GAusszhu（A2,b2'）

x21=GAuss（A2,b2'）

运行结果：

zx21=

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

0.999999999999999

1

0.999999999999995

1.00000000000001

0.999999999999979

1.00000000000004

0.999999999999917

1.00000000000017

0.999999999999666

1.00000000000067

0.999999999998666

1.00000000000267

0.999999999994664

1.00000000001067

0.999999999978657

1.00000000004269

0.999999999914629

1.00000000017074

0.999999999658526

1.00000000068293

0.999999998634227

1.00000000273121

0.999999994538909

1.00000001091685

0.999999978187649

1.00000004353933

0.999999913262824

1.00000017210841

0.999999661246936

1.00000065565109

0.999998776117961

1.0000020980835

0.999997202555339

x21=

1

1

1

1

0.999999999999998

1

0.999999999999993

1.00000000000001

0.999999999999972

1.00000000000006

0.999999999999886

1.00000000000023

0.999999999999545

1.00000000000091

0.999999999998181

1.00000000000364

0.999999999992724

1.00000000001455

0.999999999970898

1.0000000000582

0.999999999883592

1.00000000023282

0.999999999534367

1.00000000093127

0.999999998137469

1.00000000372506

0.999999992549874

1.00000001490025

0.999999970199497

1.00000005960101

0.999999880797986

1.00000023840403

0.999999523191946

1.00000095361611

0.999998092767783

1.00000381446444

0.99999237107113

1.00001525785774

0.99996948428452

1.00006103143096

0.999877937138081

1.00024412572384

0.999511748552323

1.00097650289535

0.998046994209291

1.00390601158141

0.992187976837187

1.01562404632557

0.968751907349088

1.06249618530092

0.875007629401807

1.24998474118183

0.500030517694535

1.99993896437811

-0.999877927824961

4.99975585192486

-6.99951168894947

16.9990233182979

-30.9980463981918

64.9960918427676

-126.992179871071

256.984344484284

-510.968627937136

1024.93701174855

-2046.8730469942

4096.74218797682

-8190.46875190732

16383.8750076293

-32764.5000305175

65531.0001220701

-131055.000488281

262097.001953124

-524127.007812498

1048001.03125

-2094975.12499999

4185857.49999999

-8355328.99999997

16645128.9999999

-33028126.9999999

65007744.9999998

-125821439

234866688.999999

-402628606.999999

536838144.999998

从两个程序所得结果看，选主元的高斯消去法更加快更准确。

（2）追赶法

则方程组Ax=f称为三对角方程组。

设矩阵A非奇异，A有Crout分解A=LU，其中L为下三角矩阵，U为单位上三角矩阵，记

可先依次求出L，U中的元素后，令Ux=y，先求解下三角方程组Ly=f得出y，再求解上三角方程组Ux=y。

程序：

functions=trim（A,bb）

n=size（A,1）;

s=zeros（n,1）;

b=diag（A）;

a=diag（A,-1）;

c=diag（A,1）;

d=zeros（n,1）;

u=zeros（n-1,1）;

fori=1:

n-1

d

（1）=b

（1）;

u（i）=c（i）/d（i）;

d（i+1）=b（i+1）-a（i）*u（i）;

end

y=zeros（n,1）;

y

（1）=bb

（1）/d

（1）;

fori=2:

n

y（i）=（bb（i）-a（i-1）*y（i-1））/d（i）;

end

s（n）=y（n）;

fori=n-1:

-1:

1

s（i）=y（i）-u（i）*s（i+1）;

end

end

运行结果：

x21=

1

1

1

0.999999999999999

1

0.999999999999996

1.00000000000001

0.999999999999986

1.00000000000003

0.999999999999943

1.00000000000011

0.999999999999773

1.00000000000045

0.999999999999091

1.00000000000182

0.999999999996362

1.00000000000728

0.999999999985449

1.0000000000291

0.999999999941796

1.00000000011641

0.999999999767184

1.00000000046563

0.999999999068734

1.00000000186253

0.999999996274937

1.00000000745013

0.999999985099748

1.0000000298005

0.999999940398993

1.00000011920201

0.999999761595973

1.00000047680805

0.999999046383891

1.00000190723222

0.999996185535565

1.00000762892887

0.99998474214226

1.00003051571548

0.99993896856904

1.00012206286192

0.999755874276161

1.00048825144768

0.999023497104645

1.00195300579071

0.996093988418586

1.00781202316281

0.98437595367443

.0312********

0.937503814699085

1.12499237059819

0.750015258818165

1.49996948230547

6.10356218864183e-005

2.99987792782496

-2.99975585192486

8.99951168894947

-14.9990233182979

32.9980463981918

-62.9960918427676

128.992179871071

-254.984344484284

512.968627937136

-1022.93701174855

2048.8730469942

-4094.74218797682

8192.46875190732

-16381.8750076293

32766.5000305175

-65529.0001220701

131057.000488281

-262095.001953124

524129.007812498

-1047999.03125

2094977.12499999

-4185855.49999999

8355330.99999997

-16645126.9999999

33028128.9999999

-65007742.9999998

125821441

-234866686.999999

402628608.999999

-536838142.999998

实验二

实验步骤：

（1）程序：

1、平方根法：

functionx=Cholesky（A,b）

n=size（A）;

n=n

（1）;

x=A^-1*b;

fork=1:

n

A（k,k）=sqrt（A（k,k））;

A（k+1:

n,k）=A（k+1:

n,k）/A（k,k）;

forj=k+1:

n;

A（j:

n,j）=A（j:

n,j）-A（j:

n,k）*A（j,k）;

end

end

forj=1:

n-1

b（j）=b（j）/A（j,j）;

b（j+1:

n）=b（j+1:

n）-b（j）*A（j+1:

n,j）;

end

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

A=A';

forj=n:

-1:

2

b（j）=b（j）/A（j,j）;

b（1:

j-1）=b（1:

j-1）-b（j）*A（1:

j-1,j）;

end

b

（1）=b

（1）/A（1,1）;

disp（'平方根法的解即为'）;

2、改进的平方根法：

functionb=GJCholesky（A,b）

n=size（A）;

n=n

（1）;

v=zeros（n,1）;

forj=1:

n

fori=1:

j-1

v（i）=A（j,i）*A（i,i）;

end

A（j,j）=A（j,j）-A（j,1:

j-1）*v（1:

j-1）;

A（j+1:

n,j）=（A（j+1:

n,j）-A（j+1:

n,1:

j-1）*v（1:

j-1））/A（j,j）;

end

B=diag（A）;

D=zeros（n）;

fori=1:

n

D（i,i）=B（i）;

A（i,i）=1;

end

A=tril（A）;

forj=1:

n-1

b（j）=b（j）/A（j,j）;

b（j+1:

n）=b（j+1:

n）-b（j）*A（j+1:

n,j）;

end

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

A=D*（A'）;

forj=n:

-1:

2

b（j）=b（j）/A（j,j）;

b（1:

j-1）=b（1:

j-1）-b（j）*A（1:

j-1,j）;

end

b

（1）=b

（1）/A（1,1）;

disp（'改进平方根法解得的解即为'）;

>>A=Sanduijiaozhen（1,10,1,100）;

>>b

（1）=11;

>>fori=2:

99

b（i）=12;

end

>>b（100）=11;

>>x1=Cholesky（A,b'）

>>x2=GJCholesky（A,b'）

实验结果：

平方根法的解即为

x1=

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

1

0.999999999999999

1

1

1

1

1

1

1

1

1

1

改进平方根法解得的解即为：

都是1

（2）要编写b

程序：

functionb=Hil（）

fork=1:

40

form=1:

40

s=0;

t=s+1/（k+m-1）;

s=t;

end

b（k,1）=s;

end

运行：

A=hilb（40）;

b=Hil（）;

x1=Cholesky（A,b）

x2=GJCholesky（A,b）

运行结果：

平方根法：

x1=

0.033613490872085

0.000998020172119

.010*********

-0.001220703125000

-0.005859375000000

.023*********

.0937********

0

.0312********

0.109375000000000

-0.187********0000

-0.281250000000000

0.062500000000000

0.375000000000000

0.187********0000

-0.054687500000000

0.621093750000000

-0.968750000000000

-0.343750000000000

-0.312500000000000

-0.156********0000

1.062500000000000

-0.578125000000000

-0.156********0000

-1.468750000000000

-0.062500000000000

0.187********0000

0.375000000000000

-1.387695312500000

-0.312500000000000

0.375000000000000

0.375000000000000

0.125000000000000

-1.113281250000000

-0.375000000000000

1.406250000000000

.0937********

0.187********0000

-0.375000000000000

.0937********

改进平方根法解得的解即为

x2=

-0.000000000826711

-0.000000566568055

0.000049078892471

-0.001256760971985

0.014552136256574

-0.086919892706862

0.260792581223125

-0.216775282770511

-0.939281695721687

2.980973059440701

-2.928333325999537

0.134********5701

-0.275199533401944

3.726638261714109

.0855********

-2.521875229205066

1.349774981607172

1.372257861047557

-1.052597643628264

0.772507230724005

-0.813662309261239

0.627925510081358

0.140108105520086

-0.318044522760314

-2.068535661603255

2.985052481854060

0.159********9163

-1.511948502515179

-1.980619075304478

3.058323084908341

0.242689914699104

-2.194798484358696

2.432840067582593

-0.820439398876644

-0.807348359627573

-0.600672815071489

.029*********

-2.346694513217951

1.970094528584222

0.312046381978007

（3）

一、用高斯两种方法计算本题第一个方程组。

运行程序：

clc;clear

A=Sanduijiaozhen（1,10,1,100）;

b

（1）=11;

fori=2:

99

b（i）=12;

end

b（100）=11;

x1=GAuss（A,b'）

x2=GAussZhu（A,b'）

两个结果都是1。

二、用高斯两种方法计算本实验第二个方程组。

运行程序：

A=hilb（40）;

b=Hil（）;

x1=GAuss（A,b）

x2=GAussZhu（A,b）

运行结果都是0，最后一位是1。

高斯消去法所得的结果与平方根法和改进的平方根法求得的结果差距很大，而且高斯消去法所得的结果大部分为零，显然平方根法和改进的平方根法求得的结果与方程的精确解比高斯消去法的更接近，更准确。

但不管是哪一类算法都只能在预定的计算步骤内或给定的精度内得到近似解，有一定的误差。

实验三

实验步骤：

（1）生成一个

，随机生成

的一条主对角线元素和二条次对角线元素，使

为严格对角占优的三对角阵和

。

程序：

clear,clc

n=300;

A=zeros（n）;

A=diag（rand*ones（n,1））+diag（rand*ones（n-1,1）,-1）+diag（rand*ones（n-1,1）,1）;

fori=1:

n

A（i,i）=sum（abs（A（i,:

）））;

end

A=A;

b=rand（n,1）;

tic

x1=A\b;

t1=toc

tic

x2=inv（A）*b;

t2=toc

tic

x3=trim（A,b）;%实验一所编追赶法函数

t3=toc

实验结果：

t1=

6.182670546696735e-004

t2=

.023*********

t3=

1.035628416524151e-0