haihongyuan.com
海量文库 文档专家
全站搜索:
您现在的位置:首页 > 小学教育 > 小学数学小学数学

数学实验4课件

发布时间:2014-01-07 10:43:15  

数学实验
---微分方程模型

常微分方程的解
? 微分方程的解析解

通解 初边值问题 常微分方程组 非刚性常微分方程 刚性常微分方程 高阶微分方程

? 微分方程的数值解

1、常微分方程解析解

dsolve

?在Matlab中,符号运算工具箱(Symbolic Math Toolbox ),提 供了功能强大的求解常微分方程的符号运算命令dsolve.

?常微分方程在Matlab中按如下规定重新表达: 符号D表示对变量的求导。Dy表示对变量y求一阶导数, 当需要求变量的n阶导数时,用Dn表示,D4y表示对变量 y求4阶导数。 任何D后所跟的字母为因变量,自变量可以指定或由 系统规则选定为缺省.

求解常微分方程的通解
?

无初边值条件的常微分方程的解就是该方程的通 解.其使用格式为: dsolve('diff_equation') dsolve(' diff_equation','var')

式中diff_equation 为待解的常微分, 第1种格式将 以变量t为自变量进行求解, 第2种格式则需定义 自变量var.

求解常微分方程的通解
du 例1求 ? 1 ? u 2 的通解. dt



输入命令:dsolve('Du=1+u^2','t')

结 果:u = tg(t+c1)
To Matlab(ff1)default

求解常微分方程的通解
?

例 1’

2 试解常微分方程 x ? y ? ( x ? 2 y ) y ' ? 0

解:编写程序如下(liti1. m): diff_equ='x^2+y+(x-2*y)*Dy=0'; y= dsolve(diff_equ,'x') %pretty(y) 结果为:x/2 + ((4*x^3)/3 + x^2 + C2)^(1/2)/2
?

x/2 -((4*x^3)/3 + x^2 + C2)^(1/2)/2

求解常微分方程的初边值问题
求解带有初边值条件的常微分方程的使用格式为: dsolve(‘diff_equation’,‘condition1,condition2,…','var'), 其中condition1,condition2,… 即为微分方程的初边值条件.

例 2 求微分方程的特解.

?d 2 y dy ? 2 ? 4 ? 29 y ? 0 ? dx dx ? y (0) ? 0, y ' (0) ? 15 ? 解:编写程序如下(liti2. m): y=dsolve('D2y+4*Dy+29*y=0', 'y(0)=0,Dy(0)=15','x')

结果为:y =3e-2xsin(5x)

求解常微分方程组
? 求解常微分方程组的命令格式为:

? dsolve(‘diff_equ1,diff_equ2,…‘,

'var') ? dsolve('diff_equ1,diff_equ2,…‘,'condition1,co ndition2,…','var')
?

第1种格式用于求解方程(组)的通解,第2种格式可 以加上初边值条件,用于具体求解。

? f ' '?3g ? sin x ? 例3 试求常微分方程组:? ? g '? f ' ? cos x
的通解和在初边值条件为 f ' (2) ? 0, f (3) ? 3, g (5) ? 1 的解。
?

解 编写程序如下(liti3.m): clc,clear equ1='D2f+3*g=sin(x)'; equ2='Dg+Df=cos(x)'; [general_f,general_g]=dsolve(equ1,equ2,'x') [f,g]=dsolve(equ1,equ2,'Df(2)=0,f(3)=3,g(5)=1','x')

练习
1、求解 y '''? y '' ? x, y(1) ? 8, y '(1) ? 7, y ''(2) ? 4 y=dsolve('D3y-D2y=x','y(1)=8,Dy(1)=7,D2y(2)=4','x') 2、求解微分方程组
? dx ? dt ? 2 x ? 3 y ? 3z ? dy ? ? ? 4 x ? 5 y ? 3z ? dt ? dz ? 4 x ? 4 y ? 2 z ? dt ?

的通解,

[x,y,z

]=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z',

'Dz=4*x-4*y+2*z', 't');

2、常微分方程数值解 非刚性常微分方程 刚性常微分方程 高阶微分方程

微分方程的数值解
建立微分方程只是解决问题的第一步,通常需 要求出方程的解来说明实际现象,并加以检验.如果 能得到解析形式的解固然是便于分析和应用的.但 是我们知道,只有线性常系数微分方程,并且自由项 是某些特殊类型的函数时,才可以肯定得到这样的 解,而绝大多数变系数方程、非线性方程都是所谓 “解不出来”的,即使看起来非常简单的方程,于是 对于用微分方程解决实际问题来说,数值解法就是 一个十分重要的手段.

刚性常微分方程(组)
在化学工程及自动控制等领域中,所涉及的常微分方程组 初值问题常常是所谓的“刚性”问题。 即,对一阶线性微分方程组 dy (*) ? Ay ? ?(x ) dx 其中 y, ? ? R m , A 为 m 阶方阵。 若矩阵 A 的特征值 ?i (i ? 1,2,?, m) 满足关系 Re ?i ? 0 (i ? 1,2,?, m) max | Re ?i |?? min | Re ?i |
1? i ? m 1? i ? m

则称方程组(*)为刚性方程组或 Stiff 方程组,称数 s ? max | Re ?i | / min | Re ?i |
1? i ? m 1? i ? m

为刚性比。 对刚性方程组,用前面所介绍的数值方法求解,都会遇到本质上 的困难,这是由数值方法本身的稳定性限制所决定的。 理论上的分析表明,求解刚性问题所选用的数值方法最好是对步 长 h 不作任何限制。

用Matlab软件求常微分方程的数值解

[t,x]=solver(’f’,ts,x0,options)
自变 量值 函数 值

ode45 ode23 ode113 ode15s ode23s

由待解 方程写 成的m文件名

ts=[t0,tf], t0、tf为自 变量的初 值和终值

函数的 初值

ode23:组合的2/3阶龙格-库塔-芬尔格算法 ode45:运用组合的4/5阶龙格-库塔-芬尔格算法
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为:options=odeset(’reltol’,rt,’abstol’,at), rt,at:分别为设定的相对误差和绝对误差.

常微分方程求解
(不能确定是否是刚性常微分方程时,首先用命令 ode45 ,然后用命令ode15S.)

ode45 ★ 微分方程高阶数值解法,基于显式龙格-库塔(4,5)法(四阶精度,五阶误 ode23 ode113

差),采用单步算法来计算 微分方程低阶数值解法,这是一个比ode45低阶的方法,基于显式龙格库塔(2,3)法 用于更高阶或大的标量计算。采用多步法。 用于解决难度适中的问题.

ode23t

ode15s ★ 与ode23相同,但要求的精度更高。采用数值差分方法。为多步法

ode23s
ode23tb

用于解决难度较大的微分方程组。对于系统中存在常量矩阵的情况也 有用。采用单步法 用于解决难度较大的问题,对于系统中存在常量矩阵的情况也有用.

? y ' ? 2x 例题1 求初值问

题 ? 的数值解. ? y ( ?1) ? 4

? y ? x2 ? c 2 解:显然,解析解存在,为 ? ? y ? x ?3 ? y (?1) ? 4 下面,求其数值解.
1建立函数文件fun1.m function dy=fun1(x,y) dy=2*x; 2主程序shuzhijie.m x0=[-1,1]; [X,Y]=ode45('fun1',x0,4); plot(X,Y,’+’)
4 3.9 3.8 3.7 3.6 3.5 3.4 3.3 3.2 3.1 3 -1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

? y1 ' ? y2 y3 ?y ' ? ?y y ? 2 1 3 例2 解微分方程组. ? ? y3 ' ? ?0.51y1 y2 ? y1 (0) ? 0, y2 (0) ? 1, y3 (0) ? 1 ?
解:1、建立m-文件rigid.m如下: function dy=rigid(t,y) dy=[y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)]
1 0.8 0.6 0.4 0.2

2、取t0=0,tf=12,输入命令: [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
0 -0.2 -0.4 -0.6 -0.8 -1 0 2 4 6 8 10 12

?

高阶微分方程的初值问题可以通过变量代换化为一阶微 分方程组初值问题。 设有m阶常微分方程初值问题 ? y ( m ) ? f ( x, y , y ' ,?, y ( m ?1) ) a ? x ? b ? ( ( y ( a ) ? y 0 , y ' ( a ) ? y 01) ,?, y ( m ?1) ( a ) ? y 0m ?1) ? 引入新变量
y1 ? y, y 2 ? y ' ,?, y m ? y ( m ?1)

化为一阶微分方程初值问题

? y '1 ? y 2 ? ? y'2 ? y3 ? ? ? ? y' ? y m ? m ?1 ? y ' m ? f ( x, y1 ,?, y m ) ?

y1 ( a ) ? y 0
( y 2 ( a ) ? y 01)

?
( y m ?1 ( a ) ? y 0m ? 2 ) ( y m ( a ) ? y 0m ?1)

例3 :考虑高阶微分方程问题

y '''? 3 y ''? y ' y ? 0, y (0) ? 0, y '(0) ? 1, y ''(0) ? ? 1
解:(1)设y1 ? y, y2 ? y ', y3 ? y '',则有

? y1 ' ? y 2 ? ? y ' 2 ? y3 ? y' ? 3 y ? y y 3 2 1 ? 3
于是初值问题变为

y1 (0) ? 0 y 2 ( 0) ? 1 y 3 ( 0 ) ? ?1

Y ' ? F (t , Y ), Y (0) ? Y0,其中,Y ? [ y1 ; y2 ; y3 ]

(2)把一阶方程组写成接受两个参数 t 和 y 返回一个 列向量的M文件djl.m: function dy=djl(t,y); dy=[y(2);y(3);3*y(3)+y(2)*y(1)]; (或者dy=zeros(3,1);dy(1)= y(2);dy(2)= y(3); dy(3)=3*y(3)+y(2)*y(1)];) (3)取t0=0,tf=1,输入命令: [T,Y]=ode45(‘djl',[0 1],[0 1 -1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+') legend(‘y’, ‘y一阶导数’,‘ y二阶导数 ')
5 0 -5 -10 -15 -20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8

y y一 阶 导 数 y二 阶 导 数

0.9

1

?

例4 求 van der Pol 方程

y ' '? ? (1 ? y 2 ) y '? y ? 0
? 0 一参数。

的数值解,这里是 ?

解 (i)化成常微分方程组。 设:

y1 ? y, y2 ? y '

则有 :

? y '1 ? y 2 ? y ' 2 ? ? (1 ? y12 ) y 2 ? y1 ?

?

(ii)书写M文件(对于

)vdp1.m:

function dy=vdp1(t,y); dy=[y(2);(1-y(1)^2)*y(2)-y(1)]; ? (iii)调用Matlab函数(liti4)。对于初值 y(0) ? 2, y' (0) ? 0 ,解为 [T,Y]=ode45('vdp1',[0 20],[2;0]); ? (iv)观察结果。利用图形输出解的结果: plot(T,Y(:,1),'-',T,Y(:,2),'--') title('Solution of van der Pol quation,mu=1'); xlabel('time t'); yl

abel('solution y'); legend('y1','y2');

Solution of van der Pol Equation,mu=1 3 y1 y2 2

1
solution y

0

-1

-2

-3

0

2

4

6

8

10 time t

12

14

16

18

20

? 例5
?

van der Pol 方程,

(刚性)

(i)书写M文件vdp1000.m: function dy=vdp1000(t,y); dy=[y(2);1000*(1-y(1)^2)*y(2)-y(1)]; (ii)观察结果 [t,y]=ode15s('vdp1000',[0 3000],[2;0]); plot(t,y(:,1),'o') title('Solution of van der Pol Equation,mu=1000'); xlabel('time t'); ylabel('solution y(:,1)');

?

Solution of van der Pol Equation,mu=1000 2 1.5 1 0.5

solution y(:,1)

0 -0.5 -1 -1.5 -2 -2.5

0

500

1000

1500 time t

2000

2500

3000


网站首页网站地图 站长统计
All rights reserved Powered by 海文库
copyright ©right 2010-2011。
文档资料库内容来自网络,如有侵犯请联系客服。zhit326@126.com