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

计算方法 上机作业.docx

上传人：b****7 文档编号：10343542 上传时间：2023-02-10 格式：DOCX 页数：44 大小：120.09KB
下载 相关 举报
计算方法 上机作业.docx_第1页
第1页 / 共44页
计算方法 上机作业.docx_第2页
第2页 / 共44页
计算方法 上机作业.docx_第3页
第3页 / 共44页
计算方法 上机作业.docx_第4页
第4页 / 共44页
计算方法 上机作业.docx_第5页
第5页 / 共44页
点击查看更多>>
下载资源
资源描述

计算方法 上机作业.docx

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

计算方法 上机作业.docx

计算方法上机作业

计算方法上机作业：

要求：

1）抄题；2）迭代公式（初值）或简单原理；3）程序，结果；（打印）4）结果分析。

上机作业1

1.分别用不动点迭代与Newton法求解方程2x-ex+3=0的正根与负根。

解：

1.不动迭代法：

迭代公式：

pn+1=g（pn），n=0,1,……

正根：

选取的迭代格式为x（k+1）=ln（2x（k）+3）

编写程序如下：

x0=0;

N=2000;Tol=10^（-4）;

k=1;

whilek<=N

x=log（2*x0+3）;

ifabs（x-x0）

k=k+1;x0=x;

end

disp（x）;disp（k）

结果：

1.9239

10

负根：

选取的迭代格式为x（k+1）=[exp（x（k））-3]/2

编写程序如下：

x0=0;

N=2000;Tol=10^（-4）;

k=1;

whilek<=N

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

ifabs（x-x0）

k=k+1;x0=x;

end

disp（x）;disp（k）

结果：

-1.3734

7

2.Newton法：

其迭代格式为x（k+1）=x（k）-f（x（k））/f’，即x=（（1-x）e^x-3）/（2-e^x）

编写程序如下：

正根：

x0=3;

N=2000;Tol=10^（-7）;

fork=1:

N

x=（（1-x0）*exp（x0）-3）/（2-exp（x0））;

ifabs（x-x0）

x0=x;

end

disp（x）;disp（k）

结果：

1.9239

6

负根：

x0=0;

N=2000;Tol=10^（-7）;

fork=1:

N

x=（（1-x0）*exp（x0）-3）/（2-exp（x0））;

ifabs（x-x0）

x0=x;

end

disp（x）;disp（k）

结果：

-1.3734

5

结果分析：

不动点迭代法的原理比较简单。

相比较之下，牛顿迭代法收敛于根的速度更快。

 

2.用Newton法求解方程x-sinx＝0的根。

再用Steffensen’smethod加速其收敛。

解：

Newton法迭代公式：

，n=0,1,……

Steffensen加速原理：

程序：

Newton法：

x0=-100,x1=100;

N=2000;Tol=10^（-9）;

fork=1:

N

x=x1-（x1-sin（x1））*（x1-x0）/（x1-sin（x1）-x0+sin（x0））;

ifabs（x-x1）

x0=x1;x1=x;

end

disp（x）;disp（k）

运行结果：

x0=

-100

1.4211e-014

2

Steffensen’smethod加速：

M文件：

functions=steff（f,x0,eps,max）

f=inline（f）;

x

（1）=x0;

disp（'kxyz'）;

fork=1:

max

y（k）=feval（f,x（k））;

z（k）=feval（f,y（k））;

x（k+1）=x（k）-（y（k）-x（k））^2/（z（k）-2*y（k）+x（k））;

ifabs（x（k+1）-x（k））

break

end

disp（sprintf（'%d%f%f%f',k,x（k）,y（k）,z（k）））;

end

s=x（k+1）;

 

命令窗口：

s=steff（'x-sin（x）',1,1.0e-7,200）

运行结果如下：

kxyz

11.0000000.1585290.000663

2-0.035793-0.000008-0.000000

s=

2.0680e-025

结果分析：

Steffensen加速极大地改善了原迭代的收敛性质，它加速了牛顿法的收敛速度，且结果更为精确，但有些时候计算量可能会有所增加。

上机作业2

分别用Jacobi、Seidel、SOR（w=1.1,1.2,1.3,1.4,1.5）方法求解方程组Ax=b，这里

x0=[0….0]’；e=10-5

解：

Jacobi迭代：

公式：

x（k+1）=（I-D-1A）x（k）+D-1b

程序：

A=[-12111111111;1-1211111111;11-121111111;111-12111111;1111-1211111;11111-121111;111111-12111;1111111-1211;11111111-121;111111111-12];

b=[-1-1-1-1-1-1-1-1-1-1]';N=500;eps=1e-5;

n=length（b）;x0=zeros（n,1）;x=zeros（n,1）;k=0;D=zeros（n,n）;

fori=1:

n

forj=1:

n

ifi==j

D（i,j）=A（i,j）;

end

end

end

BJ=eye（n）-inv（D）*A;dJ=inv（D）*b

whilek

x=BJ*x0+dJ;

ifnorm（x-x0,inf）

x0=x;k=k+1;

end

ifk==N,warning（'已达到迭代次数上限'）;end

disp（['k=',num2str（k）]）,x

运行结果如下：

dJ=

0.0833

0.0833

0.0833

0.0833

0.0833

0.0833

0.0833

0.0833

0.0833

0.0833

k=32

x=

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

故x=[0.33330.33330.33330.33330.33330.33330.33330.33330.33330.3333]’

Gauss-Seidel迭代:

公式：

x（k+1）=[-（L+D）-1U]x（k）+（L+D）-1b

程序：

A=[-12111111111;1-1211111111;11-121111111;111-12111111;1111-1211111;11111-121111;111111-12111;1111111-1211;11111111-121;111111111-12];

b=[-1-1-1-1-1-1-1-1-1-1]';N=500;ep=1e-5;

n=length（b）;x0=zeros（n,1）;

x=zeros（n,1）;k=0;

whilek

fori=1:

n

ifi==1x

（1）=（b

（1）-A（1,2:

n）*x0（2:

n））/A（1,1）;

elseifi==n

x（n）=（b（n）-A（n,1:

n-1）*x（1:

n-1））/A（n,n）;

else

x（i）=（b（i）-A（i,1:

i-1）*x（1:

i-1）-A（i,i+1:

n）*x0（i+1:

n））/A（i,i）;

end

end

end

ifnorm（x-x0,inf）

x0=x;k=k+1;end

ifk==N,warning（'已达到迭代次数上限'）;end

disp（['k=',num2str（k）]）,x

运行结果如下：

k=18

x=

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

故x=[0.33330.33330.33330.33330.33330.33330.33330.33330.33330.3333]’

SOR迭代：

公式：

x（k+1）=（D+ωL）-1[（1-ω）D-ωU]x（k）+ω（D+ωL）-1b

程序：

w=1.1

A=[-12111111111;1-1211111111;11-121111111;111-12111111;1111-1211111;11111-121111;111111-12111;1111111-1211;11111111-121;111111111-12];

b=[-1-1-1-1-1-1-1-1-1-1]'

n=length（b）;N=500;eps=1e-5;x0=zeros（n,1）;omega=1.1;x=zeros（n,1）;k=0;

whilek

fori=1:

n

ifi==1

x1

（1）=（b

（1）-A（1,2:

n）*x0（2:

n））/A（1,1）;

elseifi==n

x1（n）=（b（n）-A（n,1:

n-1）*x（1:

n-1））/A（n,n）;

else

x1（i）=（b（i）-A（i,1:

i-1）*x（1:

i-1）-A（i,i+1:

n）*x0（i+1:

n））/A（i,i）;

end

end

x（i）=（1-omega）*x0（i）+omega*x1（i）;

end

ifnorm（x0-x,inf）

k=k+1;x0=x;

end

ifk==N,warning（'已达到迭代次数上限'）;end

disp（['k=',num2str（k）]）,x

运行结果如下：

b=

-1

-1

-1

-1

-1

-1

-1

-1

-1

-1

k=15

x=

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

故x=[0.33330.33330.33330.33330.33330.33330.33330.33330.33330.3333]’

w=1.2

A=[-12111111111;1-1211111111;11-121111111;111-12111111;1111-1211111;11111-121111;111111-12111;1111111-1211;11111111-121;111111111-12];

b=[-1-1-1-1-1-1-1-1-1-1]'

n=length（b）;N=500;eps=1e-5;x0=zeros（n,1）;omega=1.2;x=zeros（n,1）;k=0;

whilek

fori=1:

n

ifi==1

x1

（1）=（b

（1）-A（1,2:

n）*x0（2:

n））/A（1,1）;

elseifi==n

x1（n）=（b（n）-A（n,1:

n-1）*x（1:

n-1））/A（n,n）;

else

x1（i）=（b（i）-A（i,1:

i-1）*x（1:

i-1）-A（i,i+1:

n）*x0（i+1:

n））/A（i,i）;

end

end

x（i）=（1-omega）*x0（i）+omega*x1（i）;

end

ifnorm（x0-x,inf）

k=k+1;x0=x;

end

ifk==N,warning（'已达到迭代次数上限'）;end

disp（['k=',num2str（k）]）,x

运行结果如下：

b=

-1

-1

-1

-1

-1

-1

-1

-1

-1

-1

k=11

x=

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

故x=[0.33330.33330.33330.33330.33330.33330.33330.33330.33330.3333]’

w=1.3

A=[-12111111111;1-1211111111;11-121111111;111-12111111;1111-1211111;11111-121111;111111-12111;1111111-1211;11111111-121;111111111-12];

b=[-1-1-1-1-1-1-1-1-1-1]'

n=length（b）;N=500;eps=1e-5;x0=zeros（n,1）;omega=1.3;x=zeros（n,1）;k=0;

whilek

fori=1:

n

ifi==1

x1

（1）=（b

（1）-A（1,2:

n）*x0（2:

n））/A（1,1）;

elseifi==n

x1（n）=（b（n）-A（n,1:

n-1）*x（1:

n-1））/A（n,n）;

else

x1（i）=（b（i）-A（i,1:

i-1）*x（1:

i-1）-A（i,i+1:

n）*x0（i+1:

n））/A（i,i）;

end

end

x（i）=（1-omega）*x0（i）+omega*x1（i）;

end

ifnorm（x0-x,inf）

k=k+1;x0=x;

end

ifk==N,warning（'已达到迭代次数上限'）;end

disp（['k=',num2str（k）]）,x

运行结果如下：

b=

-1

-1

-1

-1

-1

-1

-1

-1

-1

-1

k=9

x=

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

故x=[0.33330.33330.33330.33330.33330.33330.33330.33330.33330.3333]’

w=1.4

A=[-12111111111;1-1211111111;11-121111111;111-12111111;1111-1211111;11111-121111;111111-12111;1111111-1211;11111111-121;111111111-12];

b=[-1-1-1-1-1-1-1-1-1-1]'

n=length（b）;N=500;eps=1e-5;x0=zeros（n,1）;omega=1.4;x=zeros（n,1）;k=0;

whilek

fori=1:

n

ifi==1

x1

（1）=（b

（1）-A（1,2:

n）*x0（2:

n））/A（1,1）;

elseifi==n

x1（n）=（b（n）-A（n,1:

n-1）*x（1:

n-1））/A（n,n）;

else

x1（i）=（b（i）-A（i,1:

i-1）*x（1:

i-1）-A（i,i+1:

n）*x0（i+1:

n））/A（i,i）;

end

end

x（i）=（1-omega）*x0（i）+omega*x1（i）;

end

ifnorm（x0-x,inf）

k=k+1;x0=x;

end

ifk==N,warning（'已达到迭代次数上限'）;end

disp（['k=',num2str（k）]）,x

运行结果如下：

b=

-1

-1

-1

-1

-1

-1

-1

-1

-1

-1

k=11

x=

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

故x=[0.33330.33330.33330.33330.33330.33330.33330.33330.33330.3333]’

w=1.5

A=[-12111111111;1-1211111111;11-121111111;111-12111111;1111-1211111;11111-121111;111111-12111;1111111-1211;11111111-121;111111111-12];

b=[-1-1-1-1-1-1-1-1-1-1]'

n=length（b）;N=500;eps=1e-5;x0=zeros（n,1）;omega=1.5;x=zeros（n,1）;k=0;

whilek

fori=1:

n

ifi==1

x1

（1）=（b

（1）-A（1,2:

n）*x0（2:

n））/A（1,1）;

elseifi==n

x1（n）=（b（n）-A（n,1:

n-1）*x（1:

n-1））/A（n,n）;

else

x1（i）=（b（i）-A（i,1:

i-1）*x（1:

i-1）-A（i,i+1:

n）*x0（i+1:

n））/A（i,i）;

end

end

x（i）=（1-omega）*x0（i）+omega*x1（i）;

end

ifnorm（x0-x,inf）

k=k+1;x0=x;

end

ifk==N,warning（'已达到迭代次数上限'）;end

disp（['k=',num2str（k）]）,x

运行结果如下：

b=

-1

-1

-1

-1

-1

-1

-1

-1

-1

-1

k=15

x=

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

0.3333

故x=[0.33330.33330.33330.33330.33330.33330.33330.33330.33330.3333]’

三种迭代方法对比如下表：

 

表1：

三种迭代方法对比表

迭代方法

迭代次数

近似解

Jacobi迭代法

32

x=[0.33330.33330.33330.33330.33330.33330.33330.33330.33330.3333]’

Gauss-Seidel迭代法

18

SOR迭代法

w=1.1

15

w=1.2

11

w=1.3

9

w=1.4

11

w=1.5

15

结果分析：

由以上计算过程可以看出，这三个迭代法都是收敛的，但不同的是Jacobi迭代法的收敛速度是最慢的，SOR迭代法收的敛速度是最快，且随着松弛因子w的变化而不同，从表中可以看出，当w=1.3的时候收敛速度最快。

方程的近似解为x=[0.33330.33330.33330.33330.33330.33330.33330.33330.33330.3333]’。

上机作业3

用Newton法与最速下降法求方程组

x2+2y2-2=0

x2=y

在（0.8,0.7）附近的根。

解：

Newton法：

迭代公式：

x（k+1）=x（k）-[F’（x（k））]-1F（x（k））

程序：

x0=[0.8;0.7];

fork=1:

10

z=-[2*x0

（1）,4*x0

（2）;2*x0

（1）,-1]\[x0

（1）^2+2*x0

（2）^2-2;x0

（1）^2-x0

（2）];

x=x0+z;

ifnorm（z）<1.0e-6;

disp（x）;disp（k）

return;

end

x0=x;

end

运行结果：

0.8836

0.7808

4

最速下降法：

原理：

对于方程组

构造函数

，

，求

的零最小值；即求

，使得

=0。

也是方程组的解。

M文件：

function[r,k]=FastDown（F,x0,h,eps）

formatlong;

ifnargin==3

eps=1.0e-8;

end

n=length（x0）;

x0=transpose（x0）;

k=1;

tol=1;

whiletol>eps

fx=subs（F,findsym（F）,x0）;

J=zeros（n,n）;

fori=1:

n

x1=x0;

x1（i）=x1（i）+h;

J（:

i）=（subs（F,findsym（F）,x1）-fx）/h;

end

lamda=fx/sum（diag（transpose（J）*J））;

r=x0-J*lamda;

fr=subs（F,findsym（F）,r）;

tol=dot（fr,fr）;

x0=r;

k=k+1;

if（k>500）

disp（'迭代步数太多,可能不收敛'）;

break;

end

end

formatshort;

命令窗口：

symsxy;

f=[x^2+2*y^2-2;x^2-y];

[r,k]=FastDown（f,[0.8,0.7],1.0e-6）

运行结果如下：

r=

0.8836

0.7807

 

k=

19

结果分析：

Newton法具有较快的收敛速度，缺点是对初值要求过高，计算量较大；最速下降法对任意初值总是收敛的，缺点是收敛速度较慢，通常是越靠近解，逼近速度越慢。

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

当前位置：首页 > 高等教育 > 军事

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

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