上机作业1.docx

资源描述

上机作业1.docx

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

上机作业1.docx

上机作业1

2015级研究生

《计算方法》作业

2015年11月

上机作业1

数值试验3

3-1试验目的：

考察不动点迭代法的局部收敛性

试验内容：

至少采用3种迭代法，迭代100次，考察收敛性，改变初值符号，再做迭代。

分析收敛与发散的原因。

（1）迭代原理：

若实数

满足

，

称为函数

的一个不动点，迭代

称为不动点迭代，

称为迭代函数。

由不动点方程建立迭代法

，其中

称为初值，需要预先给定。

方程

分别对应下列不同形式的不动点方程：

取

，按

迭代，并分析收敛性。

（2）不动点迭代法代码

编制函数文件

Iteratepro.m

functiony=iteratepro（x）

x1=g（x）;

n=1;

while（norm（x1-x）>=1.0e-5）&（n<=100）

x=x1;

x1=g（x）;n=n+1;

end

编制函数文件

不动点方程

（1）

functiony=g（x）

y=（exp（x）-3）/2;

不动点方程

（2）

functiony=g（x）

y=exp（x）-x-3;

不动点方程（3）

functiony=g（x）

y=x-（2*x-exp（x）+3）/（2-exp（x））;

（3）迭代结果汇总表：

近似解

初值x=-0.5

迭代次数n

初值x=0.5

迭代次数n

格式

（1）

-1.1967

-1.3734

格式

（2）

-1.3734

格式（3）

-1.3734

（4）结果分析：

（1）三种迭代格式均收敛，近似解为

。

（2）当设定的初值接近真实解，所用的迭代次数较小。

（3）不同迭代格式，效果不一样，迭代格式（3）实质是newton迭代法，其迭代效果明显快于其他格式。

3-2试验目的：

考察Newton法求根的收敛速度

试验内容：

应用Newton法求解试验3-1中的方程，并与实验3-1中收敛的迭代法进行比较，考察收敛速度。

精确到0.00001。

（1）Newton迭代法原理:

在求解非线性方程

时，它的困难在于

是非线性函数，为克服这一困难，考虑它的线性展开。

设当前点为

，在

处的

展开式为

令

，解其方程得到

（2）程序代码

编制函数文件

newton.m

functiony=newton（x0）

x1=x0-fc（x0）/df（x0）;

n=1;

while（abs（x1-x0）>=1.0e-5）&（n<=100）

x0=x1;

x1=x0-fc（x0）/df（x0）;n=n+1;

end

对于试验3-1编制函数文件

fc.m

functiony=fc（x）

y=2*x-exp（x）+3;

df.m

functiony=df（x）

y=2-exp（x）;

（3）运行结果

结果汇总表：

近似解

初值x=-0.5

迭代次数n

初值x=0.5

迭代次数n

格式

（1）

-1.1967

-1.3734

格式

（2）

-1.3734

Newton迭代

-1.3734

结果分析：

（1）三者结果均收敛，解

。

（2）Newton迭代法迭代次数最小，收敛速度明显优于格式

（1），

（2）。

3-3试验目的：

掌握求重根的方法

试验内容：

分别用Newton法和不动点迭代法求解方程

，考察收敛速度，再用求重根的两种方法求方程的根，精确到0.00001。

（1）实验原理同3-1,3-2；

（2）程序代码

Newton迭代法程序

函数文件

fc.m

functiony=fc（x）

y=x-sin（x）;

df.m

functiony=df（x）

y=1-cos（x）;

程序代码

newton.m

functiony=newton（x0）

x1=x0-fc（x0）/df（x0）;

n=1;

whileabs（x1-x0）>=1.0e-5

x0=x1;

x1=x0-fc（x0）/df（x0）;n=n+1;

end

不动点迭代法

编制函数文件

Iteratepro.m

functiony=iteratepro（x）

x1=g（x）;

n=1;

whilenorm（x1-x）>=1.0e-6

x=x1;

x1=g（x）;n=n+1;

end

函数文件

g.m

functiony=g（x）

y=sin（x）;

运行结果

近似解

初值x=1

迭代次数n

牛顿法

1.9449e-04

不动点法

0.0182

k=410

（3）为了提高牛顿法求重根的收敛阶采用以下俩种方法：

方法一：

程序：

function[p1,err,y]=newtonroot（f,df,p0,eps,max1）

fork=1:

max1

p1=p0-2*（feval（'f',p0）/feval（'df',p0））；

err=abs（p1-p0）;

p0=p1;

p1,err,k,y=feval（'f',p1）

if（err

break;

end

取初值p=0.5,当迭代次数k=9满足精度要求,近似解为2.4925e-05

方法二：

程序：

function[p1,err,y]=newtonroot（f,df,ddf,p0,eps,max1）

fork=1:

max1

p1=p0-（feval（'f',p0）*feval（'df',p0））/（（feval（'df',p0））^2-feval（'f',p0）*feval（'ddf',p0））

err=abs（p1-p0）;

p0=p1;

p1,err,k,y=feval（'f',p1）

if（err

break;

end

取初值p=0.5,k=3执行结果为：

p1=1.5022e-08err=2.2741e-08

ans=1.5022e-08

以上四种方法中简单迭代法的收敛速度最慢最后一种收敛速度最快。

3-4试验目的：

体验Steffensen’smethod加速技巧

试验内容：

先用Newton法求解方程

，再用Steffensen’smethod求解，比较迭代步数，精确到0.00001。

（1）简单原理：

其基本思想是：

对给定的初值p0（0）,首先应用不动点迭代法pn+1=g（pn）计算俩步，然后用Atiken公式加速。

（2）程序代码

Newton法迭代

程序：

function[p1,err,y]=newtonroot（f,df,p0,eps,max1）

fork=1:

max1

p1=p0-feval（'f',p0）/feval（'df',p0）;

err=abs（p1-p0）;

p0=p1;

p1,err,k,y=feval（'f',p1）

if（err

break;

end

End

定义迭代函数：

存储为f.m

functiony=f（x）

y=x-tan（x）

end

定义迭代函数：

存储为df.m

functiong=df（x）

g=1-1/（cos（x））^2

end

Atiken公式加速程序：

function[n,err,p1,p2]=steffen（g,p0,tol,max1）

n=0,p

（1）=p0;

whilen<=max1

fork=2:

p（k）=feval（'g',p（k-1））;

end

p1=p

（1）-（p

（2）-p

（1））^2/（p（3）-2*p

（2）+p

（1））;

err=abs（p1-p

（1））;

n=n+1;

（1）=p1

n,err;

if（err

break

end

（3）结果与分析

牛顿法

取初值p0=0.5，迭代次数k=20满足精度要求.

p1=1.6018e-04err=8.0089e-05k=20

Ans=1.6018e-04

Atiken公式加速

取初值p0=0.5，迭代次数k=20满足精度要求。

err=8.6054e-05n=20p=0.1776×1.0e-03

二者迭代步数相同。

从结果上看牛顿法更精确一些。

上机作业2

数值试验5-1

实验目的：

熟悉Jacobi、Seidel、Sor迭代法，了解松弛因子对收敛速度的影响。

实验内容：

分别用Jacobi、Seidel、Sor

迭代法求解下面的方程组，并作结果分析。

初值

，精度要求：

（1）

（2）

（1）实验原理

Jacobi迭代法求解方程组

（1）

程序代码：

先建立M文件，保存文件名为jacobi.m。

functiony=jacobi（a,b,x0）

D=diag（diag（a））;

U=-triu（a,1）;

L=-tril（a,-1）;

B=D\（L+U）;

f=D\b;

y=B*x0+f;n=1;

whilenorm（y-x0）>=1.0e-5

x0=y;

y=B*x0+f;n=n+1;

end

Seidel迭代法

D=diag（diag（A））;

D0=inv（D）;

L=tril（A）-D;

U=triu（A）-D;

M=L+D;

M0=inv（M）;

B=-M0*U;

d=M0*b;

x0=[00000]';

x=B*x0+d;

whilenorm（x-x0）>1e-5*norm（x）

x0=x;

x=B*x0+d;

i=i+1;

end

Sor迭代式求解方程组

程序代码：

建立文件名为sor.m的M文件，具体如下：

functiony=sor（a,b,w,x0）

D=diag（diag（a））;

U=-triu（a,1）;

L=-tril（a,-1）;

B=（D-w*L）\（（1-w）*D+w*U）;

f=（D-w*L）\b*w;

y=B*x0+f;n=1;

whilenorm（y-x0）>=1.0e-5

x0=y;

y=B*x0+f;

n=n+1;

End

2程序计算结果

三种迭代法求解方程组

（1）结果列表

Sor

Jacobi

Seidel

0.8

1.1

1.2

1.3

1.4

1.5

—

0.3906

NaN

0.3906

0.1678

NaN

0.1678

0.0996

NaN

0.0996

0.1754

NaN

0.1754

0.0447

-Inf

0.0447

迭代次数

8357

三种迭代法求解方程组

（2）结果列表

Sor

Jacobi

Seidel

0.8

1.1

1.2

1.3

1.4

1.5

—

0.4720

NaN

0.4720

0.1330

NaN

0.1330

-0.1302

NaN

-0.1302

0.1628

NaN

0.1628

0.1701

NaN

0.1701

迭代次数

8058

2260

3结果分析：

（1）Seidel迭代法收敛速度最快，迭代次数最少。

（2）当松弛因子

越大，其收敛速度越慢，迭代次数越多

（3）当松弛因子选的合适时（

在1.0附近），Sor迭代法的收敛速度要比Jacobi迭代法快点。

上机作业3

数值试验5-2

实验目的：

掌握Newton法与最速下降法求解非线性方程组，观察各自的优势。

实验内容：

（1）分别用Newton法和最速下降法求解下面非线性方程组，初值

精度要求

（2）改变初值

，再用如上两种方法求解，得到什么结果？

（3）采用初值

，先用最速下降法求解3步，再用Newton迭代法求解，得到什么结果？

对以上运算结果作分析。

（1）简单原理：

Newton法:

当方程组的解

可表示为

称为求解非线性方程组的newton法迭代格式。

最速下降法简单原理：

在包含

的某个区域内任取一初始点

，然后在此点沿着

下降最快的负梯度方向，取一个合适的步长，求得下一个

，依次循环，直到

来终止计算，此时

就是满足要求的近似解。

（2）MATLAB程序

牛顿法解非线性方程组程序

symsx1x2x3

f1=3*x1-cos（x2*x3）-（1/2）;

f2=（x1）^2-81*（x2+0.1）^2+sin（x3）+1.06;

f3=exp（-x1*x2）+20*x3+（（10*pi-3）/3）;

g1=diff（f1,x1）;g2=diff（f1,x2）;g3=diff（f1,x3）;

g4=diff（f2,x1）;g5=diff（f2,x2）;g6=diff（f2,x3）;

g7=diff（f3,x1）;g8=diff（f3,x2）;g9=diff（f3,x3）;

x0=[0.1;0.1;-0.1];

Tol=10^（-5）;

%x0=[20;20;20];

F2=[subs（f1,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）;

subs（f2,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）;

subs（f3,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）];

F3=[subs（g1,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g2,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g3,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）;

subs（g4,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g5,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g6,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）;

subs（g7,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g8,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g9,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）];

k=1;

z=-inv（F3）;

y=z*F2;

whilek<=20

x=x0+y;

p=norm（x,inf）;q=norm（x0,inf）;

if（（p-q）/p

break

else

x0=x;

F2=[subs（f1,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）;

subs（f2,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）;

subs（f3,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）];

F3=[subs（g1,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g2,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g3,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）;

subs（g4,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g5,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g6,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）;

subs（g7,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g8,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）subs（g9,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）];

z=-inv（F3）;

y=z*F2;

x=x0+y;

k=k+1;

end

disp（x）;disp（k）;

最速下降法解非线性方程程序

symsx1x2x3t

f1=3*x1-cos（x2*x3）-（1/2）;

f2=（x1）^2-81*（x2+0.1）^2+sin（x3）+1.06;

f3=exp（-x1*x2）+20*x3+（（10*pi-3）/3）;

F=（f1）^2+（f2）^2+（f3）^2;

g1=diff（F,x1）;

g2=diff（F,x2）;

g3=diff（F,x3）;

G11=[g1;g2;g3];

G1=-G11;

x0=[0.1;0.1;-0.1];

G0=subs（G1,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）;

z0=x0+（t*G0）;

F1=subs（F,{x1,x2,x3},{z0

（1）,z0

（2）,z0（3）}）;

FD=diff（F1,t）;

xs0=solve（FD,t）;

FD_xs0=vpa（subs（FD,t,xs0）,6）;

F1_xs0=vpa（subs（F1,t,xs0）,6）;

t=xs0;

k=1;

Tol=10^（-5）;

whilek<=20

x=x0+（t*G0）;

p=norm（x,inf）;q=norm（x0,inf）;

if（（p-q）/p

break

else

x0=x;

G0=subs（G1,{x1,x2,x3},{x0

（1）,x0

（2）,x0（3）}）;

x=x0+t*G0;

k=k+1;

end

disp（x）;disp（k）;

（3）程序运行的结果

（1）将初值定为

分别运行上述两个程序，所得结果分别为：

牛顿法：

0.5000

0.0000

-0.5236

最速下降法：

0.12928443961888410941500498462202

.0770********

-0.52527958444225176826882692504685

（2）改变初值为

，分别运行上述两个程序，所得结果分别为：

牛顿法：

NaN

最速下降法：

20.219127745313790277417091321903

2.158********27343723864664313333

19.999648528574337249240101667477

（3）采用初值

，先用最速下降法求解3步，再用牛顿迭代法求解。

首先将第一个程序运行3次后，将结果作为第二个程序的初值再做迭代，运行结果如下：

0.4274

-0.0318

-0.4736

结果分析

Newton法的收敛速度较快，并且用相邻两次迭代的无穷范数相对误差来控制迭代次数，由程序运行的结果可知，改变初值后Newton法迭代的结果发散，由此可见Newton法对初值的要求较高。

而最速下降法，对初值的要求不高，它对任意的初值都是收敛的，因此，为了提高迭代的效率，应先用最速下降法迭代几次，来确定初值，再用牛顿迭代法进行迭代。

上机作业4

1.用规范的幂法与反幂法求矩阵A的按模最大，最小特征值与对应的特征向量。

e=10^-5

（1）幂法算法步骤

（1）取初始向量

（例如取

）,置精度要求

，置

。

（2）计算

（3）若

，则停止计算（

作为绝对值最大特征值

，

作为相应的特征向量）否则置

，转

（2）。

程序代码

function[l,v,s]=mifa（A,x0,eps）

ifnargin==2

eps=1.0e-5;

end

v=x0;

M=5000;

m=0;

l=0;

fork=1:

y=A*v;

m=max（y）;

v=y/m;

if（abs（m-l）

l=m;

s=k;

return

else

ifk==m

disp（'迭代步数太多，收敛速度太慢'）;

l=m;

s=M;

else

l=m;

end

运行结果

>>A=[4111;13-11;1-120;1102];

>>x0=[1111]';

>>[l,v,s]=mifa（A,x0）

结果如下：

特征向量

，特征值为

。

（2）反幂法

计算

的按模最小的特征值

的问题就是计算

的按模最大的特征值问题。

对于

应用幂法迭代（又称反幂法），可求得矩阵

的主特征值

，从而求得

的按模最小的特征值

。

反幂法的迭代公式为：

任取初始向量

，构造向量序列

迭代向量

可以通过解线性方程组

求得。

程序代码

function[l,v]=fanmifa（A,x0,eps）

ifnargin==2

eps=1.0e-5;

end

v=x0;

M=500;

m=0;

l=0;

fork=1:

y=A*v;

m=max（y）;

v=y/m;

if（abs（m-l）

l=1/m;

return

else

ifk==m

disp（'迭代步数太多，收敛速度太慢'）;

l=1/m;

s=M;

else

l=m;

end

运行结果

>>A=[4111;13-11;1-120;1102];

>>x0=[1111]';

>>[l,v]=fanmifa（A,x0）

结果如下：

特征向量

，特征值为

。

（3）结论

幂法程序可以用来计算矩阵绝对值最大的特征值及相应的特征向量。

幂法的缺点是开始的时候并不知道矩阵是否有单一的主特征值，也不知道如何选择

展开阅读全文