如何使用MATLAB求解微分方程(组)

Where ?
工程控制 航空航天
ODE
医学生理 金融分析
Where ?
算法开发 数据分析
数值计算 MAT LAB
数据可视化
When ?
当对问题进行建模后，有常微分方程 需要求解时。
在生物建模中，经常需要求解常微分 方程。如药物动力学的房室模型的建模 仿真。
How ?
数 值 解
数值解？刚性方程？
果有初始条件，则求出特解。 用字符串表示常微分方程，自变量缺省时为t,导数用
D表示微分。y的2阶导数用D2y表示，依此类推。
如何调用？
[T,Y,TE,YE,IE]=solver('odefun',tspan,y0,options)
其中solver为ode23、ode45、ode113、ode15s、ode23s、
通过列写方程组，建立odefun.m文 件即方程组文件 function dC=odefun(t,C)
dC=[0.1*C(3)+37.5-2.52*C(1); 1.68*C(1)-0.01*C(2); 0.004*C(2)-0.1*C(3)];
end 程序关键部分列写如下： [t,C]=ode15s('fun1',[0,30],[20.3,341 0.5,136.4])
数值解是采用如有限元方 法、 数值逼近方法、插值方 法等计算方法得到的解。只 能利用数值计算的结果， 而 不能随意给出自变量并求出 计算值。
当无法藉由微积分技巧求 得解析解时，这时便只能利 用数值分析的方式来求得其 数值解了。实际情况下，常 微分方程往往只能求解出其
数值解。
在数学中，刚性方程是指一 个微分方程，其数值分析的解 只有在时间间隔很小时才会稳 定，只要时间间隔略大，其解 就会不稳定。
E.g.3结果
C1(t)/(ug/L)
20.3 20.295
20.29 20.285
0
3410.6
体 内 无 机 碘 浓 度 C1(t)
5
10
15
20
25
30
t/d
甲 状 腺 中 有 机 碘 浓 度 C2(t)
C2(t)/(ug/L)
xdot(2)=1/C*x(1);
e方n式d 二：无额外参数
function dC=odefun(t,C)
dC=[ 0.1*C(3)+37.5-2.52*C(1); %直接书写列矩阵
1.68*C(1)-0.01*C(2);
0.004*C(2)-0.1*C(3)];
end
Examples
E.g.1 求下列微分方程组的通解 x'=2x-3y+3z y'=4x-5y+3z z'=4x-4y+2z
差,输出参数,事件等，可缺省。
使用ODE?时如何编 写微分方程 ？
方式一：带额外参数，使用时需对参数进行赋值
function odefun(t,x,flag,R,L,C) %用flag说明R、L、C为变
量
xdot=zeros(2,1);
%表明其为列向量
xdot(1)=-R/L*x(1)-1/L*x(2)+1/L*f(t);
y(1)=x(2);
y1
y2
y(2)= -t*x(1)+exp(t)*x(2)+3*sin(2*t);
end
1000
求解程序关键步骤
[t,y]=ode45('odefun',[3.9 4.0],[2 8])
y
500
0
3.9
3.92
3.94
3.96
3.98
4
4.02
t
Examples
E.g.3 求解人体中不同形式的碘浓度的三房室模型。已知碘在三房室之间的 转换速率为：k21=0.84/d，k01=1.68/d，k32=0.01/d，k13=0.08/d， k03=0.02/d，f10=150g/d，f30=0g/d。初始条件为：x1(0)=81.2g， x2(0)=6821g，x3(0)=682g，v1=4L, v2=2L, v3=5L。仿真时间为30d。
目前很难去精确地去定义哪 些微分方程是刚性方程，但是 大体的想法是：这个方程的解
包含有快速变化的部分。
如何调用？
y=dsolve('e1,e2,...','c1,c2,...','v')
其中'e1,e2,...'为微分方程或微分方程组； 'c1,c2,...',是初始条件或边界条件； 'v'是独立变量，默认的独立变量是't'； y 返回解析解。如果没有初始条件，则求出通解，如
Topic: 如何使用MATLAB求 解常微分方程（组）
TMU_BME_2013
a.What ?
微分方程指描述未知函数的导数与自变 量之间的关系的方程。未知函数是一元函 数的微分方程称作常微分方程。未知函数 是多元函数的微分方程称作偏微分方程。
MATLAB（matrix&laboratory）意为矩 阵工厂（矩阵实验室）.MATLAB是美国 MathWorks公司出品的商业数学软件，提 供高级技术计算语言和交互式环境，主要 包括MATLAB和Simulink两大部分。
Examples
E.g.2 求解y''== - t*y + e^t*y' +3sin(2t)。 已知3.9<t<4.0，y'(0)=2，y''(0)=8
函数文件odefun.m的建立
程序执行结果
function y=odefun(t,x) y=zeros(2,1); %列向量
y' '=-t*y + et*y' +3sin2t 1500
输入命令： [x,y,z]=dsolve(‘Dx=2*x-3*y+3*z’,’Dy=4*x-5*y+3*z’,
‘Dz=4*x4*y+2*z’,’t’) x=simple(x); y=simple(y); z=simple(z); %化简结果
运行结果为： x=(c1-c2+c3+c2e -3t-c3e-3t)e2t y= - c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z=(-c1 e-4t+c2e-4t+c1-c2+c3)e2t
ode23t、ode23tb 函数；
odefun 是函数句柄；
tspan 微分定义区间；
y0
为初值行矩阵；
T
值是t序列（为列向量）；
Y
值是微分方程的解YHale Waihona Puke Baidu各点t的值（为列向量）；
TE
表示事件发生时间，可缺省；
YE
表示事件解决时间，可缺省；
IE
表示事件消失时间，可缺省；
options 是求解参数设置,可以用odeset在计算前设定误