当前位置:文档之家› 清华大学高等数值分析_第二次实验作业

清华大学高等数值分析_第二次实验作业

清华大学高等数值分析_第二次实验作业
清华大学高等数值分析_第二次实验作业

清华大学高等数值分析课程作业

第二次实验 作业

第一题:构造例子特征值全部在右半平面时,观察基本的Arnoldi 方法和GMRES 方法的数值性态,和相应重新启动算法的收敛性。 答:

1、计算初始条件

1) 矩阵A 的生成

根据实Schur 分解,构造矩阵如下形式

11112222/2/2/2/2n n

A n n n n ?-?? ? ? ?- ?= ? ? ?

- ?

???

其中,A 由n/2个块形成,每个对角块具有如下形式,对应一对特征向量i i i αβ+

i

i i i A α

ββα-??

= ???

、 这里,取n=1000,得到矩阵A 。经过验证,A 的特征值分布均在右半平面,如下图所示

50

100

150

200250300

350

400

450

500

-500-400-300-200-1000100200

300400500

复平面中A 的特征值分布情况

实部 Im(x)

虚部 R e (x )

特征值

2) b 的初值为 b=(1,1,1,1..1)T 3) 迭代初值为 x 0=0 4) 停机准则为 ε=10-6

2、基本的Arnoldi 和GMRES 方法

代入前面提到的初始值A 、b 、x0,得到的收敛结果如下

100200

300400500600

10-7

10

-6

10

-5

10

-4

10

-3

10

-2

10

-1

10

两种基本算法的||r k ||收敛曲线 (阶数n=1000)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

结果讨论:从图中可以看出,基本的Arnoldi 方法经过554步收敛,基本的GMRES 方法经过535步收敛。这是由于GMRES 具有残差最优性,会略快于Arnoldi 方法,但是,由于两种方法的基本原理近似,GMRES 方法不会实质性的提速。

此外,从收敛曲线上看,由于特征值均处在右半平面,收敛曲线平滑,收敛速度(收敛因子)比较均匀。

3、重启动的GMRES 和Arnoldi 算法

对上述A 、b 、x0使用重启动的Arnoldi 和GMRES 算法。变化m 经过多次实验,得到如下结果: 1. 总迭代步长与m 之间的关系

50

100

150

9001000110012001300140015001600

170018001900两种重启算法总 迭代步长 与m 的关系曲线比较

总迭代次数

m

重启GMRES 重启Aronldi

从表中数据可以看出,结果相差并不大

2. 重启次数与m 之间的关系

50

100150

050

100

150

200

250

两种重启算法 重启次数 与m 的关系曲线比较

重启次数

m

重启GMRES 重启Aronldi

3. 计算时间与m 之间的关系

050

100150

3

3.5

4

4.5

5

5.5

重启的Arnoldi 算法 求解运行时间 与m 的关系曲线

求解运行时间(s )

m

50

100150

34

5

6

7

8

9

10

重启的GMRES 算法 求解运行时间 与m 的关系曲线

求解运行时间(s )

m

从上述结果中看到在m=40左右时,同时有求解时间和迭代总步长最小的结果。当m=40时,得到的收敛曲线结果如下图所示:

510

1520253035

10-7

10

-6

10

-5

10

-4

10

-3

10

-2

10

-1

10

两种重启算法的||r k ||收敛曲线 (m=40)

重启次数

||r k ||/||b ||

重启GMRES 重启Aronldi

200400

600800100012001400

10-7

10

-6

10

-5

10

-4

10

-3

10

-2

10

-1

10

两种重启算法的||r k ||收敛曲线 (m=40)

总迭代次数

||r k ||/||b ||

重启GMRES 重启Aronldi

160180200220

240260280300320340

10

-2

两种重启算法的||r k ||收敛曲线 (m=40)

总迭代次数

||r k ||/||b ||

重启GMRES 重启Aronldi

4. 结论

从上述结果中分析可得如下结论:

1) 对于这两种算法来说,重启动算法的总迭代步长要大于基本的算法,这个很容易理解,这

是由于在重启动时,会丢掉一些之前算过的信息,从另一个点我们认为更好(其实不一定好)的初值x0重新开始计算。

2) 这两种重启动算法,在两次重启动的交接点处都会产生不连续。GMRES 在两次重启动之间不

满足残差单调不增,然而在单次重启的内部计算过程中,则满足单调不增的原则。Arnodli 过程由于不具有残差二范数最小的性质,所以在以上两个位置处都可能会发生跳动。 3) 一般来说,重启次数随着m 的增大不断减小,当m 非常大的时候,就过度到了基本的不启动

算法。

4) 随着m 的增大,迭代总次数并不是单调下降的。重新启动的Arnoldi 方法在m=40时迭代次数

最少,重新启动的GMRES 方法m=30或40时迭代次数最少。

5) 计算代价(用运算时间估算)与m 不成线性关系,存在一个最优值,当m 大到一定程度,虽

然重启次数很小,但是每步代价很大,随着m 的进一步增加,代价逐渐增大。 6) 重启的算法如果收敛的话,代价可能会远远小于不重启的算法。但是最合适的m 的选取因问

题不同而不同,不好把握。

7) 对于同样的矩阵A ,GMRES 方法的收敛速度一般比相应的Arnoldi 方法快。这是由于其残差具

有最优性,而且在迭代步数相同的情况下,GMRES 方法的相对残差比相应的Arnoldi 方法的相对残差要小。

第二题。对于1中的矩阵,将特征值进行平移,使得实部有正有负,和1的结果进行比较,方法的收敛速度会如何?基本的Arnoldi 算法有无峰点?若有,基本的GMRES 算法相应地会怎样? 答:

1. 使用如下方法对第一题中的A 矩阵进行平移

'A A I μ=-

则得到的矩阵的特征值为'

i i λλμ=-。

2. 对第一题的A ,分别取μ=1.5、2.5、5.5、20.5、50.5、100.5、200.5、400.5、600.5、900.5、950.5,得到的

负特征值分别为1、2、5、20、50、100、200、400、600、900、950个。分别求解上述矩阵,得到的结果如下

100

200

300400500

600

700

10-8

10

-6

10

-4

10

-2

10

10

2

两种基本算法的||r k ||收敛曲线 (负特征值个数为1个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

100200

300400500600700

10-8

10

-6

10

-4

10

-2

10

10

2

两种基本算法的||r k ||收敛曲线 (负特征值个数为2个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

100

200

300400500

600

700

10-6

10

-4

10

-2

10

10

2

10

4

两种基本算法的||r k ||收敛曲线 (负特征值个数为5个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

100200

300400500600700

10-8

10

-6

10

-4

10

-2

10

10

2

两种基本算法的||r k ||收敛曲线 (负特征值个数为20个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

100

200

300

400500600

700

800

900

10-8

10

-6

10

-4

10

-2

10

10

2

两种基本算法的||r k ||收敛曲线 (负特征值个数为50个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

100200300

400500600700800

10-8

10

-6

10

-4

10

-2

10

10

2

10

4

两种基本算法的||r k ||收敛曲线 (负特征值个数为100个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

100

200

300

400500600

700

800

900

10-8

10

-6

10

-4

10

-2

10

10

2

10

4

两种基本算法的||r k ||收敛曲线 (负特征值个数为200个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

50

100150

200250

10-8

10

-6

10

-4

10

-2

10

10

2

10

4

两种基本算法的||r k ||收敛曲线 (负特征值个数为400个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

20

40

60

80

100120140

160

180

200

10-7

10

-610-510-410-3

10

-210-1100101

两种基本算法的||r k ||收敛曲线 (负特征值个数为600个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

204060

80100120140160180

10-7

10

-6

10-5

10-4

10-3

10

-2

10-1

100

101

两种基本算法的||r k ||收敛曲线 (负特征值个数为900个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

20

40

60

80100120

140

160

180

10-7

10

-610-510-410-3

10

-210-1100101

两种基本算法的||r k ||收敛曲线 (负特征值个数为950个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

204060

80100120140160180

10-7

10

-6

10-5

10-4

10-3

10

-2

10-1

100

101

两种基本算法的||r k ||收敛曲线 (负特征值个数为1000个)

迭代次数

||r k ||/||b ||

基本的Arnoldi 算法基本的GMRES 算法

当负特征值个数逐渐增多时,迭代次数与负特征值个数m 的关系曲线如下图

100200300400

5006007008009001000

100200300400500600

700800900重启的Arnoldi 算法总 迭代步长 与m 的关系曲线

总迭代次数

m

基本的Arnoldi 算法基本的GMRES 算法

注:上图标题有误,应为 两种基本算法迭代步长与负特征值个数m 的关系曲线

3. 结果讨论

1) 现象:随着负特征值的增多,Arnoldi 过程跳动越严重,但是整体是收敛的!这是由于负特征值的存

在,使得海森贝格矩阵H 发生近似奇异而发生近似中断而引起的。然而,GMRES 的残差始终单调不减,当Aronldi 出现尖峰时,GMRES 的残差不变。

2) 当负半平面特征值个数比较少时,迭代次数明显变多,并在负半平面特征值个数为200左右时,迭代

次数最多,耗时最长,收敛速度最慢; 3) 当负半平面特征值个数变得比较多时,迭代次数减少,并且将比第一题的收敛速度更快得多,只需200

多步迭代即可收敛。 4) 针对选取A 的特殊性,

第三题、对1中的例子固定特征值的实部,变化虚部,比较收敛性。 答:

1) 首先,固定A 的特征值的实部,变化虚部。

由前面可知,A 的表达式为

11112222/2/2/2/2n n

A n n n n ?-?? ? ? ?- ?= ? ? ?

- ?

???

其中,A 由n/2个块形成,每个对角块具有如下形式,对应一对特征向量i i i αβ+

i

i i

i A α

ββα-??

= ???

将上式中的

'i i k ββ=

改变k 的大小,可以增加虚部的大小。

本次实验,取k 分别为0.5、1、1.5、2.5、4、5,分别得到的收敛曲线如下

由前面可知,A 的表达式为 使用如下方法

则可以在A 的基础上固定A 的实特征值,而增加

Matlab 程序

基本的Arnoldi 过程

%Arnold Method

%By liu libin @2011-12-2

%Email:llb11@https://www.doczj.com/doc/f819192466.html,

function [ x,Errs,flag]=Arnoldllb(A,b,x0,tol,mtime,cl) flag=0; H=[]; Errs=[]; r0=b-A*x0; nr0=norm(r0,2); V=r0/nr0; [nnn,n]=size(A); for j=1:mtime

%disp(['k=',num2str(j)]) %clc

%disp([num2str(100*j/mtime),'%'])

w=A*V(:,j);

for i=1:j

H(i,j)=V(:,i)'*w;

w=w-H(i,j)*V(:,i);

end

H(j+1,j)=norm(w,2) ;

if H(j+1,j)<1e-13

disp('3????D??')

end

V=[V,w/H(j+1,j)];

if j>=2

%**********devide H=QR by Givens Method

%{

gm=nr0*[1;zeros(j-1,1)];

Rm=H(1:j,1:j);

Qm=eye(j,j);

for iqr=1:j-1

Ci=Rm( iqr ,iqr)/((Rm(iqr,iqr)^2+Rm(iqr+1,iqr)^2)^0.5) ; Si=Rm(iqr+1,iqr)/((Rm(iqr,iqr)^2+Rm(iqr+1,iqr)^2)^0.5) ; Rm_temp=Rm;

Rm_temp( iqr ,:)= Ci*Rm(iqr,:)+Si*Rm(iqr+1,:);

Rm_temp(iqr+1,:)=-Si*Rm(iqr,:)+Ci*Rm(iqr+1,:);

Rm=Rm_temp;

%compute Qm

Qm_temp=Qm;

Qm_temp( iqr ,:)= Ci*Qm(iqr,:)+Si*Qm(iqr+1,:);

Qm_temp(iqr+1,:)=-Si*Qm(iqr,:)+Ci*Qm(iqr+1,:);

Qm=Qm_temp;

gm_temp=gm;

gm_temp( iqr )= Ci*gm(iqr)+Si*gm(iqr+1);

gm_temp(iqr+1)=-Si*gm(iqr)+Ci*gm(iqr+1);

gm=gm_temp;

end

ym=zeros(j,1);

ym( j )=gm( j )/Rm(j,j);

ym(j-1)=(gm(j-1)-ym(j)*Rm(j-1,j))/Rm(j-1,j-1);

for is=j-2:-1:1

ym(is)=(gm(is)-Rm(is,is+1:j)*ym(is+1:j))/Rm(is,is);

end %for is

%}

gm=nr0*[1;zeros(j-1,1)];

ym=H(1:j,1:j)\gm;

x=x0+V(:,1:j)*ym;

Errm=abs(H(j+1,j)*ym(j))/(norm(r0,2));

Errs=[Errs,Errm];

if cl==1

if Errm

disp('*******Converge!*******')

flag=1;

break;

end%for Err

end

end%for j>2

end% for j

Errs=[1,Errs];

%ok

end%for func

基本的GMRES过程

%GMRES Method

%By liu libin @2011-12-2

%Email:llb11@https://www.doczj.com/doc/f819192466.html,

function [x,Errs,flag,name]=GMRESllb(A,b,x0 ,Err,mtime,InterEnabled)

name='GMRES';

flag=0;

H=[];

Errs=[];

r0=b-A*x0;

nr0=norm(r0,2);

V=r0/nr0;

[nnn,n]=size(A);

for j=1:mtime

disp([num2str(100*j/mtime),'%'])

w=A*V(:,j) ;

for i=1:j

H(i,j)=V(:,i)'*w;

w=w-H(i,j)*V(:,i) ;

end

H(j+1,j)=norm(w,2);

V=[V,w/H(j+1,j)];

if j>=2

%**********devide H=QR by Givens Method

Rm=H(1:j+1,1:j);

Qm=eye(j+1,j+1);

for iqr=1:j

Ci=Rm( iqr ,iqr)/((Rm(iqr,iqr)^2+Rm(iqr+1,iqr)^2)^0.5) ;

Si=Rm(iqr+1,iqr)/((Rm(iqr,iqr)^2+Rm(iqr+1,iqr)^2)^0.5) ;

Rm_temp=Rm;

Rm_temp( iqr ,:)= Ci*Rm(iqr,:)+Si*Rm(iqr+1,:);

Rm_temp(iqr+1,:)=-Si*Rm(iqr,:)+Ci*Rm(iqr+1,:);

Rm=Rm_temp;

%compute Qm

Qm_temp=Qm;

Qm_temp( iqr ,:)= Ci*Qm(iqr,:)+Si*Qm(iqr+1,:);

Qm_temp(iqr+1,:)=-Si*Qm(iqr,:)+Ci*Qm(iqr+1,:);

Qm=Qm_temp;

end

%{

gm=Qm(1:j,1)*norm(r0,2);

ym=zeros(j ,1);

ym( j )=gm( j )/Rm(j,j);

ym(j-1)=(gm(j-1)-ym(j)*Rm(j-1,j))/Rm(j-1,j-1);

for is=j-2:-1:1

ym(is)=(gm(is)-Rm(is,is+1:1:j)*ym(is+1:j))/Rm(is,is);

end %for is

%}

ym=Rm(1:j,1:j)\Qm(1:j,1)*norm(r0,2);

x=x0+V(:,1:j)*ym ;

Errm=abs(Qm(j+1,1))*norm(r0,2)/norm(b,2);

Errs=[Errs,Errm];

if InterEnabled==1

if Errm

disp('*******Converge!*******')

flag=1;

break;

end%for Err

end

%Rm*ym-gm

end%for j>2

end% for j

%ok

end%for func

重启算法

重启的Arnoldi算法

%Arnold Method

%By liu libin @2011-12-2

%Email:llb11@https://www.doczj.com/doc/f819192466.html,

function [x,Errs,rkSteps,flag]=ArnoldMllb(A,b,x0,tol,mtime,MaxM)

flag=0;

H=[];

Errs=[];

Maxi=1;r0=b-A*x0;

rkSteps=[1];

normr0=(norm(A,1)*norm(x0)+norm(b));%norm(r0,2);

while flag==0

H=[];

r0=b-A*x0;

nr0=norm(r0,2);

V=r0/nr0;

[nnn,n]=size(A);

for j=1:mtime

%disp(['k=',num2str(j)])

%clc

%disp([num2str(100*j/mtime),'%'])

w=A*V(:,j);

for i=1:j

H(i,j)=V(:,i)'*w;

w=w-H(i,j)*V(:,i);

end

H(j+1,j)=norm(w,2) ;

if H(j+1,j)<1e-13

disp('3????D??')

end

V=[V,w/H(j+1,j)];

if j>=2

gm=nr0*[1;zeros(j-1,1)];

ym=H(1:j,1:j)\gm;

x=x0+V(:,1:j)*ym;

Errm=abs(H(j+1,j)*ym(j))/norm(b);

Errs=[Errs,Errm];

if Errm

disp('******Arnoldi m·收敛!*******') flag=1;

break;

end%for Err

end%for j>2

end% for j

x0=x;

rkSteps=[rkSteps,Errm];

if Maxi==MaxM

break;

end

Maxi=Maxi+1 ;

end%for while

%ok

end%for func

重启的GMRES算法

%GMRES Method

%By liu libin @2011-12-2

%Email:llb11@https://www.doczj.com/doc/f819192466.html,

function [x,Errs,rkSteps,flag]=GMRESMllb(A,b,x0 ,Err,mtime,MaxM) flag=0;

H=[];

Errs=[];

rkSteps=[];

while flag==0

r0=b-A*x0;

nr0=norm(r0,2);

V=r0/nr0;

[nnn,n]=size(A);

for j=1:mtime

%disp([num2str(100*j/mtime),'%'])

w=A*V(:,j) ;

for i=1:j

H(i,j)=V(:,i)'*w;

w=w-H(i,j)*V(:,i) ;

end

H(j+1,j)=norm(w,2);

V=[V,w/H(j+1,j)];

if j>=2

%**********devide H=QR by Givens Method

Rm=H(1:j+1,1:j);

Qm=eye(j+1,j+1);

for iqr=1:j

Ci=Rm( iqr ,iqr)/((Rm(iqr,iqr)^2+Rm(iqr+1,iqr)^2)^0.5) ;

Si=Rm(iqr+1,iqr)/((Rm(iqr,iqr)^2+Rm(iqr+1,iqr)^2)^0.5) ;

Rm_temp=Rm;

Rm_temp( iqr ,:)= Ci*Rm(iqr,:)+Si*Rm(iqr+1,:);

Rm_temp(iqr+1,:)=-Si*Rm(iqr,:)+Ci*Rm(iqr+1,:);

Rm=Rm_temp;

%compute Qm

Qm_temp=Qm;

Qm_temp( iqr ,:)= Ci*Qm(iqr,:)+Si*Qm(iqr+1,:);

Qm_temp(iqr+1,:)=-Si*Qm(iqr,:)+Ci*Qm(iqr+1,:);

Qm=Qm_temp;

end

%{

gm=Qm(1:j,1)*norm(r0,2);

ym=zeros(j ,1);

ym( j )=gm( j )/Rm(j,j);

ym(j-1)=(gm(j-1)-ym(j)*Rm(j-1,j))/Rm(j-1,j-1);

for is=j-2:-1:1

ym(is)=(gm(is)-Rm(is,is+1:1:j)*ym(is+1:j))/Rm(is,is); end %for is

%}

ym=Rm(1:j,1:j)\Qm(1:j,1)*norm(r0,2);

x=x0+V(:,1:j)*ym ;

Errm=abs(Qm(j+1,1))*norm(r0,2)/norm(b);

Errs=[Errs,Errm];

if Errm

disp('*******????????~~~~·?·¨ê?á2£?*******')

flag=1;

break;

end%for Err

%Rm*ym-gm

end%for j>2

end% for j

x0=x;

rkSteps=[rkSteps,Errm];

end

%ok

end%for func

重启算法调用函数扫描m

clc

global A

global b

x0=zeros(length(b),1);

mtimeAll=[10:10:150];

%mtimeAll=40;

tol=1e-6;

MaxM=300;

GMK=[];

ARK=[];

GMKAll=[];

ARKAll=[];

ARTime=[];

GMTime=[];

for i=1:length(mtimeAll)

mtime=mtimeAll(i);

tic

[x,Errs,rkSteps,flag]=ArnoldMllb(A,b,x0,tol,mtime,MaxM);

t1=toc;

ARTime=[ARTime,t1];

ARKAll=[ARKAll,length(Errs)];

ARK=[ARK,length(rkSteps)];

tic

figure(4)

semilogy([1:length(Errs)],Errs,'b.-')

stringsCG=['重启Arnodli算法的||r_k||收敛曲线 (m=']; stringsCG=[stringsCG,num2str(mtime),')'];

title(stringsCG)

xlabel('总迭代次数')

ylabel('||r_k||/||b||')

hold on

grid on

figure(5)

semilogy([1:length(rkSteps)],rkSteps,'b.-')

stringsCG=['重启Arnodli算法的||r_k||收敛曲线 (m=']; stringsCG=[stringsCG,num2str(mtime),')'];

title(stringsCG)

xlabel('重启次数')

ylabel('||r_k||/||b||')

hold on

grid on

[x,Errs,rkSteps,flag]=GMRESMllb(A,b,x0,tol,mtime,MaxM); t2=toc;

GMTime=[GMTime,t2];

GMKAll=[GMKAll,length(Errs)];

GMK=[GMK,length(rkSteps)];

i/length(mtimeAll)*100

figure(4)

semilogy([1:length(Errs)],Errs,'r.-')

stringsCG=['两种重启算法的||r_k||收敛曲线 (m='];

stringsCG=[stringsCG,num2str(mtime),')'];

title(stringsCG)

xlabel('总迭代次数')

ylabel('||r_k||/||b||')

legend('重启GMRES','重启Aronldi')

hold on

grid on

figure(5)

semilogy([1:length(rkSteps)],rkSteps,'r.-')

stringsCG=['两种重启算法的||r_k||收敛曲线 (m='];

stringsCG=[stringsCG,num2str(mtime),')'];

title(stringsCG)

xlabel('重启次数')

ylabel('||r_k||/||b||')

legend('重启GMRES','重启Aronldi')

hold on

grid on

end

figure(12)

plot(mtimeAll,ARKAll,'bo-','MarkerFaceColor','b') title('重启的Arnoldi算法总迭代步长与m的关系曲线') ylabel('总迭代次数')

xlabel('m')

grid on;

figure(11)

plot(mtimeAll,ARK,'bo-','MarkerFaceColor','b') title('重启的Arnoldi算法重启次数与m的关系曲线')

ylabel('重启次数')

xlabel('m')

grid on;

figure(16)

plot(mtimeAll,GMKAll,'ro-','MarkerFaceColor','r') title('重启的GMRES算法总迭代步长与m的关系曲线')

ylabel('总迭代次数')

xlabel('m')

grid on;

figure(17)

plot(mtimeAll,GMK,'ro-','MarkerFaceColor','r') title('重启的GMRES算法重启次数与m的关系曲线')

ylabel('重启次数')

xlabel('m')

grid on;

figure(18)

plot(mtimeAll,GMKAll,'ro-','MarkerFaceColor','r') hold on

plot(mtimeAll,ARKAll,'b^-','MarkerFaceColor','b') title('两种重启算法总迭代步长与m的关系曲线比较')

ylabel('总迭代次数')

xlabel('m')

legend('重启GMRES','重启Aronldi')

grid on;

figure(19)

plot(mtimeAll,GMK,'ro-','MarkerFaceColor','r') hold on

plot(mtimeAll,ARK,'b^-','MarkerFaceColor','b')

title('两种重启算法重启次数与m的关系曲线比较')

ylabel('重启次数')

legend('重启GMRES','重启Aronldi')

xlabel('m')

grid on;

%时间比较

figure(20)

plot(mtimeAll,GMTime,'bo-','MarkerFaceColor','b') title('重启的GMRES算法求解运行时间与m的关系曲线') ylabel('求解运行时间(s)')

xlabel('m')

grid on;

figure(21)

plot(mtimeAll,ARTime,'bo-','MarkerFaceColor','b') title('重启的Arnoldi算法求解运行时间与m的关系曲线') ylabel('求解运行时间(s)')

xlabel('m')

grid on;

李庆扬数值分析第五版习题复习资料清华大学出版社

第一章 绪论 1.设0x >,x 的相对误差为δ,求ln x 的误差。 解:近似值* x 的相对误差为* **** r e x x e x x δ-= = = 而ln x 的误差为()1 ln *ln *ln ** e x x x e x =-≈ 进而有(ln *)x εδ≈ 2.设x 的相对误差为2%,求n x 的相对误差。 解:设()n f x x =,则函数的条件数为'() | |() p xf x C f x = 又1 '()n f x nx -=Q , 1 ||n p x nx C n n -?∴== 又((*))(*)r p r x n C x εε≈?Q 且(*)r e x 为2 ((*))0.02n r x n ε∴≈ 3.下列各数都是经过四舍五入得到的近似数,即误差限不超过最后一位的半个单位,试指 出它们是几位有效数字:*1 1.1021x =,*20.031x =, *3385.6x =, *456.430x =,* 57 1.0.x =? 解:* 1 1.1021x =是五位有效数字; *20.031x =是二位有效数字; *3385.6x =是四位有效数字; *456.430x =是五位有效数字; *57 1.0.x =?是二位有效数字。 4.利用公式(2.3)求下列各近似值的误差限:(1) ***124x x x ++,(2) ***123x x x ,(3) **24/x x . 其中**** 1234,,,x x x x 均为第3题所给的数。 解:

*4 1* 3 2* 13* 3 4* 1 51()1021()1021()1021()1021()102 x x x x x εεεεε-----=?=?=?=?=? *** 124***1244333 (1)()()()() 1111010102221.0510x x x x x x εεεε----++=++=?+?+?=? *** 123*********123231132143 (2)() ()()() 111 1.10210.031100.031385.610 1.1021385.610222 0.215 x x x x x x x x x x x x εεεε---=++=???+???+???≈ ** 24**** 24422 *4 33 5 (3)(/) ()() 11 0.0311056.430102256.43056.430 10x x x x x x x εεε---+≈ ??+??= ?= 5计算球体积要使相对误差限为1,问度量半径R 时允许的相对误差限是多少? 解:球体体积为343 V R π= 则何种函数的条件数为 2 3'4343 p R V R R C V R ππ===g g (*)(*)3(*)r p r r V C R R εεε∴≈=g 又(*)1r V ε=Q

清华大学数值分析A第一次作业

7、设y0=28,按递推公式 y n=y n?1? 1 100 783,n=1,2,… 计算y100,若取≈27.982,试问计算y100将有多大误差? 答:y100=y99?1 100783=y98?2 100 783=?=y0?100 100 783=28?783 若取783≈27.982,则y100≈28?27.982=0.018,只有2位有效数字,y100的最大误差位0.001 10、设f x=ln?(x? x2?1),它等价于f x=?ln?(x+ x2?1)。分别计算f30,开方和对数取6位有效数字。试问哪一个公式计算结果可靠?为什么? 答: x2?1≈29.9833 则对于f x=ln x?2?1,f30≈?4.09235 对于f x=?ln x+2?1,f30≈?4.09407 而f30= ln?(30?2?1) ,约为?4.09407,则f x=?ln?(x+ x2?1)计算结果更可靠。这是因为在公式f x=ln?(x? x2?1)中,存在两相近数相减(x? x2?1)的情况,导致算法数值不稳定。 11、求方程x2+62x+1=0的两个根,使它们具有四位有效数字。 答:x12=?62±622?4 2 =?31±312?1 则 x1=?31?312?1≈?31?30.98=?61.98 x2=?31+312?1= 1 31+312?1 ≈? 1 ≈?0.01613

12.(1)、计算101.1?101,要求具有4位有效数字 答:101.1?101= 101.1+101≈0.1 10.05+10.05 ≈0.004975 14、试导出计算积分I n=x n 4x+1dx 1 的一个递推公式,并讨论所得公式是否计算稳定。 答:I n=x n 4x+1dx 1 0= 1 4 4x+1x n?1?1 4 x n?1 4x+1 dx= 1 1 4 x n?1 1 dx?1 4 x n?1 4x+1 dx 1 = 1 4n ? 1 4 I n?1,n=1,2… I0= 1 dx= ln5 1 记εn为I n的误差,则由递推公式可得 εn=?1 εn?1=?=(? 1 )nε0 当n增大时,εn是减小的,故递推公式是计算稳定的。

清华大学高等数值计算(李津)实践题目一(共轭梯度CG法,Lanczos算法与MINRES算法)

高等数值计算实践题目一 1. 实践目的 本次计算实践主要是在掌握共轭梯度法,Lanczos 算法与MINRES 算法的基础上,进一步探讨这3种算法的数值性质,主要研究特征值特征向量对算法收敛性的影响。 2. 实践过程 (一)生成矩阵 (1)作5个100阶对角阵i D 如下: 1D 对角元:1,1,...,20,1+0.1(-20),21,...,100j j d j d j j ==== 2D 对角元:1,1,...,20,1+(-20),21,...,100j j d j d j j ==== 3D 对角元:,1,...,80,81,81,...,100j j d j j d j ==== 4D 对角元:,1,...,40,41,41,...,60,41+(60),61,...,100j j j d j j d j d j j =====-= 5D 对角元:,1,...,100j d j j == 记i D 的最大模特征值和最小模特征值分别为1i λ和i n λ,则i D 特征值分布有如下特点: 1D 的特征值有较多接近于i n λ,并且1/i i n λλ较小, 2D 的特征值有较多接近于i n λ,并且1/i i n λλ较大, 3D 的特征值有较多接近于1i λ,并且1/i i n λλ较大, 4D 的特征值有较多接近于中间模特征值,并且1/i i n λλ较大, 5D 的特征值均匀分布,并且1/i i n λλ较大 (2)随机生成10个100阶矩阵j M : (100(100))j M fix rand = 并作它们的QR 分解,得j Q 和j R ,这样可得50个对称的矩阵T ij j i j A Q DQ =,其中i D 的对角元就是ij A 的特征值,若它们都大于0,则ij A 正定,j Q 的列就是相应的特征向量。结合(1)可知,ij A 都是对称正定阵。

数值分析实验报告176453

实验报告 插值法 数学实验室 数值逼近 算法设计 级 ____________________________ 号 ____________________________ 名 _____________________________ 实验项目名称 实验室 所属课程名称 实验类型 实验日期

实验概述: 【实验目的及要求】 本次实验的目的是熟练《数值分析》第二章“插值法”的相关内容,掌握三种插 多项式插值,三次样条插值,拉格朗日插值,并比较三种插值方法的 优劣。 本次试验要求编写牛顿多项式插值,三次样条插值,拉格朗日插值的程序编码,并 去实现。 【实验原理】 《数值分析》第二章“插值法”的相关内容,包括:牛顿多项式插值,三次样条插值, 拉格朗日 插值的相应算法和相关性质。 【实验环境】(使用的软硬件) 软件: MATLAB 2012a 硬件: 电脑型号:联想 Lenovo 昭阳E46A 笔记本电脑 操作系统: Win dows 8专业版 处理器:In tel ( R Core ( TM i3 CPU M 350 @2.27GHz 2.27GHz 实验内容: 【实验方案设计】 第一步,将书上关于三种插值方法的内容转化成程序语言,用 MATLA B 现; 第二步,分别用牛顿多项式插值,三次样条插值,拉格朗日插值求解不同的问题。 【实验过程】(实验步骤、记录、数据、分析) 实验的主要步骤是:首先分析问题,根据分析设计 MATLA 程序,利用程序算出 问题答案,分析所得答案结果,再得出最后结论。 实验一: 已知函数在下列各点的值为 试用4次牛顿插值多项式 P 4( x )及三次样条函数 S ( x )(自然边界条件)对数据进行插值。 用图给出{( X i , y i ), X i =0.2+0.08i , i=0 , 1, 11, 10 } , P 4 ( x )及 S ( x )。 值方法:牛顿 在MATLAB 件中

数值分析实验报告模板

数值分析实验报告模板 篇一:数值分析实验报告(一)(完整) 数值分析实验报告 1 2 3 4 5 篇二:数值分析实验报告 实验报告一 题目:非线性方程求解 摘要:非线性方程的解析解通常很难给出,因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。利用二分法求解给定非线性方程的根,在给定的范围内,假设f(x,y)在[a,b]上连续,f(a)xf(b) 直接影响迭代的次数甚至迭代的收敛与发散。即若x0 偏离所求根较远,Newton法可能发散的结论。并且本实验中还利用利用改进的Newton法求解同样的方程,且将结果与Newton法的结果比较分析。 前言:(目的和意义) 掌握二分法与Newton法的基本原理和应用。掌握二分法的原理,验证二分法,在选对有根区间的前提下,必是收

敛,但精度不够。熟悉Matlab语言编程,学习编程要点。体会Newton使用时的优点,和局部收敛性,而在初值选取不当时,会发散。 数学原理: 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法:二分法和Newton法。 对于二分法,其数学实质就是说对于给定的待求解的方程f(x),其在[a,b]上连续,f(a)f(b) Newton法通常预先要给出一个猜测初值x0,然后根据其迭代公式xk?1?xk?f(xk) f'(xk) 产生逼近解x*的迭代数列{xk},这就是Newton法的思想。当x0接近x*时收敛很快,但是当x0选择不好时,可能会发散,因此初值的选取很重要。另外,若将该迭代公式改进为 xk?1?xk?rf(xk) 'f(xk) 其中r为要求的方程的根的重数,这就是改进的Newton 法,当求解已知重数的方程的根时,在同种条件下其收敛速度要比Newton法快的多。 程序设计: 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式,如下 function y=f(x);

《MATLAB与数值分析》第一次上机实验报告

电子科技大学电子工程学院标准实验报告(实验)课程名称MATLAB与数值分析 学生姓名:李培睿 学号:2013020904026 指导教师:程建

一、实验名称 《MATLAB与数值分析》第一次上机实验 二、实验目的 1. 熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算 操作。(用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序) 2. 熟练掌握算术符号操作和基本运算操作,包括矩阵合并、向量合并、符号 转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。(用.m文件编写进行符号因式分解和函数求反的程序) 3. 掌握Matlab函数的编写规范。 4、掌握Matlab常用的绘图处理操作,包括:基本平面图、图形注释命令、 三维曲线和面的填充、三维等高线等。(用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形,并给出充分的图形注释) 5. 熟练操作MATLAB软件平台,能利用M文件完成MATLAB的程序设计。 三、实验内容 1. 编程实现以下数列的图像,用户能输入不同的初始值以及系数。并以x, y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图,允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。不允许使用sort 函数。 四、实验数据及结果分析 题目一: ①在Editor窗口编写函数代码如下:

数值分析实验报告

实验一 误差分析 实验(病态问题) 实验目的:算法有“优”与“劣”之分,问题也有“好”与“坏”之别。对数值方法的研究而言,所谓坏问题就是问题本身对扰动敏感者,反之属于好问题。通过本实验可获得一个初步体会。 数值分析的大部分研究课题中,如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决,当然一般要付出一些代价(如耗用更多的机器时间、占用更多的存储空间等)。 问题提出:考虑一个高次的代数多项式 )1.1() ()20()2)(1()(20 1∏=-=---=k k x x x x x p Λ 显然该多项式的全部根为1,2,…,20共计20个,且每个根都是单重的。现考虑该多项式的一个扰动 )2.1(0 )(19=+x x p ε 其中ε是一个非常小的数。这相当于是对()中19x 的系数作一个小的扰动。我们希望比较()和()根的差别,从而分析方程()的解对扰动的敏感性。 实验内容:为了实现方便,我们先介绍两个Matlab 函数:“roots ”和“poly ”。 roots(a)u = 其中若变量a 存储n+1维的向量,则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,,,+n a a a Λ,则输出u 的各分量是多项式方程 01121=+++++-n n n n a x a x a x a Λ 的全部根;而函数 poly(v)b = 的输出b 是一个n+1维变量,它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“poly ”是两个互逆的运算函数。 ;000000001.0=ess );21,1(zeros ve = ;)2(ess ve = ))20:1((ve poly roots +

清华大学贾仲孝老师高等数值分析报告第二次实验

高等数值分析第二次实验作业

T1.构造例子特征值全部在右半平面时, 观察基本的Arnoldi 方法和GMRES 方法的数值性态, 和相应重新启动算法的收敛性. Answer: (1) 构造特征值均在右半平面的矩阵A : 根据实Schur 分解,构造对角矩阵D 由n 个块形成,每个对角块具有如下形式,对应一对特 征值i i i αβ± i i i i i S αββα-?? = ??? 这样D=diag(S 1,S 2,S 3……S n )矩阵的特征值均分布在右半平面。生成矩阵A=U T AU ,其中U 为 正交阵,则A 矩阵的特征值也均在右半平面。不妨构造A 如下所示: 2211112222 /2/2/2/2N N A n n n n ?-?? ? ? ?- ? = ? ? ? - ? ?? ? 由于选择初值与右端项:x0=zeros(2*N,1);b=ones(2*N,1); 则生成矩阵A 的过程代码如下所示: N=500 %生成A 为2N 阶 A=zeros(2*N); for a=1:N A(2*a-1,2*a-1)=a; A(2*a-1,2*a)=-a; A(2*a,2*a-1)=a; A(2*a,2*a)=a; end U = orth(rand(2*N,2*N)); A1 = U'*A*U; (2) 观察基本的Arnoldi 和GMRES 方法 编写基本的Arnoldi 函数与基本GMRES 函数,具体代码见附录。 function [x,rm,flag]=Arnoldi(A,b,x0,tol,m) function [x,rm,flag]=GMRES(A,b,x0,tol,m) 输入:A 为方程组系数矩阵,b 为右端项,x0为初值,tol 为停机准则,m 为人为限制的最大步数。 输出:x 为方程的解,rm 为残差向量,flag 为解是否收敛的标志。 外程序如下所示: e=1e-6; m=700;

数值分析实验报告

数值分析实验报告 姓名:周茹 学号: 912113850115 专业:数学与应用数学 指导老师:李建良

线性方程组的数值实验 一、课题名字:求解双对角线性方程组 二、问题描述 考虑一种特殊的对角线元素不为零的双对角线性方程组(以n=7为例) ?????????? ?????? ? ???? ?d a d a d a d a d a d a d 766 55 44 3 32 211??????????????????????x x x x x x x 7654321=?????????? ? ???????????b b b b b b b 7654321 写出一般的n (奇数)阶方程组程序(不要用消元法,因为不用它可以十分方便的解出这个方程组) 。 三、摘要 本文提出解三对角矩阵的一种十分简便的方法——追赶法,该算法适用于任意三对角方程组的求解。 四、引言 对于一般给定的d Ax =,我们可以用高斯消去法求解。但是高斯消去法过程复杂繁琐。对于特殊的三对角矩阵,如果A 是不可约的弱对角占优矩阵,可以将A 分解为UL ,再运用追赶法求解。

五、计算公式(数学模型) 对于形如????? ?? ????? ??? ?---b a c b a c b a c b n n n n n 111 2 2 2 11... ... ...的三对角矩阵UL A =,容易验证U 、L 具有如下形式: ??????? ????? ??? ?=u a u a u a u n n U ...... 3 3 22 1 , ?? ????? ? ?? ??????=1 (1) 1132 1l l l L 比较UL A =两边元素,可以得到 ? ?? ??-== = l a b u u c l b u i i i i i i 111 i=2, 3, ... ,n 考虑三对角线系数矩阵的线性方程组 f Ax = 这里()T n x x x x ... 2 1 = ,()T n f f f f ... 2 1 = 令y Lx =,则有 f Uy = 于是有 ()?????-== --u y a f y u f y i i i i i 1 1 11 1 * i=2, 3, ... ,n 再根据y Lx =可得到

数值分析实验报告

实验一、误差分析 一、实验目的 1.通过上机编程,复习巩固以前所学程序设计语言及上机操作指令; 2.通过上机计算,了解误差、绝对误差、误差界、相对误差界的有关概念; 3.通过上机计算,了解舍入误差所引起的数值不稳定性。 二.实验原理 误差问题是数值分析的基础,又是数值分析中一个困难的课题。在实际计算中,如果选用了不同的算法,由于舍入误差的影响,将会得到截然不同的结果。因此,选取算法时注重分析舍入误差的影响,在实际计算中是十分重要的。同时,由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差,因此算法的好坏会影响到数值结果的精度。 三.实验内容 对20,,2,1,0 =n ,计算定积分 ?+=10 5dx x x y n n . 算法1:利用递推公式 151--=n n y n y , 20,,2,1 =n , 取 ?≈-=+=1 00182322.05ln 6ln 51dx x y . 算法2:利用递推公式 n n y n y 51511-= - 1,,19,20 =n . 注意到 ???=≤+≤=10 10202010201051515611261dx x dx x x dx x , 取 008730.0)12611051(20120≈+≈y .: 四.实验程序及运行结果 程序一: t=log(6)-log(5);

n=1; y(1)=t; for k=2:1:20 y(k)=1/k-5*y(k-1); n=n+1; end y y =0.0884 y =0.0581 y =0.0431 y =0.0346 y =0.0271 y =0.0313 y =-0.0134 y =0.1920 y =-0.8487 y =4.3436 y =-21.6268 y =108.2176 y =-541.0110 y =2.7051e+003 y =-1.3526e+004 y =6.7628e+004 y =-3.3814e+005 y =1.6907e+006 y =-8.4535e+006 y =4.2267e+007 程序2: y=zeros(20,1); n=1; y1=(1/105+1/126)/2;y(20)=y1; for k=20:-1:2 y(k-1)=1/(5*k)-(1/5)*y(k); n=n+1; end 运行结果:y = 0.0884 0.0580 0.0431 0.0343 0.0285 0.0212 0.0188 0.0169

数值分析实验2014

数值分析实验(2014,9,16~10,28) 信计1201班,人数34人 数学系机房 数值分析 计算实习报告册 专业 学号 姓名 2014~2015年第一学期

实验一 数值计算的工具 Matlab 1.解释下MATLAB 程序的输出结果 程序: t=0.1 n=1:10 e=n/10-n*t e 的结果:0 0 -5.5511e-017 0 0 -1.1102e-016 -1.1102e-016 0 0 0 2.下面MATLAB 程序的的功能是什么? 程序: x=1;while 1+x>1,x=x/2,pause(0.02),end 用迭代法求出x=x/2,的最小值 x=1;while x+x>x,x=2*x,pause(0.02),end 用迭代法求出x=2*x,的值,使得2x>X x=1;while x+x>x,x=x/2,pause(0.02),end 用迭代法求出x=x/2,的最小值,使得2x>X 3.考虑下面二次代数方程的求解问题 02=++c bx ax 公式a ac b b x 242-+-=是熟知的,与之等价地有a c b b c x 422-+-=,对于 1,100000000,1===c b a ,应当如何选择算法。 应该用a ac b b x 242-+-=计算,因为b 做分母 4.函数)sin(x 有幂级数展开...! 7!5!3sin 7 53+-+-=x x x x x 利用幂级数计算x sin 的MATLAB 程序为

function s=powersin(x) s=0; t=x; n=1; while s+t~=s; s=s+t ; t=-x^2/((n+1)*(n+2))*t ; n=n+2; end t1=cputime; pause(10); t2=cputime; t0=t2-t1 (a)解释上述程序的终止准则。 当s+t=s ,终止循环。 (b)对于2/21,2/11 ,2/πππ=x 计算的进度是多少?分别计算多少项? X=pi/2时,s =1.0000 x=11pi/2时,s=-1.0000 x=21pi/2时,s =0.9999 Cputime 分别是0.1563 0.0469 0.0156 5.考虑调和级数∑∞ =11 n n ,它是微积分中的发散级数,在计算机上计算该级数的部 分和,会得到怎么样的结果,为什么? function s=fun(n) s=0; t=1/n; for i=1:n s=s+1/i; end 当n=100时s =5.1874 当n=80时s =4.9655 当n=50时,s =4.4992 当n=10时,s =2.9290 6.指数函数的级数展开...! 3!213 2++++=x x x e x ,如果对于0

数值分析实验报告总结

数值分析实验报告总结 随着电子计算机的普及与发展,科学计算已成为现代科 学的重要组成部分,因而数值计算方法的内容也愈来愈广泛和丰富。通过本学期的学习,主要掌握了一些数值方法的基本原理、具体算法,并通过编程在计算机上来实现这些算法。 算法算法是指由基本算术运算及运算顺序的规定构成的完 整的解题步骤。算法可以使用框图、算法语言、数学语言、自然语言来进行描述。具有的特征:正确性、有穷性、适用范围广、运算工作量少、使用资源少、逻辑结构简单、便于实现、计算结果可靠。 误差 计算机的计算结果通常是近似的,因此算法必有误差, 并且应能估计误差。误差是指近似值与真正值之差。绝对误差是指近似值与真正值之差或差的绝对值;相对误差:是指近似值与真正值之比或比的绝对值。误差来源见表 第三章泛函分析泛函分析概要 泛函分析是研究“函数的函数”、函数空间和它们之间 变换的一门较新的数学分支,隶属分析数学。它以各种学科

如果 a 是相容范数,且任何满足 为具体背景,在集合的基础上,把客观世界中的研究对象抽 范数 范数,是具有“长度”概念的函数。在线性代数、泛函 分析及相关的数学领域,泛函是一个函数,其为矢量空间内 的所有矢量赋予非零的正长度或大小。这里以 Cn 空间为例, Rn 空间类似。最常用的范数就是 P-范数。那么 当P 取1, 2 ,s 的时候分别是以下几种最简单的情形: 其中2-范数就是通常意义下的距离。 对于这些范数有以下不等式: 1 < n1/2 另外,若p 和q 是赫德尔共轭指标,即 1/p+1/q=1 么有赫德尔不等式: II = ||xH*y| 当p=q=2时就是柯西-许瓦兹不等式 般来讲矩阵范数除了正定性,齐次性和三角不等式之 矩阵范数通常也称为相容范数。 象为元素和空间。女口:距离空间,赋范线性空间, 内积空间。 1-范数: 1= x1 + x2 +?+ xn 2-范数: x 2=1/2 8 -范数: 8 =max oo ,那 外,还规定其必须满足相容性: 所以

清华大学杨顶辉数值分析第6次作业

清华大学杨顶辉数值分析第6次作业

9.令*()(21),[0,1]n n T x T x x =-∈,试证*{()}n T x 是在[0,1]上带权 2 ()x x x ρ= -****0123(),(),(),()T x T x T x T x . 证明: 1 1 **2 1 1 * *20 12 2 1**20 ()()()(21)(21)211()()()()()211()22 ()()1()1()()()()()1n m n m n m n m n m n n m n m x T x T x dx x T x dx x x t x x T x T x dx t T t dt t t t T t dt t T x x x T x T x dx t T t t ρρρ---=---=-=++-= --= -???? ?令,则 由切比雪夫多项式1 01=02 m n dt m n m n ππ ≠??? =≠??==??? 所以*{()}n T x 是在[0,1]上带权2 ()x x x ρ= - *00*11* 2 2 2 2*33233()(21)1()(21)21 ()(21)2(21)188()(21)4(21)3(21)3248181 T x T x T x T x x T x T x x x x T x T x x x x x x =-==-=-=-=--=-=-=---=-+- 14.已知实验数据如下: i x 19 25 31 38 44 i y 19.0 32.3 49.0 73.3 97.8 用最小二乘法求形如2y a bx =+的经验公式,并求均方误差 解: 法方程为

数值分析实验题目及解答

内容包括: 实验题目1:算法的数值稳定性实验 实验题目2:LU分解实验 实验题目3:三次样条插值外推样条实验 实验题目4:第二类Fredholm 积分方程实验实验题目5:M级显式R_K法

实验题目:算法的数值稳定性实验 实验内容:计算积分()1 0()d 1515n x I n x a x ==+? (n=1,2,…,20) 易得到下面递推公式 ()()1 1I n aI n n =--+ 并有估计式 ()() ()() 1 1 111I n a n a n << +++ 计算方法: 算法一:采用下面递推公式计算: ()()1 1I n aI n n =--+ ()1,2,,20 n = 取初值()116 0ln ln 15a I a +== 算法二: 采用下面递推公式计算: ()()111I n I n a n ??-= -+???? ()20,19,,1 n =

结果分析:(分析哪个好哪个不好,原因是什么) 我觉得算法二比较好, 原因一:根据式 ()() ()() 1 1 111I n a n a n << +++得知,I(n)不可能小于 零,而算法一的计算结果有部分结果小于零。原因二:对算法一记初始误差 ε0=/I 0-I(0)/>0; 则εn =/I n -I(n)/=a/I n-1-I(n-1)/=a n *ε0 由此可知,当n=20时, ε20把ε0放大了a 20倍,其结果造成严重的。 而对于算法二^ ^ 11n n a εε-= ,…, ^ ^ 01 n n a εε=,尽管有初始误差^ 20ε,但随着计算的进程,这个误差的影响不断减小。 附:源程序:(把源程序附上) 算法一程序: >> format long >> a=15;I=log(16/15); for n=1:20 n I=-a*I+1/n end 算法二程序: >> format long >> a=15;I=31/10080; >> for n=20:-1:1 n I I=1/a*(-I+1/n); End

数值分析实验报告

学生实验报告实验课程名称 开课实验室 学院年级专业班 学生姓名学号 开课时间至学年学期

if(A(m,k)~=0) if(m~=k) A([k m],:)=A([m k],:); %换行 end A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k))*A(k, k:c); %消去end end x=zeros(length(b),1); %回代求解 x(n)=A(n,c)/A(n,n); for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n))/A(k,k); end y=x; format short;%设置为默认格式显示,显示5位 (2)建立MATLAB界面 利用MA TLAB的GUI建立如下界面求解线性方程组: 详见程序。 五、计算实例、数据、结果、分析 下面我们对以上的结果进行测试,求解:

? ? ? ? ? ? ? ? ? ? ? ? - = ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? - - - - - - 7 2 5 10 13 9 14 4 4 3 2 1 13 12 4 3 3 10 2 4 3 2 1 x x x x 输入数据后点击和,得到如下结果: 更改以上数据进行测试,求解如下方程组: 1 2 3 4 43211 34321 23431 12341 x x x x ?? ???? ?? ???? ?? ???? = ?? ???? - ?? ???? - ???? ?? 得到如下结果:

数值分析实验报告资料

机电工程学院 机械工程 陈星星 6720150109 《数值分析》课程设计实验报告 实验一 函数插值方法 一、问题提出 对于给定的一元函数)(x f y =的n+1个节点值(),0,1,,j j y f x j n ==。试用Lagrange 公式求其插值多项式或分段二次Lagrange 插值多项式。 数据如下: (1 求五次Lagrange 多项式5L ()x ,计算(0.596)f ,(0.99)f 的值。(提示:结果为(0.596)0.625732f ≈, (0.99) 1.05423f ≈) 实验步骤: 第一步:先在matlab 中定义lagran 的M 文件为拉格朗日函数 代码为: function[c,l]=lagran(x,y) w=length(x); n=w-1; l=zeros(w,w); for k=1:n+1 v=1; for j=1:n+1 if(k~=j) v=conv(v,poly(x(j)))/(x(k)-x(j)); end end l(k,:)=v; end c=y*l; end

第二步:然后在matlab命令窗口输入: >>>> x=[0.4 0.55 0.65 0.80,0.95 1.05];y=[0.41075 0.57815 0.69675 0.90 1.00 1.25382]; >>p = lagran(x,y) 回车得到: P = 121.6264 -422.7503 572.5667 -377.2549 121.9718 -15.0845 由此得出所求拉格朗日多项式为 p(x)=121.6264x5-422.7503x4+572.5667x3-377.2549x2+121.9718x-15.0845 第三步:在编辑窗口输入如下命令: >> x=[0.4 0.55 0.65 0.80,0.95 1.05]; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718 *x-15.0845; >> plot(x,y) 命令执行后得到如下图所示图形,然后 >> x=0.596; >> y=121.6264*x.^5-422.7503*x.^4+572.5667*x.^3-377.2549*x.^2+121.9718 *x-15.084 y =0.6257 得到f(0.596)=0.6257 同理得到f(0.99)=1.0542

清华大学杨顶辉数值分析第6次作业

9.令*()(21),[0,1]n n T x T x x =-∈,试证*{()}n T x 是在[0,1] 上带权()x ρ=的正交多项式,并求****0123(),(),(),()T x T x T x T x . 证明: 1 1 * *0 1 1 * *011**0 ()()()(21)(21)211()()()()()2()()()()()()()()n m n m n m n m n m n n m n m x T x T x dx x T x dx t x x T x T x dx t T t dt t T t dt T x x T x T x dx t T t ρρρ---=--=-== = ???? ?令,则 由切比雪夫多项式1 01=02 m n dt m n m n ππ ≠??? =≠??==??? 所以*{()}n T x 是在[0,1] 上带权()x ρ= *00*11* 22 2 2*33233()(21)1()(21)21 ()(21)2(21)188()(21)4(21)3(21)3248181 T x T x T x T x x T x T x x x x T x T x x x x x x =-==-=-=-=--=-=-=---=-+- 14.已知实验数据如下: 用最小二乘法求形如2y a bx =+的经验公式,并求均方误差 解: 法方程为

22222(1,)(1,1)(1,)(,)(,1)(,)a y x b x y x x x ?????? =???? ?????? ?? 即 5 5327271.453277277699369321.5a b ??????=???????????? 解得 0.972579 0.050035a b =?? =? 拟合公式为20.9725790.050035y x =+ 均方误差 2 4 2 2 0[]0.015023i i i y a bx σ==--=∑ 21.给出()ln f x x =的函数表如下: 用拉格朗日插值求ln 0.54的近似值并估计误差(计算取1n =及2n =) 解:1n =时,取010.5,0.6x x == 由拉格朗日插值定理有 1 100.60.5 0.693147 0.510826 0.50.(60.60.51.82321)0 1.()6047()52 j j j x x x L x f x l x ==------=-=∑ 所以1ln0.54(0.54)0.620219L ≈=- 误差为ln 0.54(0.620219)= 0.004032ε=-- 2n =时,取0120.4,0.5,0.6x x x === 由拉格朗日插值定理有

清华大学高等数值计算(李津)实践题目二(SVD计算及图像压缩)(包含matlab代码)

第1部分 方法介绍 奇异值分解(SVD )定理: 设m n A R ?∈,则存在正交矩阵m m V R ?∈和n n U R ?∈,使得 T O A V U O O ∑??=?? ?? 其中12(,, ,)r diag σσσ∑=,而且120r σσσ≥≥≥>,(1,2, ,)i i r σ=称为A 的 奇异值,V 的第i 列称为A 的左奇异向量,U 的第i 列称为A 的右奇异向量。 注:不失一般性,可以假设m n ≥,(对于m n <的情况,可以先对A 转置,然后进行SVD 分解,最后对所得的SVD 分解式进行转置,就可以得到原来的SVD 分解式) 方法1:传统的SVD 算法 主要思想: 设()m n A R m n ?∈≥,先将A 二对角化,即构造正交矩阵1U 和1V 使得 110T B n U AV m n ?? =?? -?? 其中1200n n B δγγδ??? ???=?????? 然后,对三角矩阵T T B B =进行带Wilkinson 位移的对称QR 迭代得到:T B P BQ =。 当某个0i γ=时,B 具有形状12B O B O B ?? =? ??? ,此时可以将B 的奇异值问题分解为两个低阶二对角阵的奇异值分解问题;而当某个0i δ=时,可以适当选取'Given s 变换,使得第i 行元素全为零的二对角阵,因此,此时也可以将B 约化为两个低 阶二对角阵的奇异值分解问题。 在实际计算时,当i B δε∞≤或者() 1j j j γεδδ-≤+(这里ε是一个略大于机器精度的正数)时,就将i δ或者i γ视作零,就可以将B 分解为两个低阶二对角阵的奇异值分解问题。

MATLAB与数值分析实验报告一

MATLAB与数值分析实验报告 报告人:秦旸照 学号: 2015020901033 时间: 2016.4.8 电子科技大学电子工程学院

一、实验目的 实验一:MATLAB软件平台与程序设计实验 二、实验原理 1.熟练掌握矩阵的生成、加、减、乘、除、转置、行列式、逆、范数等运算操作。(用.m文件和Matlab函数编写一个对给定矩阵进行运算操作的程序) 2. 熟练掌握算术符号操作和基本运算操作,包括矩阵合并、向量合并、符号转换、展开符号表达式、符号因式分解、符号表达式的化简、代数方程的符号解析解、特征多项式、函数的反函数、函数计算器、微积分、常微分方程的符号解、符号函数的画图等。(用.m文件编写进行符号因式分解和函数求反的程序) 3. 掌握Matlab函数的编写规范。 4.掌握Matlab常用的绘图处理操作,包括:基本平面图、图形注释命令、三维曲线和面的填充、三维等高线等。(用.m文件编写在一个图形窗口上绘制正弦和余弦函数的图形,并给出充分的图形注释) 5. 熟练操作MATLAB软件平台,能利用M文件完成MATLAB的程序设计。 三、实验方案 1. 编程实现以下数列的图像,用户能输入不同的初始值以及系数。并以 x,y为坐标显示图像 x(n+1) = a*x(n)-b*(y(n)-x(n)^2); y(n+1) = b*x(n)+a*(y(n)-x(n)^2) 2. 编程实现奥运5环图,允许用户输入环的直径。 3. 实现对输入任意长度向量元素的冒泡排序的升序排列。 不允许使用sort函数。 四、实验结果 1. 编程实现以下数列的图像,用户能输入不同的初始值以及系数。并以 x,y为坐标显示图像

清华大学高等数值分析作业李津1——矩阵基础

20130917题目 求证:在矩阵的LU 分解中,1 11n n T n ij i j j i j L I e e α-==+??=- ??? ∑∑ 证明: 在高斯消去过程中,假设0jj a ≠ ,若a=0,可以通过列变换使得前面的条件成立,这里不考虑这种情况。 对矩阵A 进行LU 分解,()() () ()()1 11 1111L M n M M M n ---=-=??-………… , 其中()1n T n ij i j i j M j I e e α=+??=+ ??? ∑ ,i e 、j e 为n 维线性空间的自然基。 ()M j 是通过对单位阵进行初等变换得到, 通过逆向的变换则可以得到单位阵,由此很容易得到()M j 的逆矩阵为1n T n ij i j i j I e e α=+??- ???∑。故111n n T n ij i j n j i j L I e e I α-==+?? ??=- ? ? ????? ∏∑ 上式中的每一项均是初等变换,从右向左乘,则每乘一次相当于对右边的矩阵进行一次 向下乘法叠加的初等变换。由于最初的矩阵为单位阵,变换从右向左展开,因而每一次变换不改变已经更新的数据,既该变换是从右向左一列一列更新数据,故 11n n T n ij i j j i j L I e e α==+??=- ??? ∑∑。 数学证明:1n T ij i j i j e e α=+?? ???∑具有 ,0 00n j j A -?? ??? 和1,1000n j n j B -+-+?? ?? ? 的形式,且有 +1,-11,10000=000n j j n j n j A B --+-+???? ?????? ? 而1 1n n T ij i j j k i j e e α-==+?? ??? ∑∑具有1,1000n k n k B -+-+?? ???的形式,因此: 1 311111211121==n n n n n n T T T n ij i j n ij i j n ik i k j i j j i j k n i k n n T n i i n ik i i i k L I e e I e e I e e I e e I e ααααα---==+==+=-=+==+??????????????=---?? ? ? ? ? ? ? ? ???????????????????????=-- ? ? ?????∏∑∏∑∑∑∑∑……11211n n n T T k n ik i k k k i k e I e e α--===+????=- ?? ?????? ∑∑∑#

相关主题
相关文档 最新文档