收藏
下载资源下载资源 加入VIP,免费下载

矩阵分析与数值分析实验报告.docx

上传人：b****7 文档编号：9732277 上传时间：2023-02-06 格式：DOCX 页数：23 大小：566.44KB
下载 相关 举报
矩阵分析与数值分析实验报告.docx_第1页
第1页 / 共23页
矩阵分析与数值分析实验报告.docx_第2页
第2页 / 共23页
矩阵分析与数值分析实验报告.docx_第3页
第3页 / 共23页
矩阵分析与数值分析实验报告.docx_第4页
第4页 / 共23页
矩阵分析与数值分析实验报告.docx_第5页
第5页 / 共23页
点击查看更多>>
下载资源
资源描述

矩阵分析与数值分析实验报告.docx

《矩阵分析与数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《矩阵分析与数值分析实验报告.docx（23页珍藏版）》请在冰豆网上搜索。

矩阵分析与数值分析实验报告.docx

矩阵分析与数值分析实验报告

《矩阵分析与数值分析》实验报告

 

院系：

姓名：

学号：

所在班号：

任课老师：

 

一．设

，分别编制从小到大和从大到小的顺序程序计算

并指出有效位数。

程序如下：

functionsum3

j=input（'请输入求和个数"j":

'）;

A=0;

B=0;

doubleB;

doubleA;

forn=2:

j

m=n^2-1;

t=1./m;

A=A+t;

end

disp（'从小到大：

'）

s=A

forn=j:

-1:

2

m=n^2-1;

t=1./m;

B=B+t;

end

disp（'从大到小:

'）

s=B

运行结果：

>>sum3

请输入求和个数"j":

100

从小到大：

s=0.740049504950495

从大到小:

s=0.740049504950495

>>sum3

请输入求和个数"j":

10000

从小到大：

s=0.749900004999506

从大到小:

s=0.749900004999500

>>sum3

请输入求和个数"j":

1000000

从小到大：

s=0.749999000000522

从大到小:

s=0.749999000000500

 

二、解线性方程组

1．分别Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组。

迭代法计算停止的条件为：

。

解：

（1）Jacobi迭代法程序代码：

functionjacobi（A,b,N）

clc;clear;

A=[-2100;1-210;01-21;001-2];

b=[-1000]';

N=100;

n=size（A,1）;

D=diag（diag（A））;

L=tril（-A,-1）;

U=triu（-A,1）;

Tj=inv（D）*（L+U）;

cj=inv（D）*b;

tol=1e-06;

k=1;

formatlong

x=zeros（n,1）;

whilek<=N

x（:

k+1）=Tj*x（:

k）+cj;

disp（k）;disp（'x='）;disp（x（:

k+1））;

ifnorm（x（:

k+1）-x（:

k））

disp（'Theprocedurewassuccessful'）

disp（'Condition||x^（k+1）-x^（k）||

disp（k）;disp（'x='）;disp（x（:

k+1））;

break

end

k=k+1;

end

结果输出

Theprocedurewassuccessful

Condition||x^（k+1）-x^（k）||

60

x=

0.79999879906731

0.59999842795870

0.39999805685009

0.199********505

（2）Gauss-Seidel迭代法程序代码：

functiongauss_seidel（A,b,N）

clc;clear;

A=[-2100;1-210;01-21;001-2];

b=[-1000]';

N=100;

n=size（A,1）;

D=diag（diag（A））;

L=tril（-A,-1）;

U=triu（-A,1）;

Tg=inv（D-L）*U;

cg=inv（D-L）*b;

tol=1e-06;

k=1;

x=zeros（n,1）;

whilek<=N

x（:

k+1）=Tg*x（:

k）+cg;

disp（k）;disp（'x='）;disp（x（:

k+1））;

ifnorm（x（:

k+1）-x（:

k））

disp（'Theprocedurewassuccessful'）

disp（'Condition||x^（k+1）-x^（k）||

disp（k）;disp（'x='）;disp（x（:

k+1））;

break

end

k=k+1;

end

结果输出

Theprocedurewassuccessful

Condition||x^（k+1）-x^（k）||

31

x=

0.79999921397935

0.59999897108561

0.39999916759077

0.199********539

 

2.用Gauss列主元消去法、QR方法求解如下方程组：

（1）Gauss列主元消去法

程序代码：

functionx=Gaussmain（A,b）

clc;clear;

formatlong

A=[1212;253-2;-2-235;1323];

b=[47-10]';

N=length（A）;

x=zeros（N,1）;

y=zeros（N,1）;

c=0;

d=0;

A（:

N+1）=b;

fork=1:

N-1

fori=k:

4

ifc

d=i;c=abs（A（i,k））;

end

end

y=A（k,:

）;

A（k,:

）=A（d,:

）;

A（d,:

）=y;

fori=k+1:

N

c=A（i,k）;

forj=1:

N+1

A（i,j）=A（i,j）-A（k,j）*c/A（k,k）;

end

end

end

b=A（:

N+1）;

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

fork=N-1:

-1:

1

x（k）=（b（k）-A（k,k+1:

N）*x（k+1:

N））/A（k,k）;

end

结果输出

ans=

18.00000000000000

-9.57142857142857

6.00000000000000

-0.42857142857143

（2）QR方法：

程序代码

functionQR（A,b）

clc;clear;

formatlong

A=[1212;253-2;-2-235;1323];

b=[47-10]';

[Q,R]=qr（A）

y=Q'*b;

x=R\y

结果输出

Q=

-0.31622776601684-0.04116934847963-0.751646028002830.57735026918962

-0.63245553203368-0.49403218175557-0.150********056-0.57735026918963

0.63245553203368-0.74104827263336-0.22549380840085-0.00000000000000

-0.31622776601684-0.452862833275940.601316822402260.57735026918963

R=

-3.16227766016838-6.00832755431992-0.948683298050512.84604989415154

0-2.42899156029822-4.65213637819829-4.158********272

00-0.67648142520255-0.52615221960200

.0414********

x=

17.99999999999989

-9.57142857142851

5.99999999999997

-0.42857142857143

三、非线性方程的迭代解法

1．用Newton迭代法求方程

的根，计算停止的条件为：

；

编程如下：

functionnewton（f,df,x,a,a0）

symsx

f=input（'pleaseenteryourequation:

'）

a0=input（'pleaseenteryoux（0）:

'）;

df=diff（f）

e=1e-6;

a1=a0+1;

N=0;

whileabs（a1-a0）>e

a=a0-subs（f,a0）/subs（df,a0）;

a1=a0;

a0=a;

N=N+1;

end

fprintf（'a=%0.6f',a）

N

运行结果：

>>newton

pleaseenteryourequation:

exp（x）+2^（-x）+2*cos（x）-6

f=exp（x）+2^（-x）+2*cos（x）-6

pleaseenteryoux（0）:

2

df=exp（x）-2^（-x）*log

（2）-2*sin（x）

a=1.829384

N=4

2．利用Newton迭代法求多项式

的所有实零点，注意重根的问题。

由于不知道重根的个数，采用试探法求重根。

function[X]=newton2（）

clc;clear;

symsx;

delta=1e-06;

f=inline（'x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951'）;

df=inline（'4*x^3-16.2*x^2+21.12*x-8.954'）;

u=inline（'（x^4-5.4*x^3+10.56*x^2-8.954*x+2.7951）/（4*x^3-16.2*x^2+21.12*x-8.954）'）;

du=inline（'1-（x^4-27/5*x^3+264/25*x^2-4477/500*x+27951/10000）/（4*x^3-81/5*x^2+528/25*x-4477/500）^2*（12*x^2-162/5*x+528/25）'）;

j=1;

fori=0:

1:

3

x0=i;

while1

x1=x0-feval（u,x0）/feval（du,x0）;

err=abs（x1-x0）;

x0=x1;

y=feval（f,x0）;

if（err

X（j）=x1;

j=j+1;

break

end

end

end

X

%运行结果：

%X=1.10001.10002.10001.1000

四、数值积分

分别用复化梯形公式和Romberg公式计算积分

的近似值，要求误差不超过

，并给出节点个数。

1）复化梯形公式

程序代码：

function[I,step]=Trapezoid（f,a,b,eps）

f='1/x';a=2;b=8;eps=1.0e-5;

n=1;

h=（b-a）/2;

I1=2;

I2=（subs（sym（f）,findsym（sym（f））,a）+subs（sym（f）,findsym（sym（f））,b））/h;

whileabs（I2-I1）>eps

n=n+1;

h=（b-a）/n;

I1=I2;

I2=0;

fori=0:

n-1

x=a+h*i;

x1=x+h;

I2=I2+（h/2）*（subs（sym（f）,findsym（sym（f））,x）+...

subs（sym（f）,findsym（sym（f））,x1））;

end

end

I=I2

step=n

结果输出

I=

1.38654458753674

step=

53

2）Romberg公式

程序代码：

functionapprox=romberg（f,a,b,TOL）

clc;clear;

f=inline（'1/x'）;

a=2;b=8;

TOL=1e-05;

A=zeros（20,20）;

A（1,1）=（b-a）*（feval（f,a）+feval（f,b））/2;

h=（b-a）/2;

A（2,1）=A（1,1）/2+h*feval（f,a+h）;

A（2,2）=（4*A（2,1）-A（1,1））/3;

errest=abs（A（2,2）-A（1,1）/2）;

i=2;

while（errest>TOL）

i=i+1;

h=h/2;

sum=0.0;

forj=1:

2:

2^（i-1）-1

sum=sum+feval（f,a+j*h）;

end;

A（i,1）=A（i-1,1）/2+h*sum;

forj=2:

i

power=4^（j-1）;

A（i,j）=（power*A（i,j-1）-A（i-1,j-1））/（power-1）;

end;

errest=abs（A（i,i）-A（i-1,i-1））/2^（i-1）;

end;

if（nargout==0）

s=sprintf（'\t\tapproximatevalueofintegral:

\t%.12f\n',A（i,i））;

s=sprintf（'%s\t\terrorestimate:

\t\t\t\t\t%.4e\n',s,errest）;

s=sprintf（'%s\t\tnumberoffunctionevaluations:

\t%d\n',s,2^（i-1）+1）;

disp（s）

else

approx=A（i,i）;

end

运行结果：

approximatevalueofintegral:

1.386297441871

errorestimate:

8.7527e-006

numberoffunctionevaluations:

17

五、插值与逼近

1．给定

上的函数

，请做如下的插值逼近：

构造等距节点分别取

，

，

的Lagrange插值多项式；

构造分段线性取

的Lagrange插值多项式；

取Chebyshev多项式

的零点：

，

作插值节点构造

的插值多项式

和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。

程序代码

functionLagrange

clc;clear;closeall;

fori=1:

3;

ifi==1

N=5;

elseifi==2

N=8;

else

N=10;

end

f=inline（'1/（1+25*x^2）'）;

x1=zeros（1,N+1）;

a=-1;

b=1;

fori=1:

N+1

x1（i）=a+（i-1）*（b-a）/N;

y1（i）=feval（f,x1（i））;

end

symsx

ff=0;

fori=1:

N+1

f=1;

forj=1:

i-1

f=f*（x-x1（j））/（x1（i）-x1（j））;

end

forj=i+1:

N+1

f=f*（x-x1（j））/（x1（i）-x1（j））;

end

ff=f*y1（i）+ff;

f=1;

end

ff=collect（ff,x）;

ff=vpa（ff,4）;

y=ff;

p=ezplot（y,[a,b]）;grid

YLIM（[-0.10.6]）;

ifN==5

set（p,'Color','black'）;

set（p,'LineStyle','--'）;

lagrange_5=y

elseifN==8

set（p,'Color','r'）;

set（p,'LineStyle','--'）;

lagrange_8=y

else

set（p,'Color','b'）;

set（p,'LineStyle','--'）

lagrange_10=y

end

holdon;

xlabel（'x'）;ylabel（'y'）;

title（'y=p（x）'）;holdon;

end

Lag_Cheb（）;

x=-1:

0.01:

1;

y=1./（1+25*x.^2）;

acu=plot（x,y）;gridon;holdon

set（acu,'Color','m'）;

set（acu,'LineStyle','-'）;

legend（'N=5','N=8','N=10','Cheb,N=10','¾«È·Öµ'）;

 

functionLag_Cheb（）

f=inline（'1/（1+25*x^2）'）;

N=10;

x1=zeros（1,N+1）;

a=-1;

b=1;

fori=1:

N+1

x1（i）=cos（（2*i-1）*pi/（2*N））;

y1（i）=feval（f,x1（i））;

end

symsx

ff=0;

fori=1:

N+1

f=1;

forj=1:

i-1

f=f*（x-x1（j））/（x1（i）-x1（j））;

end

forj=i+1:

N+1

f=f*（x-x1（j））/（x1（i）-x1（j））;

end

ff=f*y1（i）+ff;

f=1;

end

ff=collect（ff,x）;

ff=vpa（ff,4）;

yy=ff;

ff=collect（ff,x）;

yy=ff;

lagrange_chebshev_10=yy

cheb=ezplot（yy,[a,b]）;gridon

YLIM（[-0.10.6]）;

set（cheb,'Color','g'）;

set（cheb,'LineStyle',':

'）;

 

结果输出

lagrange_5=.5673+1.202*x^4-1.731*x^2

lagrange_8=1.+53.69*x^8-102.8*x^6+61.37*x^4-13.20*x^2

lagrange_10=1.-220.9*x^10+494.9*x^8-381.4*x^6+123.4*x^4-16.86*x^2

lagrange_chebshev_10=.7413-5.359*x^10-.4e-2*x^9+18.96*x^8-.1321*x^7-

25.78*x^6+.10*x^5+16.81*x^4+.5e-2*x^3-5.336*x^2+.5288e-3*x

2．已知函数值

0

1

2

3

4

5

6

7

8

9

10

2.51

3.30

4.04

4.70

5.22

5.54

5.78

5.40

5.57

5.70

5.80

和边界条件：

．　求三次样条插值函数

并画出其图形．

程序代码：

functionSpline3_1（x,y,df0,dfn）

formatshort

x=[012345678910];

y=[2.513.304.044.705.225.545.785.405.575.705.80];

plot（x,y,'g*','MarkerSize',15）;holdon;

df0=1;

dfn=0;

n=length（x）;

h=zeros（n-1,1）;

lan=zeros（n-2,1）;

mu=zeros（n-2,1）;

g=zeros（n-2,1）;

m=zeros（n,1）;

m

（1）=df0;

m（n）=dfn;

fori=1:

n-1

h（i）=x（i+1）-x（i）;

end

fori=1:

n-2

mu（i）=h（i）/（h（i+1）+h（i））;

lan（i）=h（i+1）/（h（i+1）+h（i））;

g（i）=3*（mu（i）*（y（i+2）-y（i+1））/h（i+1）+lan（i）*（y（i+1）-y（i））/h（i））;

end

A=zeros（n-2,n-2）;

A（1,1）=2;

A（1,2）=mu

（1）;

A（n-2,n-2）=lan（n-2）;

fori=2:

n-2

A（i,i）=2;

A（i,i-1）=mu（i）;

A（i-1,i）=lan（i）;

end

g

（1）=g

（1）-lan

（1）*df0;

g（n-2）=g（n-2）-mu（n-2）*dfn;

b=A\g;

fori=2:

n-1

m（i）=b（i-1）;

end

symsz

fori=1:

n-1

xx=x（i）:

0.01:

x（i+1）;

sx1=y（i）*（xx-x（i+1））.^2.*（h（i）+2*（xx-x（i）））/h（i）^3;

sx2=y（i+1）*（xx-x（i））.^2.*（h（i）-2*（xx-x（i+1）））/h（i）^3;

sx3=m（i）*（xx-x（i+1））.^2.*（xx-x（i））/h（i）^2;

sx4=m（i+1）*（xx-x（i））.^2.*（xx-x（i+1））/h（i）^2;

sx=sx1+sx2+sx3+sx4;

z1=y（i）*（z-x（i+1））.^2.*（h（i）+2*（z-x（i）））/h（i）^3;

z2=y（i+1）*（z-x（i））.^2.*（h（i）-2*（z-x（i+1）））/h（i）^3;

z3=m（i）*（z-x（i+1））.^2.*（z-x（i））/h（i）^2;

z4=m（i+1）*（z-x（i））.^2.*（z-x（i+1））/h（i）^2;

zz=z1+z2+z3+z4;

zz=collect（zz）;

vpa（zz,4）

plot（xx,sx,'r'）;

holdon;

gridon;

end

程序输出

ans=2.510+.1379*z^3-.3479*z^2+z

ans=2.692-.4369e-1*z^3+.1969*z^2+.4552*z

ans=2.287+.6872e-2*z^3-.1065*z^2+1.062*z

ans=3.655-.4380e-1*z^3+.3495*z^2-.3060*z

ans=-6.080+.1083*z^3-1.476*z^2+6.995*z

ans=41.14-.2694*z^3+4.190*z^2-21.34*z

ans=-109.8+.4295*z^3-8.390*z^2+54.14*z

ans=133.0-.2784*z^3+6.475*z^2-49.91*z

ans=-57.75+.9411e-1*

展开阅读全文
相关资源
猜你喜欢
相关搜索

当前位置：首页 > 总结汇报 > 学习总结

copyright@ 2008-2022 冰豆网网站版权所有

经营许可证编号:鄂ICP备2022015515号-1