实验报告七常微分方程初值问答的数值解法Word文件下载.docx

实验报告七常微分方程初值问答的数值解法Word文件下载.docx

《实验报告七常微分方程初值问答的数值解法Word文件下载.docx》由会员分享，可在线阅读，更多相关《实验报告七常微分方程初值问答的数值解法Word文件下载.docx（15页珍藏版）》请在冰豆网上搜索。

用以下三种不同的方法求下述微分方程的数值解，取

画出解的图形，与精确值比较并进行分析。

1）欧拉法；

2）改进欧拉法；

3）龙格－库塔方法；

2-4分析应用题

考虑一个涉及到社会上与众不同的人的繁衍问题模型。

假设在时刻（单位为年），社会上有人口人，又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人。

而固定比例为的所有其他的后代也是与众不同的人。

如果对所有人来说出生率假定为常数，又如果普通的人和与众不同的人的婚配是任意的，则此问题可以用微分方程表示为：

其中变量表示在时刻社会上与众不同的人的比例，表示在时刻人口中与众不同的人的数量。

1）假定和，当步长为年时，求从到解的近似值，并作出近似解的曲线图形。

2）精确求出微分方程的解，并将你当时在分题（b）中得到的结果与此时的精确值进行比较。

【MATLAB相关函数】

⏹求微分方程的解析解及其数值的代入

dsolve（‘egn1’，‘egn2’，‘’）

subs（expr,{x,y,…},{x1,y1,…}）

其中‘egn’表示第个方程，‘’表示微分方程中的自变量，默认时自变量为。

subs命令中的expr、x、y为符合型表达式，x、y分别用数值x1、x2代入。

>

symsxyz

subs（'

x+y+z'

{x,y,z},{1,2,3}）

ans=

6

symsx

x^2'

x,2）

4

s=dsolve（‘’,‘’,‘’）

ans=

subs（s,x,2）

-0.3721

⏹右端函数的自动生成

f=inline（‘expr’,’var1’,‘var2’,……）

其中’expr’表示函数的表达式，’var1’,‘var2’表示函数表达式中的变量，运行该函数，生成一个新的函数表达式为f（var1,var2,……）。

f=inline（'

x+3*y'

'

x'

y'

）

f=

Inlinefunction:

f（x,y）=x+3*y

f（2,3）

11

⏹4，5阶龙格－库塔方法求解微分方程数值解

[t,x]=ode45（f,ts,x0,options）

其中f是由待解方程写成的m文件名；

x0为函数的初值；

t,x分别为输出的自变量和函数值（列向量），t的步长是程序根据误差限自动选定的。

若ts=[t0,t1,t2,…,tf]，则输出在自变量指定值，等步长时用ts=t0:

k:

tf，输出在等分点；

options用于设定误差限（可以缺省，缺省时设定为相对误差，绝对误差），程序为：

options=odeset（‘reltol’,rt,’abstol’,at），这里rt,at分别为设定的相对误差和绝对误差。

常用选项见下表。

选项名

功能

可选值

省缺值

AbsTol

设定绝对误差

正数

RelTol

设定相对误差

InitialStep

设定初始步长

自动

MaxStep

设定步长上界

MaxOrder

设定ode15s的最高阶数

1,2,3,4,5

5

Stats

显示计算成本统计

on,off

off

BDF

设定ode15s是否用反向差分

例：

解微分方程

在命令窗口执行

=（‘’,‘’,‘’）；

；

ans=

01.0000

0.0502　1.0490

0.10051.0959

0.15071.1408

3.85072.9503

3.90052.9672

3.95022.9839

4.00003.0006

plot（,,‘o-’,）%解函数图形表示

%不用输出变量，则直接输出图形

1.0000　1.7321

2.00002.2361

3.00002.6458

三.操作方法与实验步骤（包括实验数据记录和处理）

Euler法

y=euler（a,b,n,y0,f,f1,b1）

y=zeros（1,n+1）;

y

（1）=y0;

h=（b-a）/n;

x=a:

h:

b;

fori=1:

n;

y（i+1）=y（i）+h*f（x（i）,y（i））;

end

plot（x,y）

holdon

%求微分方程的精确解

x1=linspace（a,b,100）;

'

精确解为'

s=dsolve（f1,b1,'

y1=zeros（1,100）;

for

i=1:

100

y1（i）=subs（s,x,x1（i））;

end

plot（x1,y1,'

r'

title（'

红色代表精确解'

改进Euler法

y=eulerpro（a,b,n,y0,f,f1,b1）

%求微分方程的数值解

y=zeros（1,n+1）;

y

（1）=y0;

h=（b-a）/n;

x=a:

for

T1=f（x（i）,y（i））;

T2=f（x（i+1）,y（i）+h*T1）;

y（i+1）=y（i）+（h/2）*（T1+T2）;

%求微分方程的精确解

'

y1=zeros（1,100）;

fori=1:

2-2分析应用题

（1）向前欧拉法

euler（0,10,100,10,inline（'

y-20'

）,'

Dy=y-20'

y（0）=10'

）

精确解为

s=

20-10*exp（x）

1.0e+005*

Columns1through8

0.00010.00010.00010.00010.00010.00000.0000

0.0000

Columns9through16

-0.0000-0.0000-0.0001-0.0001-0.0001-0.0001-0.0002

-0.0002

Columns17through24

-0.0003-0.0003-0.0004-0.0004-0.0005-0.0005-0.0006

-0.0007

Columns25through32

-0.0008-0.0009-0.0010-0.0011-0.0012-0.0014-0.0015

-0.0017

Columns33through40

-0.0019-0.0021-0.0024-0.0026-0.0029-0.0032-0.0035

-0.0039

Columns41through48

-0.0043-0.0048-0.0053-0.0058-0.0064-0.0071-0.0078

-0.0086

Columns49through56

-0.0095-0.0105-0.0115-0.0127-0.0140-0.0154-0.0170

-0.0187

Columns57through64

-0.0206-0.0227-0.0250-0.0275-0.0302-0.0333-0.0366

-0.0403

Columns65through72

-0.0444-0.0488-0.0537-0.0591-0.0651-0.0716-0.0788

-0.0867

Columns73through80

-0.0954-0.1049-0.1154-0.1270-0.1397-0.1537-0.1691

-0.1860

Columns81through88

-0.2046-0.2251-0.2477-0.2724-0.2997-0.3297-0.3627

-0.3990

Columns89through96

-0.4389-0.4828-0.5311-0.5842-0.6427-0.7070-0.7777

-0.8555

Columns97through101

-0.9410-1.0352-1.1387-1.2526-1.3779

（2）改进欧拉法

>

eulerpro（0,10,100,10,inline（'

ans=

精确解为

s=

20-10*exp（x）

1.0e+00