大连理工大学矩阵上机.docx

资源描述

大连理工大学矩阵上机.docx

《大连理工大学矩阵上机.docx》由会员分享，可在线阅读，更多相关《大连理工大学矩阵上机.docx（47页珍藏版）》请在冰豆网上搜索。

大连理工大学矩阵上机.docx

大连理工大学矩阵上机

大连理工大学

DALIANUNIVERSITYOFTECHNOLOGY

研究生作业

课程名称：

矩阵与数值分析

研究生姓名：

学号：

作业成绩：

任课教师签名

能源与动力学院

上交作业时间：

2015年12月25日

1.设

，其精确值为

（1）编制按从大到小顺序

，计算

的通用程序；

（2）编制按从小到大顺序

，计算

的通用程序；

（3）按两种顺序分别计算

，并指出有效位数（编制时的单精度）；

（4）通过本次上机题，你明白了什么？

MATLAB程序

function[s1,s2]=s（N）

formatlongg

s1=0;s2=0;

forj=2:

s1=1.0/（j*j-1）+s1;

end

forj=N:

-1:

s2=1.0/（j*j-1）+s2;

end

计算结果

>>[s1,s2]=s（1.0e2）

s1=

0.740049504950495

s2=

0.740049504950495

>>[s1,s2]=s（1.0e4）

s1=

0.749900004999506

s2=

0.7499000049995

>>[s1,s2]=s（1.0e6）

s1=

0.749999000000522

s2=

0.7499990000005

结果分析

计算机做加减法时，运算次序会影响结果，计算和时应先安排绝对值小的数参与运算，这样能取得较高的精度。

2.秦九韶算法。

已知n次多项式

，用秦九韶算法编写通用的程序计算函数在

点的值，并计算

在23点的值。

（提示：

编写程序时，输入系数向量和点

，输出结果，多项式的次数可以通过向量的长度来判断）.

MATLAB程序

A=input（'ÇëÊäÈëÏµÊý,ÓÉ¸ß´ÎÃÝ¿ªÊ¼'）;

n=input（'ÇëÊäÈë¼ÆËã±äÁ¿µÄÖµ'）;

len=length（A）;

val=zeros（len）;

val

（1）=A

（1）;

%%printf（'len=%c',len）

fori=2:

len

%disp（val（i-1））

%disp（n）

val（i）=val（i-1）*n+val（i）;

end

%printf（'?

¨¢?

%f',val（len））

val（len）

计算结果

请输入系数，由高次幂开始[73-511]

请输入计算变量的值23

ans=

85169

3.分别用Gauss消元法和列主元消去法编程求解方程组Ax=b，其中

。

Gauss消元法

functionx=gauss（a,b）

n=length（b）;

fork=1:

n-1

ifa（k,k）==0

fprintf（'Error:

thepivotelementequaltozero!

\n',k）;

return;

end

index=[k+1:

n];

m=-a（index,k）/a（k,k）;

a（index,index）=a（index,index）+m*a（k,index）;

b（index）=b（index）+m*b（k）;

end

x=zeros（n,1）;

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

fori=n-1:

-1:

x（i）=（b（i）-a（i,[i+1:

n]）*x（[i+1:

n]））/a（i,i）;

end

计算结果

>>a=[31,-13,0,0,0,-10,0,0,0;-13,35,-9,0,-11,0,0,0,0;0,-9,31,-10,0,0,0,0,0;0,0,-10,79,-30,0,0,0,-9;0,0,0,-30,57,-7,0,-5,0;0,0,0,0,-7,47,-30,0,0;0,0,0,0,0,-30,41,0,0;0,0,0,0,-5,0,0,27,-2;0,0,0,-9,0,0,0,-2,29]

>>b=[-15;27;-23;0;-20;12;-7;7;10]

x=gau（a,b）

-0.289233816015754

0.345435715779115

-0.712811731086879

-0.220608510570529

-0.430400432704022

0.154308739838311

.0578********

0.201053894823681

0.290228661879745

列主元消去法

function[x]=gauss（a,b）

n=length（a）;

x=zeros（n,1）;

a=[ab];

fork=1:

n-1

max=k;

fori=k+1:

ifa（i,k）>a（max,k）

max=i;

end

temp=a（k,k:

n+1）;

a（k,k:

n+1）=a（max,k:

n+1）;

a（max,k:

n+1）=temp;

fori=k+1:

a（i,k）=-a（i,k）/a（k,k）;

a（i,k+1:

n+1）=a（i,k+1:

n+1）+a（i,k）*a（k,k+1:

n+1）;

end

x（n,1）=a（n,n+1）/a（n,n）;

fori=n-1:

-1:

sum=0;

forj=i+1:

sum=sum+x（j,1）*a（i,j）;

end

x（i,1）=（a（i,n+1）-sum）/a（i,i）;

end

计算结果

>>b=[-15;27;-23;0;-20;12;-7;7;10]

>>[x]=gauss（a,b）

-0.289233816015754

0.345435715779115

-0.712811731086879

-0.220608510570529

-0.430400432704022

0.154308739838311

.0578********

0.201053894823681

0.290228661879745

4.编程求解题3中矩阵A的LU分解及列主元的LU分解（求出L,U和P），并用LU分解的方法求A的逆矩阵及A的行列式.

LU分解

function[L,U]=mylu（a）

[n,n]=size（a）;

fori=1:

L（i,i）=1;

end

U（1,1）=a（1,1）/L（1,1）;

ifL（1,1）*U（1,1）==0

fprintf（'Factorizationimpossible'）

else

forj=2:

U（1,j）=a（1,j）;

L（j,1）=a（j,1）/U（1,1）;

end

fori=2:

n-1

sum=0;

fork=1:

i-1

sum=sum+L（i,k）*U（k,i）;

end

U（i,i）=（a（i,i）-sum）/L（i,i）;

ifL（i,i）*U（i,i）==0

fprintf（'Factorizationimpossible'）

else

forj=i+1:

h=0;

s=0;

fork=1:

i-1

h=h+L（i,k）*U（k,j）;

s=s+L（j,k）*U（k,i）;

end

U（i,j）=1/L（i,i）*（a（i,j）-h）;

L（j,i）=1/U（i,i）*（a（j,i）-s）;

end

sum=0;

fork=1:

n-1

sum=sum+L（n,k）*U（k,n）;

U（n,n）=（a（n,n）-sum）/L（n,n）;

end

计算结果

>>[L,U]=mylu（a）

Columns1through3

100

-0.41935483870967710

0-0.3045851528384281

00-0.353872899362565

000

Columns4through6

000

100

-0.39755492585681310

0-0.1569436359494821

00-0.653976718762685

0-0.112102597106773-0.0175453757582425

-0.119266477757044-0.0833908821051642-0.0142268147798013

Columns7through9

000

100

-0.024618525643364810

-0.0199621375629645-0.09231647025909681

Columns1through3

31-130

029.5483870967742-9

0028.2587336244541

000

Columns4through6

00-10

0-11-4.19354838709677

-10-3.35043668122271-1.27729257641921

75.4612710063743-31.185********5-0.451999227351748

044.6019996774714-7.17969451931716

0045.8731926371318

000

Columns7through9

000

00-9

0-5-3.57799433271131

-30-0.784718179747412-0.561543439982356

21.3806984371195-0.513187420344639-0.367236336322372

026.4130849214705-2.41999576495225

0027.3895043327769

LU分解计算a的行列式和逆矩阵

>>det（U）

ans=

6.1817e+13

>>inv（U）*inv（L）

ans=

0.03900.01600.00580.00360.00730.01760.01290.00140.0012

0.01590.03780.01310.00650.01210.00970.00710.00240.0022

0.00490.01160.03810.00770.00690.00390.00280.00150.0025

0.00080.00200.00640.01810.01050.00330.00240.00240.0058

0.00050.00110.00360.01010.02430.00700.00510.00480.0035

0.00010.00030.00100.00280.00680.04190.03060.00130.0010

0.00010.00020.00070.00210.00500.03060.04680.00100.0007

0.00010.00020.00080.00230.00480.00140.00100.03820.0033

0.00030.00060.00200.00580.00360.00110.00080.00340.0365

列主元LU分解

function[l,u,p]=lielu（a）

[m,n]=size（a）;

ifm~=n

error（'¾ØÕó²»ÊÇ·½Õó'）

return

end

ifdet（a）==0

error（'¾ØÕó²»ÄÜ±»Èý½Ç·Ö½â'）

end

u=a;p=eye（m）;l=eye（m）;

fori=1:

forj=i:

t（j）=u（j,i）;

fork=1:

i-1

t（j）=t（j）-u（j,k）*u（k,i）;

end

a=i;b=abs（t（i））;

forj=i+1:

ifb

b=abs（t（j））;

a=j;

end

ifa~=i

forj=1:

c=u（i,j）;

u（i,j）=u（a,j）;

u（a,j）=c;

end

forj=1:

c=p（i,j）;

p（i,j）=p（a,j）;

p（a,j）=c;

end

c=t（a）;

t（a）=t（i）;

t（i）=c;

end

u（i,i）=t（i）;

forj=i+1:

u（j,i）=t（j）/t（i）;

end

forj=i+1:

fork=1:

i-1

u（i,j）=u（i,j）-u（i,k）*u（k,j）;

end

l=tril（u,-1）+eye（m）;

u=triu（u,0）;

end

计算结果

>>[l,u,p]=lielu（a）

1.000000000000

-0.41941.00000000000

0-0.30461.0000000000

00-0.35391.000000000

000-0.39761.00000000

0000-0.15691.0000000

00000-0.65401.000000

0000-0.1121-0.0175-0.02461.00000

000-0.1193-0.0834-0.0142-0.0200-0.09231.0000

31.0000-13.0000000-10.0000000

029.5484-9.00000-11.0000-4.1935000

0028.2587-10.0000-3.3504-1.2773000

00075.4613-31.1856-0.452000-9.0000

000044.6020-7.17970-5.0000-3.5780

0000045.8732-30.0000-0.7847-0.5615

00000021.3807-0.5132-0.3672

000000026.4131-2.4200

0000000027.3895

100000000

010000000

001000000

000100000

000010000

000001000

000000100

000000010

000000001

5.编制程序求解矩阵A的cholesky分解，并用程序求解方程组Ax=b，其中

，

MATLAB程序

A=[21-51

1-507

021-1

16-1-4

];

B=[13,-9,6,0];

%disp（A）

dim=size（A,1）;

L=zeros（dim,dim）;

y=[]

x=zeros（dim,1）;

fori=1:

dim

forj=1:

dim

ifi>=（j+1）

L（i,j）=（A（i,j）-sum（L（i,1:

（j-1））.*L（j,1:

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

else

L（j,j）=sqrt（A（j,j）-sum（L（j,1:

j-1）.*L（j,1:

j-1）））;

end

L_T=L.';

fori=1:

dim

ifi==1

（1）=B

（1）/L（1,1）;

else

y（i）=（B（i）-sum（L（i,1:

（i-1））.*y（1:

（i-1））））/L（i,i）;

end

fori=dim:

-1:

ifi==dim

x（i）=y（i）/L（i,i）;

else

x（i）=（y（i）-sum（L（（i+1）:

dim,i）.*x（（i+1）:

dim）））/L（i,i）;

end

%disp（'（（（（（（（（'）

%disp（x（i））

end

disp（'x'）

disp（x）

计算结果

52.2500

-38.7500

30.7500

-52.7500

6.用追赶法编制程序求解方程组Ax=b，其中

，

MATLAB程序

A=[4,2,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6];

a=[032-1];c=[213];b=[4-256];d=[62105];

n=length（b）;

u0=0;y0=0;a

（1）=0;

（1）=b

（1）-a

（1）*u0;

（1）=（d

（1）-y0*a

（1））/L

（1）;

（1）=c

（1）/L

（1）;

fori=2:

（n-1）

L（i）=b（i）-a（i）*u（i-1）;

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

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

end

L（n）=b（n）-a（n）*u（n-1）;

y（n）=（d（n）-y（n-1）*a（n））/L（n）;

x（n）=y（n）;

fori=（n-1）:

-1:

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

end

计算结果

ans=

7.已知

编程求解矩阵A的QR分解.

MATLAB程序

function[Q,R]=qrhouseholder（A）

dim=size（A,1）;

R=A;

Q=eye（dim）;

fori=1:

（dim-1）

x=R（i:

dim,i）;

y=[1;zeros（dim-i,1）];

Ht=householder（x,y）;

H=blkdiag（eye（i-1）,Ht）;

Q=Q*H;

R=H*R;

end

function[H,rho]=householder（x,y）

x=x（:

）;

y=y（:

）;

iflength（x）~=length（y）

error（'TheColumnVectorsXandYMustHaveTheSameLength!

'）;

end

rho=-sign（x

（1））*norm（x）/norm（y）;

y=rho*y;

v=（x-y）/norm（x-y）;

I=eye（length（x））;

H=I-2*v*v';

End

计算结果

1.00001.000000

-1.00003.0000-0.50000.5000

-2.00002.00001.50000.5000

-2.00002.0000-0.50002.5000

>>[Q,R]=qrhouseholder（C）

-0.3162-0.7071-0.25820.5774

0.3162-0.70710.2582-0.5774

0.6325-0.0000-0.7746-0.0000

0.6325-0.00000.51640.5774

-3.16233.16230.47432.0555

0.0000-2.82840.3536-0.3536

0.00000.0000-1.54921.0328

-0.0000-0.0000-0.00001.1547

8.分别应用Jacobi迭代法和Gauss-Seidel迭代法求解如下方程组

Jacobi法主程序

function[x,k,succ]=Jacobi（A,b,eps,it_max）

n=length（A）;

k=0;

x=zeros（n,1）;

y=zeros（n,1）;

succ=1;

while1

fori=1:

y（i）=b（i）;

forj=1:

ifj~=i

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

end

ifk==it_max

succ=0;return;

end

y（i）=y（i）/A（i,i）;

end

ifnorm（y-x,2）

break;

end

x=y;k=k+1;

end

计算结果

>>A=[4,-1,1;4,8,1;-2,1,5];

>>b=[7,-21,15];

>>[x,k,succ]=Jacobi（A,b,1e-6,1000）

0.0606

-3.1111

3.6465

succ=

G-

展开阅读全文