MATLAB微分方程解析解数值解图形解解法教程

Matlab常微分方程求解问题分类
边值问题：
初值问题：
• 定解附加条件在自变量 的一端
• 一般形式为： y' f (x, y)
y(a)
y0
• 初值问题的数值解法一 般采用步进法，如 Runge-Kutta法
➢ 在自变量两端均给定附加 条件
y' f (x, y)
➢ 一般形式：y(a)y1, y(b)y2
1.根据常微分方程要求的求解精度与速度要求
求解初值问题：
y
'
y
2x y
y ( 0 ) 1
(0x1)
比较ode45和ode23的求解精度和速度
ode45和ode23的比较-1
function xODE clear all clc
format long
y0 = 1; [x1,y1] = ode45(@f,[0,1],y0); [x2,y2] = ode23(@f,[0,1],y0); plot(x1,y1,'k-',x2,y2,'b--') xlabel('x') ylabel('y')
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);
% Mass balances dCdt = [rA; rB; rC; rD; rE];
三个串联的CSTR等温反应器（例4-3）
function IsothermCSTRs clear all clc CA0 = 1.8; % kmol/m^3 CA10 = 0.4; % kmol/m^3 CA20 = 0.2; % kmol/m^3 CA30 = 0.1; % kmol/m^3 k = 0.5; % 1/min tau = 2; stoptime = 2.9; % min [t,y] = ode45(@Equations,[0 stoptime],[CA10 CA20 CA30],[],k,CA0,tau); disp(' Results:') disp(' t CA1 CA2 CA3') disp([t,y]) plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-') legend('CA_1','CA_2','CA_3') xlabel('Time (min)') ylabel('Concentration') % -----------------------------------------------------------------function dydt = Equations(t,y,k,CA0,tau) CA1 = y(1); CA2 = y(2); CA3 = y(3); dCA1dt = (CA0-CA1)/tau - k*CA1; dCA2dt = (CA1-CA2)/tau - k*CA2; dCA3dt = (CA2-CA3)/tau - k*CA3; dydt = [dCA1dt; dCA2dt; dCA3dt];

xn1 xn
f (t, y(t))dt .
(2.4)
右端积分用左矩形公式hf(xn,y(xn))近似，再以yn代替 y(xn)，yn+1代替y(xn+1)也得到欧拉公式(2.1)，局部截
断误差也是(2.3).
如果右端积分用右矩形公式hf(xn+1,y(xn+1))近似， 则得到另一个公式
yn1 yn hf ( xn1, yn1 ),
基于上述几何解释，我们从初始点P0(x0, y0)出发, 先依方向场在该点的方向推进到x=x1上一点P1，然后 再从P1点依方向场在该点的方向推进到 x=x2 上一点 P2 , 循环前进做出一条折线P0 P1 P2.
上页 下页
一般地，设已做出该折线的顶点Pn,过Pn(xn, yn)依
方向场的方向再推进到Pn+1(xn+1, yn+1)，显然两个顶
y1 y0 hf ( x0, y0 ),
y2 y1 hf (x1, y1),
上页 下页
例1 用欧拉公式求解初值问题
y
y
2x y
(0 x 1),
y(0) 1.
(2.2)
解 取步长h=0.1，欧拉公式的具体形式为
yn1
yn
h( yn
2xn ) yn
其中xn=nh=0.1n (n=0,1,,10), 已知y0 =1, 由此式可得
)],
于是有
yn1
y(k1) n1
hL 2
yn1
y(k) n1
,
使得
hL 1,
2
则当k→∞时有yn(k)1 yn1 , 这说明迭代过程(2.8)是收敛的.
上页 下页
9.2.3 改进的欧拉公式

参数说明
[T,Y] = solver(odefun,tspan,y0)
odefun 为显式常微分方程，可以用命令 inline 定义，或 在函数文件中定义，然后通过函数句柄调用。
dy 2 2 y 2 x 2x 求初值问题 的数值解，求解范 例： dx 围为 [0,0.5] y( 0 ) 1
dsolve的输出个数只能为一个 或 与方程个数相等。
只有很少一部分微分方程（组）能求出解析解。 大部分微分方程（组）只能利用数值方法求数值解。
Matlab函数数值求解
[T,Y] = solver(odefun,tspan,y0)
其中 y0 为初值条件，tspan为求解区间；Matlab在数值求解 时自动对求解区间进行分割，T (列向量) 中返回的是分割点 的值(自变量)，Y (数组) 中返回的是这些分割点上的近似解， 其列数等于因变量的个数。
数学实验
实验四
求微分方程的解
问题背景和实验目的
自牛顿发明微积分以来，微分方程在描述事物运 动规律上已发挥了重要的作用。实际应用问题通过 数学建模所得到的方程，绝大多数是微分方程。 由于实际应用的需要，人们必须求解微分方程。 然而能够求得解析解的微分方程十分有限，绝大多 数微分方程需要利用数值方法来近似求解。 本实验主要研究如何用 Matlab 来计算微分方程 （组）的数值解，并重点介绍一个求解微分方程的 基本数值解法－－Euler折线法。
Runge-Kutta 方法
Euler 法与 R-K法误差比较
Matlab 解初值问题
用 Maltab自带函数 解初值问题 求解析解：dsolve 求数值解：
ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb

科学计算与MATLAB
2012
第七讲 常微分方程数值解法
内容提要
引言 欧拉近似方法 龙格-库塔(R-K)方法 MATLAB的常微分方程函数 小结
1、引言
物理学所研究的各种物质运动中，有许多物质运动的 过程是用常微分方程来描述的。
例如,质点的加速运动，简谐振动等。
F m dv dt
d2x 2x2 0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
作如下近似
yn y(xn )
得：
yn1 yn f yn , xn h
2.1.4 欧拉法误差
利用泰勒级数得：
y xn1 y(xn h)
y(xn )
hy(xn )
1 2
h2
y(xn )
y(xn )
x2 …. xn ….
y(x0) y(x1) y(x2) …. y(xn) ….
y0
y1 y2 …. yn ….
在xn节点上，微分方程可以写为
y(xn1) y(xn ) f y(xn ) , xn h
作如下近似：
yn y(tn )
则得到欧拉解法递推公式的一般形式：
yn1 yn f ( yn , xn ) h
hf
y(xn ),
xn
1 2
h2 y(xn )
作如下近似
yn y(xn )
yn1 yn f yn , xn h
局部截断误差
y
y0
dy dx
x0 x x0
y0 f ( y0 , x0 ) x x0
此切线与x=x1交点纵坐标为：
y1 y0 f ( y0 , x0 ) x1 x0

差,输出参数,事件等，可缺省。 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

d2x 2x2 0
dt 2
简朴问题可以求得解析解，多数实际问题靠数值求解 。
第4页
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立旳节点 xn 上旳近似值 yn，用以近似y(xn)
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
第17页
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
yn y(xn )
第5页
2、欧拉近似办法
2.1 简朴欧拉(L.Euler, 1707-1783)办法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一种函数值由前一种函数值得到。核 心是构造递推公式。
y0 y1 y2 yn
第6页
i 1,2,...
第36页
没有一种算法可以有效地解决所有旳 ODE 问题，因此 MATLAB 提供了多种ODE函数。
函数 ODE类
特点
阐明
型
ode45
非刚性 单步法；4，5 阶 R-K 措施；合计 大部分场合旳首选措施
截断误差为 (△x)3
ode23
非刚性 单步法；2，3 阶 R-K 措施；合计 使用于精度较低旳情形

例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 =3*exp(-2*x)*sin(5*x) 作图命令： 作图命令：ezplot(y,[1.0,4])
解析法假设导弹在t时刻的位置为pxtyt乙舰位于由于导弹头始终对准乙舰故此时直线pq就是导弹的轨迹曲线弧op处的切线即有又根据题意弧op的长度为aq3数学建模实例由12消去t整理得模型
《高等数学》 高等数学》
—上机教学（三） 上机教学（ 微分方程求解
上机目的
1、学会用 Matlab 求简单微分方程的解析解 、 求简单微分方程的解析解. 2、学会用 Matlab 求微分方程的数值解 、 求微分方程的数值解.
例3
求微分方程组的通解. 求微分方程组的通解 dx dt = 2 x − 3 y + 3 z dy = 4 x − 5 y + 3z dt dz = 4 x − 4 y &olve('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)
假设导弹在 t 时刻的位置为 P(x(t), y(t))，乙舰位于 Q(1, v0 t ) . y(t))

3月15日作业: 1.Van der Pol 方程的两种解法:1)采 用ode45命令 2)Runge-Kutta方法 2.Duffing 方程的求解（Runge-Kutta方法，计算步长h=0.005，计
算时间t0=0.0,tN=100）
要求：写出程序体，打印所绘图形，图形标题用个 人的名字。
Duffing 方程Fra bibliotek子程序：RK_sub.m function ydot = vdpol (t, y) ydot=zeros(size(y)); ydot(1) = y(2); ydot(2) = -y(2)*(y(1)^2-1)-y(1); 或写为：
ydot = [y(1) ;-y(2)*(y(1)^2-1)-y(1)];
K4 = ZCX_sub (t0 + h, y0 + hK3, P) Y1 = y0 + (h/6) (K1 + 2K2 + 2K3 + K4) t1, Y1 (输出 t1, y1) next i 输出数据或图形
(t ) AY (t ) P (t ) Y
Matlab 程序(子程序：ZCX_sub.m)
多步法-Admas方法
计算
的近似值 y(tn1 )
时只用到 n1
y
，是自开始方法 n n
t ,y

Runge-Kutta法是常微分方程的一种经典解法 MATLAB 对应命令：ode45
四阶Runge-Kutta公式
h yn 1 yn (k1 2k 2 2k3 k 4 ) 6 k1 f (t n , yn ) 1 h k 2 f (t n h, yn k1 ) 2 2 1 h k3 f (t n h, yn k 2 ) 2 2 k 4 f (t n h, yn hk3 )

6
但能求解析解的常微分方程是有限的，大多数的常
微分方程是给不出解析解的.
例
2 2 y x y
这个一阶微分方程就不能用初等函数及其积
分来表达它的解。
例
y y y (0) 1
x x y e , e 的值仍需插值方法来计算. 的解
7
8
事实上，从实际问题当中抽象出来的微分方程，通
15
同理，在x= xn 处，用差商代替导数： y( xn1 ) y( xn ) y( xn1 ) y( xn ) y( xn ) xn1 xn h 由 得 若记
y( xn ) f ( xn , yn )
y( xn1 ) y( xn ) hf ( xn , yn ) y( xn ) yn , y( xn1 ) yn1
26
Two Body PБайду номын сангаасoblem
(t ) u (t ) / r (t ) 3 u (t ) v (t ) / r (t ) 3 v
r (t ) u (t ) 2 v(t ) 2
u (t ) v(t ) y (t ) u (t ) v ( t )
则上式可记为
yn1 yn hf ( xn , yn )
此 即 为 求 解 初 值 问 题 的 Euler 方 法 ， 又 称 显 式 16 Euler方法。
例: 用Euler方法求解常微分方程初值问题
y 2 y 2 y x y (0) 0. (0 x3 )
值积分、数值微分、泰勒展开等离散化方法，对初 值问题
y f ( x, y ) y ( x0 ) y 0
y