欧拉法解常微分方程.docx

欧拉法解常微分方程.docx

《欧拉法解常微分方程.docx》由会员分享，可在线阅读，更多相关《欧拉法解常微分方程.docx（15页珍藏版）》请在冰豆网上搜索。

欧拉法解常微分方程

纟沙理工久

数学与计算科学学院

实验报告

实验项目名称Eular方法求解一阶常微分方程数值解

所属课程名称偏微分方程数值解

实验类型验证性

班

学

姓

成

实验日期2015-3-26

级

号

名

绩

一、实验概述：

【实验目的】

熟练掌握应用显性Eular法和隐式Eular法求解一般一阶常微分方程的近似数值解。

【实验原理】

虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程。

求解从实际问题当中归结出来的微分方程主要靠数值解法。

欧拉方法是一类重要的数值解法。

这类方法回避解y（x）的函数表达式，而是寻求它在一系

列离散节点上的近似值，相邻的两个节点的间距称作步长。

假定步长为定数。

欧拉方法是一类离散化方法，这类方法将寻求解y（x）的分析问题转化为计算离

散值值的代数问题，从而使问题获得了实质性的简化。

然而随之带来的困难是，由于数据量往往很大，差分方法所归结出的可能是个大规模的代数方程组。

【实验环境】

1.硬件环境

22软件环境

MATLAB7.0

、实验内容：

【实验过程】（实验步骤）

（一）实验任务

描述某种化学反应过程的方程，利用显性和隐形Eualar方法求解下列一阶线性

微分方程组的近似数值解:

4

0.04y110y1y2

0.04%104y-iy231O7y2

yi（0）5（0）0,ya（0）0

（二）求解过程

Eular方法：

一阶线性微分方程初值问题

y'f（x,y）,axb

y（a）y。

ax0x.,....xnb

（1）

Xnxnh,h为步长

方程离散化：

差分和差商

y'gy1y0y1y0

x1x0h

愀必）y/0

h

y1y°hf（x°,y°）

（2）

yn1ynhf（X°,y°）

通过初始值y，依据递推公式

（2）逐步算出Y1,Y2,....,yn就为显性的Eular方

法。

隐形Eular方法:

（3）

y1y°hf（X1,yJyn1ynhf（Xn1,yn1）

公式（3）即为隐式Eular公式

（3）程序算法

1.利用显式Eular法方求解

利用MATLA进行求解，编写脚本文件如下:

文件名：

hql.m

%显性Eular方法

f0=1;g0=O;zO=O

delta=0.01;

time=1;

t=0:

delta:

time;

f=zeros（size（t））;

g=zeros（size（t））;

z=zeros（size（t））;

f1=zeros（size（t））;

g1=zeros（size（t））;

z1=zeros（size（t））;

f

（1）=f0;

g

（1）=g0;

z

（1）=z0;

fori=2:

length（t）

f1（i-1）=-0.04*f（i-1）+10000*f（i-1）*g（i-1）;

f（i）=f（i-1）+f1（i-1）*delta;

g1（i-1）=0.04*f（i-1）-10000*f（i-1）*g（i-1）-3*10A7*g（i-1）A2;

g（i）=g（i-1）+g1（i-1）*delta;

z1（i-1）=3*10A7*g（i-1）A2;

z（i）=z（i-1）+z（i-1）*delta;

Fun=f+g+zend

figure

plot（t,f.

'o'

）;

xlabel（

't'）;

ylabel（

'y1'

）;

title（

't-y1

变化图'）

figure

plot（t,g,

'o'

）;

xlabel（

't'）;

ylabel（

'y2'

）;

title（

't-y2

变化图'）

figure

plot（t,z,

'o'

）;

xlabel（

't'）;

ylabel（

'y3'

）;

title（

't-y3

变化图'）

figure

plot（t,Fun）;

xlabel（

't'）;

ylabel（

'y1+y2+y3'）;

title（

't-y1+y2+y3变化图'）

【实验结论】

A步长h=0.001时进行数据测试。

结果如下:

迭代第一次时,

Columns12through22

结果与方程描述内容相符。

迭代第二次时，

结果与方程描述内容基本相符。

迭代三次时,

A

ODe-a

1J.

Fun

Columns1throueh.11

匚olumns12through22

1.0000

1.oooo

i.aooo

0.&959

0.9803

结果与方程描述内容基本相符。

迭代1000次时,

0.0000

①OOOO-

山OOM

-7.LEh5

-r-Hf

Culmins23thxcr_i呂h33

KJi

3dtJ

]thxQLUhIL

atdtSAfitCisluHTiA34thtou£h44

0.0wo

1.oooo

■。

上DM

-'hoc»o

-a.OCQQi

■mown

-O,OGOQ

T・0000

riH

KaX

ZlaR

NaX

Nall

Eli&N

XsCl

JU誌

Euil

II曲

JIJ

El

模拟结果已经严重脱离事实，故当选择

delta为0.001时，该迭代方法不收敛。

时间与个变量直接的变化关系如图所示:

（k总:

匸:

越巫:

-1-

-?

■6

3I」_』L』L■!

L11

X00010020030U&050C60070C800901

.-233

xW

-3

^nrirLA

-€

iil」I」・」L1

0.010.020030J4ffl.05Q.06007Q.OSC.090.1

从上述图形可以明显看出，在迭代的不断进行时，各变量与时间的变化越来越大，且严重脱离了方程所描述的现实意义。

B.当选择h=0.00000001时，模拟结果如下:

迭代第一次，

Fun=

Columns

1thraxifti18

I

ioa

0

Cclvmns

19through36

Q

QQQ

Q

与A中结果相同。

迭代第二次，

1throughL8

11

19through35

值

I.9000e-

1x101创

1

Jxl01乱

1x11dan

并未产生误差。

跌二次迭代结果明显优于一中。

跌三次迭代结果,

pg=

Column?

1'through1超

111100

匚olunns19through36

地1000次迭代结果,

Tim•

匚z-lums1

itLTC-Lfih1'J

:

L兀i：

」

L"jQJO

KQW0

JLQCOO

i.o^oq

i.:

]b皿:

L.oaao

匚^lUtXilr11

lk£.srU|li药

3.（009

L30A-3

lboqw

i.aw

1.MDO

l.'330fl

1■旳如

].303：

'

mm

I.OOCO

C-El'jrns2j.

thrnufh

:

「

A

LOlOOl

:

LMftfl

1-OO'JCi

i.d'X'O

i._0（-06

i.oaaQ

j.oooa

J-OOliD

LOfl-Jtl

i.tura

（3

“luw31

1hxaugjq

c

4

LMO''

L0003

loom

b"«0

lomo

hOQQ・

hQOM

IrlCiJ

du

IIQ应

C#lune11

ikiTfoe^

0

In

hCWQ

].3000

B.QOOO

i.awn

1.D4D0

L.anaa

1,0000

].OTW

loo^d

I.OGOO

gw

*

txC-lijins5d

rhxtue?

i

结果明显是收敛的。

时间与个变量直接的变化关系如图所示:

卜界*2十/］变叱图

2Q

从图中能够清晰看出，当h=0.00000001时，模拟结果与方程所表示的显示意义相吻合。

说

明了显性Eualr方法的收敛性是与步长的选择是相关。

这就对我们们选择步长造成了困难，由于选择的步长不合适有可能得出错误的结论。

【实验小结】（收获体会）

1、软件使用

在写MATLABS言的时候要深刻理解题的意图，整理好思绪再做题目，在我运算的过程中，h取值取得越小、越细微，曲线逼近的越好。

2、欧拉法的缺点

简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。

因此欧拉格式一般不用于实际计算。

3、实验感想

在这次上机实验中，我掌握了解决常微分方程的基本方法，同时学会使用计算机软件对两种不同方法得到的结果进行判断，对我们以后对数据进行分析很有帮助。

三、指导教师评语及成绩：

评语

评语等级

优

良

中

及格

不及格

1.实验报告按时完成，字迹清楚，文字叙述流畅，逻辑性强

2.实验方案设计合理

3.实验过程（实验步骤详细，记录完整，数据合理，分析透彻）

4实验结论正确.

成绩：

指导教师签名：

批阅日期：

附录：

源程序

程序1:

%显性Eular方法

f0=1;g0=O;zO=O

delta=0.00000001;

time=0.00001;

t=O:

delta:

time;

f=zeros（size（t））;

g=zeros（size（t））;

z=zeros（size（t））;

f1=zeros（size（t））;

g1=zeros（size（t））;

z1=zeros（size（t））;

f

（1）=f0;

g

（1）=g0;

z

（1）=z0;

fori=2:

length（t）

f1（i-1）=-0.04*f（i-1）+10000*f（i-1）*g（i-1）;

f（i）=f（i-1）+f1（i-1）*delta;

g1