如何使用MATLAB求解微分方程(组) PPT
第五讲 Matlab求解微分方程
第五讲 Matlab求解微分方程教学目的：学会用MATLAB求简单微分方程的解析解、数值解，加深对微分方程概念和应用的理解；针对一些具体的问题，如追击问题，掌握利用软件求解微分方程的过程；了解微分方程模型解决问题思维方法及技巧；体会微分方程建摸的艺术性．教学重点：利用机理分析建模微分方程模型，掌握追击问题的建模方法，掌握利用MATLAB求解数值解．教学难点：利用机理分析建模微分方程模型，通过举例，结合图形以及与恰当的假设突破教学难点．1微分方程相关函数（命令）及简介因为没有一种算法可以有效地解决所有的ODE 问题，为此，Matlab 提供了多种求解器Solver，对于不同的ODE 问题，采用不同的Solver．阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度．2 求解微分方程的一些例子2.1 几个可以直接用 Matlab 求微分方程精确解的例子：例1：求解微分方程22x xe xy dxdy-=+，并加以验证． 求解本问题的Matlab 程序为：syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2 diff(y,x)+2 *x*y-x*exp(-x^2) %line3 simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4说明：(1) 行line1是用命令定义x,y 为符号变量．这里可以不写，但为确保正确性，建议写上；(2) 行line2是用命令求出的微分方程的解：1/2*exp(-x^2)*x^2+exp(-x^2)*C1(3) 行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出：-x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 表明)(x y y =的确是微分方程的解．例2：求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解，并画出解函数的图形．求解本问题的 Matlab 程序为： syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即xe e y x+=，此函数的图形如图 1：图1 y 关于x 的函数图象2.2 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例３：求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222y xx y dx dy 的数值解，求解范围为区间[0, 0.5]．fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); x'; y'; plot(x,y,'o-') >> x' ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000>> y' ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图 2．图2 y 关于x 的函数图像3 常微分在实际中的应用 3.1 导弹追踪问题1、符号说明，w ，乙舰的速率恒为v 0；设时刻t 乙舰的坐标为((),())X t Y t ，导弹的坐标为((),())x t y t ；当零时刻，((0),(0))(1,0)X Y =，((0),(0))(0,0)x y =．建立微分方程模型．202,01(0)0,'(0)0v d yk x dx w y y ⎧⎪⎪= =<<⎨⎪⎪==⎩由微分方程模型解得11(1)(1)11(),12(1)2(1)2(1)2(1)k k x x y x k k k k k +-+--=-+-≠+-+-++代入题设的数据1/5k =，得到导弹的运行轨迹为4655555(1)(1)81224y x x =--+-+当1=x 时245=y ，即当乙舰航行到点)245,1(处时被导弹击中. 被击中时间为:00245v v y t ==. 若v 0=1, 则在t =0.21处被击中. 利用MALAB 作图如图３．clear, x=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'*')图３导弹运行轨迹（解析法） 图４两种方法对比的导弹运行轨迹2、数值方法求解．设导弹速率恒为w ，则得到参数方程为⎪⎪⎩⎪⎪⎨⎧--+-=--+-=)()()()()()(2222y Y y Y x X wdt dy x X y Y x X w dt dx因乙舰以速度0v 沿直线1x =运动，设01v =，5,1,w X Y t ===，因此导弹运动轨迹的参数方程为：⎪⎪⎪⎩⎪⎪⎪⎨⎧==--+-=--+-=0)0(,0)0()()()1(5)1()()1(52222y x y t y t x dtdyx y t x dt dxMATLAB 求解数值解程序如下，结果见图４ t0=0,tf=0.21;[t,y]=ode45('eq2',[t0 tf],[0 0]); X=1;Y=0:0.001:0.21;plot(X,Y,'-') plot(y(:,1),y(:,2),'*'),hold onx=0:0.01:1; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24; plot(x,y,'r')很明显数值计算的方法比较简单而且适用． 3.2 蚂蚁追击问题（1）建立平面直角坐标系．以正多边形的中心为原点, 设正多边形的一个顶点为起始点, 连接此点和原点作x 轴. 根据x 轴作出相应的y 轴; 选取足够小的t ∆进行采样．（2）在每一时刻t ，计算每只蚂蚁在下一时刻t t +∆时的坐标．不妨设甲追逐对象是乙，在时间t 时，甲的坐标为A 11(,)x y ，乙的坐标为B 22(,)x y ．甲在t t +∆时在'A 点(如图1), 其坐标为11(cos ,sin )x v t y v t θθ+∆+∆，其中2221212121cos ,sin ,()()x x y yd x x y y d dθθ--===-+-. 同理，依次计算下一只蚂蚁在t t +∆时的坐标．通过间隔t ∆进行采样，得到新一轮各个蚂蚁在一个新的正多边形位置坐标．（4）重复2）步，直到d 充分小为止．（5）连接每只蚂蚁在各时刻的位置，就得到所求的轨迹．用MALAB 求解并作图，函数zhuJi(x,y)在附录一定义，如图６ t=[1:8]; s=7*exp(t.*2*pi/length(t)*i); x=real(s); y=imag(s);zhuJi(x,y)图6当蚂蚁为7只时的图形习题1. 求微分方程0sin 2')1(2=-+-x xy y x 的通解．2. 求微分方程x e y y y x sin 5'2''=+-的通解．3. 求微分方程组'A 11,(,)A x y 22,(,)B x y dθ图1 在采样时间内，相连蚂蚁追击⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dtdx在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =的图形．4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t ∈．利用画图来比较两种求解器之间的差异．4 参考文献:[1] Mastering Matlab 6, D. Hanselman, B. Littlefield, 清华大学出版,2002[2] 赵静等编．数学建模与数学实验（第3版）．北京：高等教育出版社．2008. [3] 姜启源编. 数学模型（第二版）．北京：高等教育出版社.1993.[4] 石勇国. 蚂蚁追击问题与等角螺线. 宜宾学院学报. 2008，(6): 23-25. [5] 张伟年，杜正东，徐冰．常微分方程．北京：高等教育出版社．2006．5 附录附录一:zhuji(x,y)的M 文件 function zhuji(x,y) clf v=1; dt=0.05;x(length(x)+1)=x(1);y(length(y)+1)=y(1); plot(x,y,'*') holdfor it=1:1000for i=1:length(x)-1d=sqrt((x(i)-x(i+1))^2+(y(i)-y(i+1))^2); x(i)=x(i)+v*dt*(x(i+1)-x(i))/d; y(i)=y(i)+v*dt*(y(i+1)-y(i))/d; endplot(x,y,'.') hold onx(length(x))=x(1); y(length(y))=y(1); end。
第六章用MATLAB求解微分方程及微分方程组
运行结果：u = tan(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') 运行结果为 : y =3e-2xsin（5x）
不同求解器Solver的特点
求解器 Solver ode23s
OD E类 型
刚 性
特点
说明
刚 ode23tb 性
一步法；2阶 当精度较低时， Rosebrock算法； 计算时间比 低精度 ode15s短 当精度较低时， 梯形算法；低精 计算时间比 度 ode15s短
d 2x dx 2 1000(1 x 2 ) x 0 例 4 dt dt x(0) 2; x' (0) 0 解: 令 y1=x，y2=y1’
y Me k1t 即
模型2：快速饮酒后，体液中酒精含量的变化率
dx k2 x k1Mek1t dt x(0) 0
用Matlab求解模型2：
syms x y k1 k2 M t x=dsolve('Dx+k2*x=k1*M*exp(-k1*t)','x(0)=0','t') 运行结果： M*k1/(-k1+k2)*exp(-k2*t+t*(-k1+k2))-exp(-k2*t)*M*k1/(-k1+k2) 即：
（1）
2. 由于弹头始终对准乙舰，故导弹的速度平行于乙舰与导弹头位置的差向量，
用MATLAB求解微分方程及微分方程组
例 3 求微分方程组的通解. dx dt 2 x 3 y 3 z dy 4 x 5 y 3z dt dz 4 x 4 y 2 z dt
任取k1、k2的一组初始值：k0=[2,1];
输入命令： k=lsqcurvefit('curvefun1',k0,t,c) 运行结果为： k =[ 1.3240 作图表示求解结果： t1=0:0.1:18; f=curvefun1(k,t1); plot(t,c,'ko',t1,f,'r-')
90 80 70 60 50 40 30 20 10 0
0.2573]
0
2
46Leabharlann 81012
14
16
18
模型2：慢速饮酒时，体液中酒精含量的变化率
dx k2 x a dt x(0) 0
其中
M a T
M为饮酒的总量，T为饮酒的时间
则有；
a x (1 e k 2 t ) k2
5 5 ) 处时被导弹击中. 当 x 1时 y ，即当乙舰航行到点 (1, 24 24 y 5 被击中时间为: t . 若 v0=1, 则在 t=0.21 处被击中. v0 24v0
轨迹图如下
例： 饮酒模型
模型1：快速饮酒后，胃中酒精含量的变化率
dy k1 y dt y (0) M
5 5 ) 处时被导弹击中. 当 x 1时 y ，即当乙舰航行到点 (1, 24 24 y 5 被击中时间为 : t . 若 v0=1, 则在 t=0.21 处被击中. v0 24v0
matlab常微分方程 ppt课件
odeplot 随计算过程，显示解向
量
属性名 OutputSel
Refine
2020/4/10
参数设置
取值
含义
若不使用缺省设置，则
有效值： 正整数向 量 缺省值：[]
OutputFcn所表现的是那 些正整数指定的解向量 中的分量的曲线或数据 。若为缺省值时，则缺 省地按上面情形进行操
作
有效值： 正整数k>1 缺省值：k
y1 (0 ) y 0
y0
y
2
(
0
)
y1
yn
(
0
)
yn
• （3）根据（1）与（2）的结果，编写能计 算
• 导数的M-函数文件odefile。
• （4）将文件odefile与初始条件传递给求解
• 器Solver中的一个，运行后就可得到ODE
• 的、在指定时间区间上的解列向量y（其中包 2含020/4y/10及不同阶的导数）。
参数设置
属性名 NormControl Events
取值
含义
有效值： on、off 缺省值：
off
为‘on’时，控制解向量 范数的相对误差，使每 步计算中，满足： norm(e)<=max(RelTol*n orm(y),AbsTol)
有效值： 为‘on’时，返回相应的 on、off 事件记录
2020/4/10
例3
x' x 2
x(0)
1
创建函数function2, 保存在function2.m中
function f=function2(t,x)
f=-x.^2;
在命令窗口中执行 >> [t,x]=ode45('function2',[0,1],1); >> plot(t,x,'-',t,x,'o'); >> xlabel('time t0=0,tt=1'); >> ylabel('x values x(0)=1');
如何使用MATLAB求解微分方程(组)ppt课件
差,输出参数,事件等，可缺省。 9
使用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);
2
Where ?
工程控制
ODE
医学生理
航空航天
金融分析
3
Where ?
算法开发 数据分析
数值计算 MAT LAB
数据可视化
4
When ?
当对问题进行建模后，有常微分方程 需要求解时。
在生物建模中，经常需要求解常微分 方程。如药物动力学的房室模型的建模 仿真。
5
方法 ode23
ode45
数 ode113
当无法藉由微积分技巧求 得解析解时，这时便只能利 用数值分析的方式来求得其 数值解了。实际情况下，常 微分方程往往只能求解出其
数值解。
在数学中，刚性方程是指一 个微分方程，其数值分析的解 只有在时间间隔很小时才会稳 定，只要时间间隔略大，其解 就会不稳定。
目前很难去精确地去定义哪 些微分方程是刚性方程，但是 大体的想法是：这个方程的解
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
matlab求解偏微分方程组
matlab求解偏微分方程组摘要：I.引言- 介绍偏微分方程组及其应用领域- 阐述使用Matlab 求解偏微分方程组的重要性II.Matlab 的基本操作与设置- 安装与配置Matlab- 熟悉Matlab 的基本语法与操作III.偏微分方程组的类型及表示方法- 分类介绍偏微分方程组：椭圆型、抛物型和双曲型- 常用偏微分方程组的表示方法IV.使用Matlab 求解偏微分方程组- 使用Matlab 内置函数求解偏微分方程组- 利用有限差分法、有限元法等数值方法求解偏微分方程组- 分析求解结果及误差V.偏微分方程组在实际问题中的应用- 举例说明偏微分方程组在物理、工程等领域的应用- 展示Matlab 在解决实际问题中的优势VI.结论- 总结Matlab 求解偏微分方程组的方法与技巧- 展望偏微分方程组在数值计算领域的未来发展正文：偏微分方程组广泛应用于物理、工程、生物学等领域，是描述现实世界中许多复杂现象的基本工具。
Matlab 作为一款功能强大的数学软件，为求解偏微分方程组提供了便利。
本文将介绍如何使用Matlab 求解偏微分方程组，并探讨其在实际问题中的应用。
首先，我们需要安装与配置Matlab 软件。
在安装过程中，请确保正确配置Matlab 的环境变量，以便在后续使用过程中能够顺利调用相关函数。
熟悉Matlab 的基本语法与操作也是十分必要的，这将有助于我们更好地使用Matlab 解决实际问题。
偏微分方程组根据其特征可分为椭圆型、抛物型和双曲型。
在实际应用中，我们通常会遇到各种不同类型的偏微分方程组，因此熟练掌握不同类型方程组的表示方法十分重要。
使用Matlab 求解偏微分方程组的方法有很多。
Matlab 内置了许多求解偏微分方程组的函数，例如`fsolve`、`dsolve`等。
此外，我们还可以利用有限差分法、有限元法等数值方法求解偏微分方程组。
在实际求解过程中，我们需要根据问题的特点选择合适的求解方法，并对求解结果进行误差分析，以确保结果的准确性。
Matlab微分方程的应用1PPT课件
2021/3/12
13
变分法（calculus of variations）：是处理函数的函数 的数学领域，和处理数的函数的普通微积分相对。譬如， 这样的泛函可以通过未知函数的积分和它的导数来构造。 变分法最终寻求的是极值函数：它们使得泛函取得极大 或极小值。
变分法的关键定理是欧拉－拉格朗日方程。它对应于 泛函的临界点。在寻找函数的极大和极小值时，在一个 解附近的微小变化的分析给出一阶的一个近似。它不能 分辨是找到了最大值或者最小值（或者都不是）。 变分原理 （variational principle）：把一个物理学问 题（或其他学科的问题）用变分法化为求泛函极值（或 驻值）的问题，后者就称为该物理问题 （或其他学科的 问题）的变分原理。
y(x)bxxln x(x1)
a d a 2a
2021/3/12
9
5 总结与讨论
y
方法 利用微分不等式建模；
有时只需求近似解。 b
模型讨论
1）视点移动时升起
o a dd
x
曲线如何求得？
2）怎样减少地面的坡度？调整参数、相邻排错位。
3）衡量经济的指标？ 座位尽量多、升起曲线占据的空间尽量少等。
2021/3/12
2021/3/12
观众厅地面设计
在影视厅或报告厅，经常会为前边观众遮挡住 自己的视线而苦恼。显然，场内的观众都在朝台上 看，如果场内地面不做成前低后高的坡度模式，那 么前边观众必然会遮挡后面观众的视线。试建立数 学模型设计良好的报告厅地面坡度曲线。
2021/3/12
1
问题的假设
1) 观众厅地面的纵剖面图一致，只需求中轴线 上地面的起伏曲线即可。
2) 同一排的座位在同一等高线上。
用MATLAB求解微分方程
1. 微分方程的解析解
求微分方程（组）的解析解命令:
dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
结 果：u = tan(t-c)
解 输入命令：dsolve('Du=1+u^2','t')
STEP2
STEP1
解 输入命令: y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
导弹追踪问题
设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中？ 解法一（解析法）
由(1),(2)消去t整理得模型:
解法二(数值解)
结 果 为：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
2、取t0=0，tf=12，输入命令： [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
3、结果如图
图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.
微分方程的matlab求解数学建模PPT学习教案
第12页/共44页
返回
使用改进欧拉公式的例
步长 h= 0.1 的数值解比较
x
精表确解 向前欧拉 向后欧拉 改进欧拉
0
1
1
1
1
0.1
1.0048
1
1.0091
1.005
0.2
1.0187
1.01
1.0264
1.019
0.3
1.0408 1.029 1.0513 1.0412
0.4
1.0703 1.0561 1.0830 1.0708
y = 3/4*exp(-x)+1/4*exp(3*x)
第23页/共44页
例③ 非常系数的二阶微分方程
x ''( t ) ( 1 x 2 ( t ) x '( ) t ) x ( t ) 0 , x ( 0 ) 3 , x '( 0 ) 0
x=dsolve('D2x-(1-x^2)*Dx+x=0', ' x(0)=3,Dx(0)=0')
第22页/共44页
例② 常系数的二阶微分方程
输入
y '' 2 y ' 3 y 0 , y ( 0 ) 1 ,y '( 0 ) 0
：y=dsolve('D2y-2*Dy-3*y=0','x')
y=dsolve('D2y-2*Dy3*y=0','y(0)=1,Dy(0)=0','x')
结果：
y = C1*exp(-x)+C2*exp(3*x)
第4页/共44页
1、欧拉方法
2) 向后欧拉公式
matlab教学第二章 微分方程PPT课件
作线性组合得到平均斜率，由此得到更高阶的精
度，这就是龙格-库塔方法的基本思路。
在MATLAB软件中含有数值求解的系统函数，其
实现原理就是龙格－库塔方法。
其中参数k >0，m=18。
2020/11/12
8
2.2.2 利用平衡与增长式
许多研究对象在数量上常常表现出某种不变的特性， 如封闭区域内的能量、货币量等。利用变量间的平衡 与增长特性，可分析和建立有关变量间的相互关系。
此类建模方法的关键是分析并正确描述基本模型的右 端，使平衡式成立。
例2.2 战斗模型：两方军队交战，希望为这场战斗建 立一个数学模型，应用这个模型达到如下目的：
“速率”、“增长”、“衰变”、“边际的”等常涉及到
导数。
我们熟悉的速度公式：dy v 就是一个简单的一阶微分方
程。
dt
微分方程是指含有导数或微分的等式。 一般形式： F(x,y,y,,y(n))0
或y： (n) f(x,y,y,,y(n1)).
常用的建立微分方程的方法有：运用已知物理定律；利用 平衡与增长式；运用微元法；应用分析法。
y2(k+1)=(y2(k)+h*x1(k+1)+h)/(1+h);
y3(k+1)=(y3(k)+(h/2)*(-y3(k)+x1(k)+x1(k+1)+2))/(1+h/2);
end
x=0:0.1:1;
y=x+exp(-x);
x1=x1(1:11),y=y(1:11),y1=y1(1:11),y2=y2(1:11),y3=y3( 1:11),
202092040捕食者的存在使食饵的增长率降低假设x1降低的程度与捕食者数量x2成正比即食饵对捕食者的数量x2起到增长的作用其程度与食饵数量x1成正比即dtdxdtdxdtdxdtdx生产理论把企业仅仅抽象为一个生产函数一种投入产出关系一个追求利润最大化的黑匣子它没有讨论企业内部是如何配置资源的企业是如何组织生产的企业和市场的关系如何各自的边界在哪里
偏微分方程的matlab解法PPT课件
求解抛物型方程的例子
考虑一个带有矩形孔的金属板上的热传导问题。板的左边
保持在100 °C，板的右边热量从板向环境空气定常流动，
t t 其他边及内孔边界保持绝缘。初始
°C ，于是概括为如下定解问题;
是板的温度为0 0
d u u0 , t
u 1 0 0 ,在 左 边 界 上
u 1， 在 右 边 界 上 n u = 0， 其 他 边 界 上 n
ppt课件完整
3
先确定方程大类
ppt课件完整
4
Draw Mode
画图模式，先将处理的区域画出来，二 维，方形，圆形，支持多边形，可以手 动更改坐标，旋转rotate
例如，对于细杆导热，虽然是一维问题， 可以将宽度y虚拟出来，对应于y的边界 条件和初始条件按照题意制定
ppt课件完整
5
Boundary Mode
ppt课件完整
20
ppt课件完整
21
第七步：单击Plot菜单中Parameter选项，打开Plot Selection对话框，选中Color,Height(3D plot)和 Show mesh三项.再单击Polt按钮，显示三维图形解， 如图22.5所示.
ppt课件完整
22
第八步：若要画等值线图和矢量场图，单击plot菜单 中parameter 选项，在plot selection对话框中选中 contour 和arrow两选项。然后单击plot按钮，可显示 解的等值线图和矢量场图，如图2.6所示。
2u (2u t 2 x2
u t to 0
区域的边界顶点坐标为(-0.5,-0.8), (0.5,-0.8), (-0.5,0.8), (-0.05,0.4) ,(0.05,-0.4),
微分方程的Matlab求解ppt课件
y1' y2 y3
y2 ' y1 y3 y3' 0.51y1 y2
y1(0) 0, y2 (0) 1, y3(0) 1
解 1、建立m-文件如下：
function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2);
s.z x=simple(s.x) %简化结果
y=simple(s.y)
z=simple(s.z)
结 果 为：
x =-(-C1-C2*exp(-3*t)+C2-C3+C3*exp(-3*t))*exp(2*t)
y =(-C1*exp(-4*t)+C1+C2*exp(-4*t)+C2*exp(-3*t)-C2+C3-C3*exp(-3*t))*exp(2*t)
x(0) 0, y(0) 0
返回
2. 模型求解
(1) w=20时,建立m-文件如下: function dy=eq3(t,y) dy=zeros(2,1); dy(1)=20*(10+20*cos(t)-y(1))/sqrt ((10+20*cos(t)-y(1))^2+(20+15*sin(t)-y(2))^2); dy(2)=20*(20+15*sin(t)-y(2))/sqrt ((10+20*cos(t)-y(1))^2+(20+15*sin(t)-y(2))^2);
注意: 1、在解n个未知函数的方程组时，x0和x均为n维向量，
m-文件中的待解方程组应以x的分量形式写成.
2、使用Matlab软件求数值解时，高阶微分方程必须 变换成等价的一阶微分方程组.
Matlab求解微分方程(组)及偏微分方程(组)
第四讲 Matlab 求解微分方程（组）理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为：X=dsolve(‘eqn1’,’eqn2’,…)函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解.注意，系统缺省的自变量为t2.函数dsolve 求解的是常微分方程的精确解法，也称为常微分方程的符号解.但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为solver ，其一般格式为：[T,Y]=solver(odefun,tspan,y0)说明：(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一.(2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解.(3)如果要获得微分方程问题在其他指定时间点012,,,,f t t t t 上的解，则令tspan 012[,,,]f t t t t =(要求是单调的).(4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，Matlab 提供了多种求解器solver ，对于不同的ODE 问题，采用不同的solver.表1 Matlab中文本文件读写函数说明：ode23、ode45是极其常用的用来求解非刚性的标准形式的一阶微分方程（组）的初值问题的解的Matlab常用程序，其中：ode23采用龙格-库塔2阶算法，用3阶公式作误差估计来调节步长，具有低等的精度.ode45则采用龙格-库塔4阶算法，用5阶公式作误差估计来调节步长，具有中等的精度.3．在matlab命令窗口、程序或函数中创建局部函数时，可用内联函数inline，inline函数形式相当于编写M函数文件，但不需编写M-文件就可以描述出某种数学关系.调用inline函数，只能由一个matlab表达式组成，并且只能返回一个变量，不允许[u,v]这种向量形式.因而，任何要求逻辑运算或乘法运算以求得最终结果的场合，都不能应用inline函数，inline函数的一般形式为：FunctionName=inline(‘函数内容’, ‘所有自变量列表’)例如：（求解F(x)=x^2*cos(a*x)-b ,a,b是标量；x是向量）在命令窗口输入:Fofx=inline(‘x .^2*cos(a*x)-b’ , ‘x’,’a’,’b’); g= Fofx([pi/3 pi/3.5],4,1) 系统输出为：g=-1.5483 -1.7259注意：由于使用内联对象函数inline 不需要另外建立m 文件，所有使用比较方便，另外在使用ode45函数的时候，定义函数往往需要编辑一个m 文件来单独定义，这样不便于管理文件，这里可以使用inline 来定义函数. 二．实例介绍1.几个可以直接用Matlab 求微分方程精确解的实例 例1 求解微分方程2'2x y xy xe -+=程序：syms x y; y=dsolve(‘Dy+2*x*y=x*exp(-x^2)’,’x ’)例 2 求微分方程'0x xy y e +-=在初始条件(1)2y e =下的特解并画出解函数的图形.程序：syms x y; y=dsolve(‘x*Dy+y-exp(1)=0’,’y(1)=2*exp(1)’,’x ’);ezplot(y)例 3 求解微分方程组530tdx x y e dtdy x y dt⎧++=⎪⎪⎨⎪--=⎪⎩在初始条件00|1,|0t t x y ====下的特解并画出解函数的图形.程序：syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') simple(x); simple(y)ezplot(x,y,[0,1.3]);axis auto2.用ode23、ode45等求解非刚性标准形式的一阶微分方程（组）的初值问题的数值解（近似解）例 4 求解微分方程初值问题2222(0)1dy y x xdx y ⎧=-++⎪⎨⎪=⎩的数值解，求解范围为区间[0,0.5].程序：fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1); plot(x,y,'o-')例 5 求解微分方程22'2(1)0,(0)1,(0)0d y dyy y y y dt dtμ--+===的解，并画出解的图形.分析：这是一个二阶非线性方程，我们可以通过变换，将二阶方程化为一阶方程组求解.令12,,7dyx y x dtμ===，则 121221212,(0)17(1),(0)0dx x x dtdx x x x x dt⎧==⎪⎪⎨⎪=--=⎪⎩ 编写M-文件vdp.m function fy=vdp(t,x)fy=[x(2);7*(1-x(1)^2)*x(2)-x(1)]; end在Matlab 命令窗口编写程序 y0=[1;0][t,x]=ode45(@vdp,[0,40],y0);或[t,x]=ode45('vdp',[0,40],y0); y=x(:,1);dy=x(:,2); plot(t,y,t,dy)练习与思考：M-文件vdp.m 改写成inline 函数程序？ 3.用Euler 折线法求解Euler 折线法求解的基本思想是将微分方程初值问题00(,)()dyf x y dxy x y ⎧=⎪⎨⎪=⎩ 化成一个代数(差分)方程，主要步骤是用差商()()y x h y x h +-替代微商dydx，于是00()()(,())()k k k k y x h y x f x y x h y y x +-⎧=⎪⎨⎪=⎩记1,(),k k k k x x h y y x +=+=从而1(),k k y y x h +=+于是0011(),,0,1,2,,1(,).k k k k k k y y x x x h k n y y hf x y ++=⎧⎪=+=-⎨⎪=+⎩例 6 用Euler 折线法求解微分方程初值问题22(0)1dyx y dxy y ⎧=+⎪⎨⎪=⎩的数值解（步长h 取0.4），求解范围为区间[0,2].分析：本问题的差分方程为00110,1,0.4,0,1,2,,1(,).k k k k k k x y h x x h k n y y hf x y ++===⎧⎪=+=-⎨⎪=+⎩程序：>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1; >> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y});%subs ，替换函数 x=x+h; szj=[szj;x,y]; end >>szj>> plot(szj(:,1),szj(:,2))说明：替换函数subs 例如：输入subs(a+b,a,4) 意思就是把a 用4替换掉，返回 4+b ，也可以替换多个变量，例如：subs(cos(a)+sin(b),{a,b},[sym('alpha'),2])分别用字符alpha 替换a 和2替换b ，返回 cos(alpha)+sin(2)特别说明：本问题可进一步利用四阶Runge-Kutta 法求解，Euler 折线法实际上就是一阶Runge-Kutta 法，Runge-Kutta 法的迭代公式为001112341213243(),,(22),6(,),0,1,2,,1(,),22(,),22(,).k k k k k k k k k k k k y y x x x h h y y L L L L L f x y k n h h L f x y L h h L f x y L L f x h y hL ++=⎧⎪=+⎪⎪=++++⎪⎪=⎪=-⎨⎪=++⎪⎪⎪=++⎪⎪=++⎩相应的Matlab 程序为：>> clear >> f=sym('y+2*x/y^2'); >> a=0; >> b=2; >> h=0.4; >> n=(b-a)/h+1; >> x=0; >> y=1;>> szj=[x,y];%数值解 >> for i=1:n-1l1=subs(f, {'x','y'},{x,y});替换函数 l2=subs(f, {'x','y'},{x+h/2,y+l1*h/2}); l3=subs(f, {'x','y'},{x+h/2,y+l2*h/2}); l4=subs(f, {'x','y'},{x+h,y+l3*h}); y=y+h*(l1+2*l2+2*l3+l4)/6; x=x+h; szj=[szj;x,y]; end>>szj>> plot(szj(:,1),szj(:,2))练习与思考：(1)ode45求解问题并比较差异. (2)利用Matlab 求微分方程(4)(3)''20y y y -+=的解.(3)求解微分方程''2',2(1)0,030,(0)1,(0)0y y y y x y y --+=≤≤==的特解. (4)利用Matlab 求微分方程初值问题2''''00(1)2,|1,|3x x x y xy y y ==+===的解. 提醒：尽可能多的考虑解法 三．微分方程转换为一阶显式微分方程组Matlab 微分方程解算器只能求解标准形式的一阶显式微分方程（组）问题，因此在使用ODE 解算器之前，我们需要做的第一步，也是最重要的一步就是借助状态变量将微分方程（组）化成Matlab 可接受的标准形式.当然，如果ODEs 由一个或多个高阶微分方程给出，则我们应先将它变换成一阶显式常微分方程组.下面我们以两个高阶微分方程组构成的ODEs 为例介绍如何将它变换成一个一阶显式微分方程组.Step 1 将微分方程的最高阶变量移到等式左边，其它移到右边，并按阶次从低到高排列.形式为：()'''(1)'''(1)()'''(1)'''(1)(,,,,,,,,,,)(,,,,,,,,,,)m m n n m n x f t x x x x y y y y y g t x x x x y y y y ----⎧=⎨=⎩Step 2 为每一阶微分式选择状态变量，最高阶除外'''(1)123'''(1)123,,,,,,,,,m m n m m m m n x x x x x x x x x y x y x y x y--++++========注意：ODEs 中所有是因变量的最高阶次之和就是需要的状态变量的个数，最高阶的微分式不需要给它状态变量.Step 3 根据选用的状态变量，写出所有状态变量的一阶微分表达式''''122334123''12123,,,,(,,,,,),,(,,,,,)m m n m m m nm n x x x x x x x f t x x x x xx xg t x x x x +++++======练习与思考：(1)求解微分方程组**'''3312*'''3312()()22x x x y x r r y y y x y r r μμμμμμ⎧+-=+--⎪⎪⎨⎪=+--⎪⎩其中2r =1r =*1,μμ=-1/82.45,μ=(0) 1.2,x =(0)0,y ='(0)0,x ='(0) 1.049355751y =-(2)求解隐式微分方程组''''''''''''2235x y x y x y x y xy y ⎧+=⎨++-=⎩ 提示：使用符号计算函数solve 求'''',x y ，然后利用求解微分方程的方法 四．偏微分方程解法Matlab 提供了两种方法解决PDE 问题，一是使用pdepe 函数，它可以求解一般的PDEs,具有较大的通用性，但只支持命令形式调用；二是使用PDE 工具箱，可以求解特殊PDE 问题，PDEtoll 有较大的局限性，比如只能求解二阶PDE 问题，并且不能解决片微分方程组，但是它提供了GUI 界面，从复杂的编程中解脱出来，同时还可以通过File —>Save As 直接生成M 代码.1.一般偏微分方程（组）的求解(1)Matlab 提供的pdepe 函数，可以直接求解一般偏微分方程（组），它的调用格式为：sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)@pdefun 是PDE 的问题描述函数，它必须换成标准形式：(,,)[(,,,)](,,,)m m u u u uc x t x x f x t u s x t u x t x x x-∂∂∂∂∂=+∂∂∂∂∂ 这样，PDE 就可以编写入口函数：[c,f,s]=pdefun(x,t,u,du)，m,x,t 对应于式中相关参数，du 是u 的一阶导数，由给定的输入变量可表示出c,f,s 这三个函数.@pdebc 是PDE 的边界条件描述函数，它必须化为形式：(,,)(,,).*(,,,)0up x t u q x t u f x t u x∂==∂ 于是边值条件可以编写函数描述为：[pa,qa,pb,qb]=pdebc(x,t,u,du),其中a 表示下边界，b 表示上边界.@pdeic 是PDE 的初值条件，必须化为形式：00(,)u x t u =，故可以使用函数描述为：u0=pdeic(x)sol 是一个三维数组，sol(:,:,i)表示i u 的解，换句话说，k u 对应x(i)和t(j)时的解为sol(i,j,k)，通过sol ，我们可以使用pdeval 函数直接计算某个点的函数值.(2)实例说明 求解偏微分2111222221220.024()0.17()u u F u u t xu u F u u tx ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩ 其中， 5.7311.46()x x F x e e -=-且满足初始条件12(,0)1,(,0)0u x u x ==及边界条件1(0,)0,u t x ∂=∂221(0,)0,(1,)1,(1,)0uu t u t t x∂===∂ 解：(1)对照给出的偏微分方程和pdepe 函数求解的标准形式，原方程改写为111221220.024()1.*()10.17u u F u u x u F u u u t x x ∂⎡⎤⎢⎥--⎡⎤⎡⎤⎡⎤∂∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥-∂∂∂⎣⎦⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦可见1121220.024()10,,,()10.17u F u u x m c f s F u u u x ∂⎡⎤⎢⎥--⎡⎤⎡⎤∂====⎢⎥⎢⎥⎢⎥-∂⎣⎦⎣⎦⎢⎥⎢⎥∂⎣⎦ %目标PDE 函数function [c,f,s]=pdefun(x,t,u,du) c=[1;1];f=[0.024*du(1);0.17*du(2)]; temp=u(1)-u(2);s=[-1;1].*(exp(5.73*temp)-exp(-11.46*temp)) end(2)边界条件改写为：下边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦上边界1110.*000u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦%边界条件函数function [pa,qa,pb,qb]=pdebc(xa,ua,xb,ub,t) pa=[0;ua(2)]; qa=[1;0]; pb=[ub(1)-1;0]; qb=[0;1]; end(3)初值条件改写为：1210u u ⎡⎤⎡⎤=⎢⎥⎢⎥⎣⎦⎣⎦%初值条件函数 function u0=pdeic(x) u0=[1;0]; end(4)编写主调函数 clc x=0:0.05:1; t=0:0.05:2; m=0;sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t); subplot(2,1,1) surf(x,t,sol(:,:,1)) subplot(2,1,2) surf(x,t,sol(:,:,2))练习与思考： This example illustrates the straightforward formulation, computation, and plotting of the solution of a single PDE.2()u u t x xπ∂∂∂=∂∂∂ This equation holds on an interval 01x ≤≤ for times 0t ≥. The PDE satisfies the initial condition (,0)sin u x x π= and boundary conditions(0,)0;(1,)0t uu t e t xπ-∂=+=∂ 2.PDEtool 求解偏微分方程(1)PDEtool （GUI ）求解偏微分方程的一般步骤在Matlab 命令窗口输入pdetool ，回车，PDE 工具箱的图形用户界面(GUI)系统就启动了.从定义一个偏微分方程问题到完成解偏微分方程的定解，整个过程大致可以分为六个阶段Step 1 “Draw 模式”绘制平面有界区域Ω，通过公式把Matlab 系统提供的实体模型：矩形、圆、椭圆和多边形，组合起来，生成需要的平面区域.Step 2 “Boundary 模式”定义边界，声明不同边界段的边界条件.Step 3 “PDE 模式”定义偏微分方程，确定方程类型和方程系数c,a,f,d ，根据具体情况，还可以在不同子区域声明不同系数.Step 4 “Mesh 模式”网格化区域Ω，可以控制自动生成网格的参数，对生成的网格进行多次细化，使网格分割更细更合理.Step 5 “Solve 模式”解偏微分方程，对于椭圆型方程可以激活并控制非线性自适应解题器来处理非线性方程；对于抛物线型方程和双曲型方程，设置初始边界条件后可以求出给定时刻t 的解；对于特征值问题，可以求出给定区间上的特征值.求解完成后，可以返回到Step 4，对网格进一步细化，进行再次求解.Step 6 “View 模式”计算结果的可视化，可以通过设置系统提供的对话框，显示所求的解的表面图、网格图、等高线图和箭头梯形图.对于抛物线型和双曲线型问题的解还可以进行动画演示.(2)实例说明用法求解一个正方形区域上的特征值问题：12|0u u u u λ∂Ω⎧-∆-=⎪⎨⎪=⎩ 正方形区域为：11,1 1.x x -≤≤-≤≤(1)使用PDE 工具箱打开GUI 求解方程(2)进入Draw 模式，绘制一个矩形，然后双击矩形，在弹出的对话框中设置Left=-1,Bottom=-1,Width=2,Height=2,确认并关闭对话框(3)进入Boundary 模式，边界条件采用Dirichlet 条件的默认值(4)进入PDE 模式，单击工具栏PDE 按钮，在弹出的对话框中方程类型选择Eigenmodes,参数设置c=1,a=-1/2,d=1,确认后关闭对话框(5)单击工具栏的 按钮，对正方形区域进行初始网格剖分，然后再对网格进一步细化剖分一次(6)点开solve菜单，单击Parameters选项，在弹出的对话框中设置特征值区域为[-20,20](7)单击Plot菜单的Parameters项，在弹出的对话框中选中Color、Height(3-D plot)和show mesh项，然后单击Done确认(8)单击工具栏的“=”按钮，开始求解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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。
通过列写方程组，建立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])
Topic: 如何使用MATLAB求 解常微分方程（组）
TMU_BME_2013
a.What ?
微分方程指描述未知函数的导数与自变 量之间的关系的方程。未知函数是一元函 数的微分方程称作常微分方程。未知函数 是多元函数的微分方程称作偏微分方程。
MATLAB（matrix&laboratory）意为矩 阵工厂（矩阵实验室）.MATLAB是美国 MathWorks公司出品的商业数学软件，提 供高级技术计算语言和交互式环境，主要 包括MATLAB和Simulink两大部分。
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
果有初始条件，则求出特解。 用字符串表示常微分方程，自变量缺省时为t,导数用
D表示微分。y的2阶导数用D2y表示，依此类推。
如何调用？
[T,Y,TE,YE,IE]=solver('odefun',tspan,y0,options)
其中solver为ode23、ode45、ode113、ode15s、ode23s、
输入命令： [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
差,输出参数,事件等，可缺省。
使用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=dsolve('e1,e2,...','c1,c2,...','v')
其中'e1,e2,...'为微分方程或微分方程组； 'c1,c2,...',是初始条件或边界条件； 'v'是独立变量，默认的独立变量是't'； y 返回解析解。如果没有初始条件，则求出通解，如
数值解是采用如有限元方 法、 数值逼近方法、插值方 法等计算方法得到的解。只 能利用数值计算的结果， 而 不能随意给出自变量并求出 计算值。
当无法藉由微积分技巧求 得解析解时，这时便只能利 用数值分析的方式来求得其 数值解了。实际情况下，常 微分方程往往只能求解出其
数值解。
在数学中，刚性方程是指一 个微分方程，其数值分析的解 只有在时间间隔很小时才会稳 定，只要时间间隔略大，其解 就会不稳定。
ode23t、ode23tb 函数；
odefun 是函数句柄；
tspan 微分定义区间；
y0
为初值行矩阵；
T
值是t序列（为点t的值（为列向量）；
TE
表示事件发生时间，可缺省；
YE
表示事件解决时间，可缺省；
IE
表示事件消失时间，可缺省；
options 是求解参数设置,可以用odeset在计算前设定误
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
Where ?
工程控制 航空航天
ODE
医学生理 金融分析
Where ?
算法开发 数据分析
数值计算 MAT LAB
数据可视化
When ?
当对问题进行建模后，有常微分方程 需要求解时。
在生物建模中，经常需要求解常微分 方程。如药物动力学的房室模型的建模 仿真。
How ?
数 值 解
数值解？刚性方程？
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)