大连理工大学矩阵与数值分析上机作业Word文档下载推荐.docx

大连理工大学矩阵与数值分析上机作业Word文档下载推荐.docx

《大连理工大学矩阵与数值分析上机作业Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《大连理工大学矩阵与数值分析上机作业Word文档下载推荐.docx（45页珍藏版）》请在冰豆网上搜索。

Test2.m

clearall;

clc;

n=100;

%区间

h=2*10^（-15）/n;

%步长

x=-10^（-15）:

h:

10^（-15）;

%第一种原函数

f1=zeros（1,n+1）;

fork=1:

n+1

ifx（k）~=0

f1（k）=log（1+x（k））/x（k）;

else

f1（k）=1;

end

end

subplot（2,1,1）;

plot（x,f1,'

-r'

axis（[-10^（-15）,10^（-15）,-1,2]）;

legend（'

原图'

%第二种算法

f2=zeros（1,n+1）;

d=1+x（k）;

if（d~=1）

f2（k）=log（d）/（d-1）;

f2（k）=1;

subplot（2,1,2）;

plot（x,f2,'

第二种算法'

显然第二种算法结果不准确，是因为计算机中的舍入误差造成的，当

时，

，计算机进行舍入造成

恒等于1，结果函数值恒为1。

程序：

秦九韶算法:

QinJS.m

functiony=QinJS（a,x）

%y输出函数值

%a多项式系数，由高次到零次

%x给定点

n=length（a）;

s=a

（1）;

fori=2:

n

s=s*x+a（i）;

y=s;

计算p（x）:

test3.m

x=1.6:

0.2:

2.4;

%x=2的邻域

x=2的邻域：

x

a=[1-18144-6722016-40325376-46082304-512];

p=zeros（1,5）;

fori=1:

5

p（i）=QinJS（a,x（i））;

相应多项式p值：

p

xk=1.95:

0.01:

20.5;

nk=length（xk）;

pk=zeros（1,nk）;

k=1;

nk

pk（k）=QinJS（a,xk（k））;

plot（xk,pk,'

xlabel（'

x'

ylabel（'

p（x）'

x=

1.60001.80002.00002.20002.4000

p=

1.0e-003*

-0.2621-0.000500.00050.2621

p（x）在

[1.95,20.5]上的图像

LU分解，LUDecom.m

function[L,U]=LUDecom（A）

%不带列主元的LU分解

N=size（A）;

n=N

（1）;

L=eye（n）;

U=zeros（n）;

U（1,i）=A（1,i）;

L（i,1）=A（i,1）/U（1,1）;

forj=i:

z=0;

fork=1:

i-1

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

U（i,j）=A（i,j）-z;

forj=i+1:

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

L（j,i）=（A（j,i）-z）/U（i,i）;

PLU分解，PLUDecom.m

function[P,L,U]=PLUDecom（A）

%带列主元的LU分解

[m,m]=size（A）;

U=A;

P=eye（m）;

L=eye（m）;

m

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

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

a=i;

b=abs（t（i））;

ifb<

abs（t（j））

b=abs（t（j））;

a=j;

ifa~=i

forj=1:

c=U（i,j）;

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

U（a,j）=c;

c=P（i,j）;

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

P（a,j）=c;

c=t（a）;

t（a）=t（i）;

t（i）=c;

U（i,i）=t（i）;

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

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

L=tril（U,-1）+eye（m）;

U=triu（U,0）;

（1）

（2）程序：

Test4.m

forn=5:

30

x=zeros（n,1）;

A=-ones（n）;

A（:

n）=ones（n,1）;

fori=1:

A（i,i）=1;

forj=（i+1）:

（n-1）

A（i,j）=0;

x（i）=1/i;

disp（'

当n='

disp（n）;

方程精确解:

x

b=A*x;

%系数b

利用LU分解方程组的解:

[L,U]=LUDecom（A）;

%LU分解

xLU=U\（L\b）

利用PLU分解方程组的解:

[P,L,U]=PLUDecom（A）;

%PLU分解

xPLU=U\（L\（P\b））

%求解A的逆矩阵

A的准确逆矩阵:

InvA=inv（A）

InvAL=zeros（n）;

%利用LU分解求A的逆矩阵

I=eye（n）;

n

InvAL（:

i）=U\（L\I（:

i））;

利用LU分解的A的逆矩阵:

InvAL

End

（1）只列出n=5,6,7的结果

当n=5

x=1.0000

0.5000

0.3333

0.2500

0.2000

xLU=

1.0000

xPLU=

当n=6

0.1667

当n=7

0.1429

（2）只列出n=5,6,7时A的逆矩阵的结果

InvA=

0.5000-0.2500-0.1250-0.0625-0.0625

00.5000-0.2500-0.1250-0.1250

000.5000-0.2500-0.2500

0000.5000-0.5000

0.50000.25000.12500.06250.0625

InvAL=

0.50000.25000.12500.06250.0625

当n=6

0.5000-0.2500-0.1250-0.0625-0.0313-0.0313

00.5000-0.2500-0.1250-0.0625-0.0625

000.5000-0.2500-0.1250-0.1250

0000.5000-0.2500-0.2500

00000.5000-0.5000

0.50000.25000.12500.06250.03130.0313

0.50000.25000.12500.06250.03130.0313

0.5000-0.2500-0.1250-0.0625-0.0313-0.0156-0.0156

00.5000-0.2500-0.1250-0.0625-0.0313-0.0313

000.5000-0.2500-0.1250-0.0625-0.0625

0000.5000-0.2500-0.1250-0.1250

00000.5000-0.2500-0.2500

000000.5000-0.5000

0.50000.25000.12500.06250.03130.01560.0156

0.50000.25000.12500.06250.03130.01560.0156

Cholesky分解：

Cholesky.m

functionL=Cholesky（A）

L=zeros（n）;

L（1,1）=sqrt（A（1,1））;

L（i,1）=A（i,1）/L（1,1）;

forj=2:

s1=0;

j-1

s1=s1+L（j,k）^2;

L（j,j）=sqrt（A（j,j）-s1）;

fori=j+1:

s2=0;

s2=s2+L（i,k）*L（j,k）;

L（i,j）=（A（i,j）-s2）/L（j,j）;

计算Ax=b；

Test5.m

forn=10:

20

A=zeros（n,n）;

b=zeros（n,1）;

A（i,j）=1/（i+j-1）;

b（i,1）=i;

n='

方程组原始解'

x0=A\b

利用Cholesky分解的方程组的解'

L=Cholesky（A）

x=L'

\（L\b）

只列出了n=10,11的结果

n=10

方程组原始解

x0=

1.0e+008*

-0.0000

0.0010

-0.0233

0.2330

-1.2108

3.5947

-6.3233

6.5114

-3.6233

0.8407

利用Cholesky分解的方程组的解

-1.2105

3.5939

-6.3219

6.5100

-3.6225

0.8405

n=

11

1.0e+009*

0.0000

-0.0002

0.0046

-0.0567

0.3687

-1.4039

3.2863

-4.7869

4.2260

-2.0685

0.4305

-0.0563

0.3668

-1.3972

3.2716

-4.7669

4.2094

-2.0608

0.4290

（1）House.m

functionu=House（x）

n=length（x）;

e1=eye（n,1）;

w=x-norm（x,2）*e1;

u=w/norm（w,2）;

（2）Hou_A.m

functionHA=Hou_A（A）

a1=A（:

1）;

n=length（a1）;

w=a1-norm（a1,2）*e1;

H=eye（n）-2*u*u'

HA=H*A;

（3）test6.m

A=[1234;

-12sqrt

（2）sqrt（3）;

-22exp

（1）pi;

-sqrt（10）2-37;

0275/2];

HA=Hou_A（A）

H=

0.2500-0.2500-0.5000-0.79060

-0.25000.9167-0.1667-0.26350

-0.5000-0.16670.6667-0.52700

-0.7906-0.2635-0.52700.16670

00001.0000

HA=

4.0000-2.58111.4090-6.5378

0.00000.47300.8839-1.7805

0.0000-1.05411.6576-3.8836

0.0000-2.8289-4.6770-4.1078

02.00007.00002.5000

Jacobi迭代：

Jaccobi.m

function[x,n]=Jaccobi（A,b,x0）

%--·

方程组系数阵A

方程组右端顶b

%--初始值x0

%--求解要求精确度eps

%--迭代步数控制M

返回求得的解x

返回迭代步数n

M=1000;

eps=1.0e-5;

D=diag（diag（A））;

%求A的对角矩阵

L=-tril（A,-1）;

%求A的下三角阵

U=-triu（A,1）;

%求A的上三角阵

J=D\（L+U）;

f=D\b;

x=J*x0+f

n=1;

%迭代次数

err=norm（x-x0,inf）

while（err>

=eps）

x0=x;

x=J*x0+f

n=n+1;

err=norm（x-x0,inf）

if（n>

=M）

Warning:

迭代次数太多，可能不收敛?

return;

Gauss_Seidel迭代：

Gauss_Seidel.m

function[x,n]=Gauss_Seidel（A,b,x0）

%--Gauss-Seidel迭代法解线性方程组

%--方程组系数阵A

%--方程组右端项b

%--初始值x0

%--求解要求的精确度eps

%--迭代步数控制M

%--返回求得的解x

%--返回迭代步数n

M=10000;

G=（D-L）\U;

f=（D-L）\b;

x=G*x0+f

x=G*x0+f

迭代次数太多，可能不收敛！

解方程组，test7.m

A=[5-1-3;

-124;

-3415];

b=[-2;

1;

10];

精确解'

x=A\b

迭代初始值'

x0=[0;

0;

0]

Jacobi迭代过程：

[xj,nj]=Jaccobi（A,b,x0）;

Jacobi最终迭代结果：

xj

迭代次数'

nj

Gauss-Seidel迭代过程：

[xg,ng]=Gauss_Seidel（A,b,x0）;

Gauss-Seidel最终迭代结果：

xg

ng

精确解

-0.0820

-1.8033

1.1311

迭代初始值

0

-0.4000

0.6667

err=

0.1000

-1.0333

0.4533

1.5333

...

9.6603e-006

xj=

迭代次数

nj=

281

0.3000

0.5067

-0.0360

-0.5313

0.8012

0.8313

-0.0256

-1.1151

0.9589

0.5838

9.4021e-006

xg=

ng=

Newton迭代法：

Newtoniter.m

function[x,iter,fvalue]=Newtoniter（f,df,x0,eps,maxiter）

%牛顿法x得到的近似解

%iter迭代次数

%fvalue函数在x处的值

%f，df被求的非线性方程及导函数

%x0初始值

%eps允许误差限

%maxiter最大迭代次数

fvalue=subs（f,x0）;

dfvalue=subs（df,x0）;

foriter=1:

maxiter

x=x0-fvalue/dfvalue

err=abs（x-x0）

fvalue=subs（f,x0）

dfvalue=subs（df,x0）;

if（err<

eps）||（fvalue==0）,break,end

弦截法：

secant.m

function[x,iter,fvalue]=secant（f,x0,x1,eps,maxiter）

%弦截法x得到的近似解

%f被求的非线性方程

%x0，x1初始值

fvalue0=subs（f,x0）;

fvalue=subs（f,x1）;

x=x1-fvalue*（x1-x0）/（fvalue-fvalue0）

err=abs（x-x1）

x0=x1;

x1=x;

fvalue0=subs（f,x0）;

fvalue=subs（f,x1）

eps）||（fvalue=