Matlab在常微分方程求解中的应用

ODE 指令基本用法

使用ODE指令时，必须先将要求解的ODE表示 成一个函数
输入为t(时间)及y(状态函数:
State Variables） 输出则为dy（状态变量的微分值）
ODE 指令基本用法

若ODE函数的文件为odeFile.m,则调用ODE指令 的格式如下： [t,y] = solver('odeFile',[t0,t1],y0)
Matlab在常微分方程 求解中的应用
实验目的
（1）学会用Matlab软件求解微分方程的初值问题 （2）了解微分方程数值解思想，掌握基本的微分方程数值解方法 （3）学会根据实际问题建立简单微分方程数学模型 （4）了解计算机数据仿真、数据模拟的基本方法
常微分方程
包含一个自变量和它的未知函数以及未知函数的导数的 等式 形成和发展与力学、天文学、物理学及其他自然科学技 术的发展互相促进和推动
函数的 初值
ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode45：运用组合的4/5阶龙格-库塔-芬尔格算法 用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差.
[t0,t1]
是积分的时间区间 y0代表起始条件（Initial Conditions） solver是前表所列的各种ODE指令 t是输出的时间向量 y是对应的状态变量向量
[t,y]=solver（’f’,ts,y0,options）
自变 函数 量值 值 ode45 ode23 由待解 ode113 方程写 ode15s ode23s 成的m文件名 ts=[t0 , tf]， t0、tf为自 变量的初 值和终值
结 果 为：x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = (c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t
虽然说解析解是最精确的，但是实际问 题中常要求研究常微分方程的数值解
y' F ( y, t )
y
为一个向量，代表状态变量

设 =1，ODE文件(vdp1.m)显示如下： >> type vdp1.m function dy = vdp1(t, y) mu = 1; dy = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];
结 果 为 : y =3e-2xsin（5x）
例3
求微分方程组的通解. 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')； x=simple(x) % 将x化简 y=simple(y) z=simple(z)
f x, y1 f x, y2 L y1 y2
则初值问题(1)的解 yx 存在并且唯一。
常微分方程的解析解
求微分方程（组）的解析解命令: dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’) 记号：在表达微分方程时，用字母D表示求微分，D2、D3等 表示求高阶微分。D后所跟字母为因变量，自变量可以指定 或由系统规则选定为缺省。 例如：微分方程 可以表示为D2y=0. du 1 u 2 的通解. 例1 求 dt
-
y ( xn 1 ) y ( xn ) h
启发：只要对曲线在区 间[ xn , xn 1 ]上的平均斜率提供一种 算法，就可以得到 一种计算yn 1的公式；如果设法在区 间[ xn , xn 1 ]内多预测几个点的斜率 值，然 由此可构造出由 yn计算yn 1的精度更高的计算公式
后将这些点处的斜率值 的加权平均值作为曲线 在区间[ xn , xn 1 ]上的平均斜率，
注意:
1、在解n个未知函数的方程组时，y0和y均为n维向量， m-文件中的待解方程组应以y的分量形式写成. 2、使用Matlab软件求数值解时，高阶微分方程必须 等价地变换成一阶微分方程组.
Van der Pol微分方程

1928年荷兰的范德波耳(Van der Pol)为描述LC 回路的电子管振荡器建立了著名的vanderPol方 程.它在自激振荡理论中有着重要的意义,一直作 为数学物理方程中的一个基本方程.这是一个具 有可变非线性阻尼的微分方程,代表了一类极为 典型的非线性问题.和其他非线性微分方程在数 学上无法精确求解一样,人们一直在努力寻找求 解这类方程近似解析解的方法,并乐于用Van der Pol方程来检验求解方法的有效性.
常微分方程的数值解
常微分方程中只有一些典型方程能求出初等解(用初 等函数表示的解) 。另外，有些初值问题虽然有初等解， 但由于形式太复杂不便于应用。因此，有必要探讨常微 分方程初值问题的数值解法。 以下主要介绍一阶常微分方程初值问题几种经典数值 解方法：欧拉法及改进的欧拉法。
其它方法：龙格-库塔法、阿达姆斯方法； 一阶微分方程组与高阶方程初值问题的数值解法； 二阶常微分方程值问题的差分方法等。
k 1) (k ) （k 1） 对于已给的精确度 ， 当满足 yi( y 时， 取 y y 1 i 1 i 1 i 1 ，
然后继续下一步 y i 2 的计算。
此即改进的欧拉法
在函数y( x)可微的条件下，由微分 中值定理可知，存在 0 1, 使得
y ( xn 1 ) y ( xn ) y ' ( xn h) h
解 输入命令：dsolve('Du=1+u^2','t')
d2y 0 2 dx
结
果：u = tg(t-c)
例 2 求微分方程的特解.
d 2 y dy 2 4 29 y 0 dx dx y (0) 0, y ' (0) 15
解 输入命令: y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
•线性多步法有四阶阿达姆斯外插公式和内插公式。
ODE 指令列表

MATLAB 用於求解常微分方程式的指令:
指令 ode45 方法 Explicit Runge-Kutta (4, 5) pair of Dormand-Prince Explicit Runge-Kutta (2, 3) pair of Bogacki and Shampine Variable order Adams-Bashforth-Moulton PECE solver Numerical differentiation formulas (NDFS) 应用ODE类别 Nonstiff ODE
3、使用泰勒公式 以此方法为基础，有龙格-库塔法、线性多步法等方 法。 4、数值公式的精度 当一个数值公式的截断误差可表示为O（hk+1）时 （k为正整数，h为步长），称它是一个k阶公式。 k越大，则数值公式的精度越高。 •欧拉法是一阶公式，改进的欧拉法是二阶公式。
•龙格-库塔法有二阶公式和四阶公式。
ode23
ode113 ode15s
Nonstiff ODE
Nonstiff ODE Stiff ODE
ode23s
ode23t ode23tb
Modified Rosenbrock formula of order 2
Trapezoidal rule with a “free” interpolant Implicit Runge-Kutta formula with a backward differentiation formula of order two
分格式， hi xi xi 1 称为步长，实用中常取定步长。
建立数值解法的一些途径
设 x i 1 xi h, i 0,1,2, n 1, 可用以下离散化方法求 解微分方程： y' f(x,y) y(x0 ) y 0
1、用差商代替导数 若步长h较小，则有
实际应用时，与欧拉公式结合使用：
0) y i( y i hf ( xi , y i ) 1 h ( k 1) (k ) y y [ f ( x , y ) f ( x , y i 1 i i i i 1 i 1 )] k 0,1,2, 2
得积分的步长(Step Sizes)变得很小,以便降低积分 误差至可容忍范围以內,会导致计算时间过长 专门对付Stiff系统的指令，例如ode15s、ode23s、 ode23t及ode23tb
提示

使用 Simulink 來求解常微分方程式
Simulink是和MATLAB共同使用的一套软件 可使用拖拉的方式來建立动力系统 可直接产生C语言代码或进行动画演示 功能非常强大

17世纪：初等解法 18世纪：初等解法和无穷级数方法
19世纪：解的存在性、奇点理论、定性 理论、稳定性理论
定理
设函数 f x, y 在区域 D : a x b, y
上连续，且在区域D内满足李普希兹(Lipschitz)条件，即存在 正数L，使得对于R内任意两点 x, y1 与 x, y 2 ,恒有
y ( xi 1 ) y ( xi )