#matlab第二章__常微分方程的数值解法

例题2：三自由 度质量弹簧系统
P0sin(wt) x1
k1
k2
m1
x2
k3 m2
x3
k4 m3
m 1 x 1 k 1 x 1 k 2 ( x 1 x 2 ) P 0 si t ) n 0
m 2 x 2 k 2 ( x 1 x 2 ) k 3 ( x 2 x 3 ) 0
0

0 0 m 3 x 3 0 k3 k3k4 x3 0
其中：
x1
X


x
2

x 3
m1 0 0
初值问题的数值 解法分为两大类
单步法-Runge-Kutta 方法 多步法-Admas方法
计算 y(tn1) 的近似值 y n 1 时只用到 t n , y n ，是自开始
方法
Runge-Kutta法是常微分方程的一种经典解法 MATLAB 对应命令：ode45
四阶Runge-Kutta公式
A M 0 n n 1 K M In n 1 C R 2 n 2 n , P (t) M 0 1 n P 1 (t) R 2 n 1
Y (t)A(t) Y P (t)
(2)
Y (t)A(t) Y P (t)
t1
)

y
2
(t
2
)

Y(:1),y1(t) Y(:2,)y2(t)
例题1：著名的Van der Pol方程 解法1：采用
y (y 2 1 )y y 0 ODE命令
% 主程序 （程序名：VanderPol _ex1.m） t0 = 0; tN = 20; tol = 1e-6; Y0 = [0.25; 0.0]; [t, Y]=ode45 (‘dYdt’, t0, tN, Y0, tol); subplot (121), plot (t, Y) subplot (122), plot (Y( :, 1), Y( :, 2))
Van der Pol方程
% 子程序 （程序名： dYdt.m ） function Ydot = dYdt (t, Y) Ydot=[Y(2);-Y(2)*(Y(1)^2-1)-Y(1)];
或写为
function Ydot = dYdt (t, Y) Ydot=zeros(size(Y)); Ydot(1)=Y(2); Ydot(2)=-Y(2)*(Y(1).^2-1)-Y(1)];
子程序：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)];
Duffing 方程
y 0 . 3 y y y 3 0 . 2 c 8 1 o . 2 t )s( y ( 0 ) 0 . 0 ,y ( 0 1 ) 0 . 0
五. 动力学系统的求解
1. 动力学方程
M X (t) C X (t) K(tX ) P (t),X (t) R n
主程序：RK_vanderpol.m 子程序：RK_sub.m（函数文件）
解法2：采用Runge_Kutta法编程计算
主程序：RK_vanderpol.m t0=0; tN=20; y0=[0.25; 0]; h=0.001; t = t0 : h : tN; N = length (t); j = 1; for i = 1 : N
低精度时，比ode15s有 效；或存在质量矩阵时
低精度时，比ode15s有 效；或存在质量矩阵时
二. 四 阶 Runge-Kutta 法
dy f (t, y) a t b dt y(t0) y0
对 I=[a,b]作分割
a t0 t1 tN 1 b
步长
h i h i 1 h i,i 0 ,1 , ,N 1
X (0 ) X 0 ,
X (0 ) X 0
(1)
2. 二阶方程转成一阶方程
X (t) 令： Y (t) X (t) ,
Y (t) R 2n,
Y (0) X X ((0 0))
Y (t)A(t) Y P (t)
(2)
M X (t) C X (t) K(tX ) P (t),X (t) R n
t1 = t0 + h; K1 = RK_sub(t0, y0); K2 = RK_sub(t0 + h/2, y0 + h*K1/2); K3 = RK_sub(t0 + h/2, y0 + h*K2/2); K4 = RK_sub(t0 + h, y0 + h*K3); y1 = y0 + (h/6)*(K1 + 2*K2 + 2*K3 + K4); yy1(j)=y1(1); yy2(j)=y1(2); t0=t1; y0=y1; j=j+1; end subplot (121), plot (t, yy1, t, yy2); grid subplot (122), plot (yy2, yy1); grid
微分方程的数值解法
四阶龙格—库塔法 （TheFourth-OrderRunge－
KuttaMethod）
常微分方程（Ordinary differential equations, ODE)
初值问题---给出初始值 边值问题---给出边界条件
与初值常微分方程解算有关的指令
ode23
ode45
四. Matlab对应命令：ode23,ode45
说明：
t0：初始时刻；tN：终点时刻
y0：初值； tol：计算精度
调用格式： [t, y]=ode23 (‘函数文件名’, t0, tN, y0, tol) [t, y]=ode45 (‘函数文件名’, t0, tN, y0, tol)
默认精度： ode23——1e-3
Solver解算指令的使用格式
说明：
t0：初始时刻；tN：终点时刻 Y0：初值； tol：计算精度
[t, Y]=solver (‘ODE函数文件名’, t0, tN, Y0, tol);
ode45
输出宗量形式
y1 (t0 )
Y


y1
(t1
)

y
1
(t
2
)
y2 (t0 )
y
2
(
y n1

yn

h 6
(k1
2k2

2k3

k4)
k1 f (tn , yn )
k2

f (tn
1 2
h,
yn

h 2
k1)
k3

f (tn
1 2
h,
yn

h 2
k
2
)
k 4 f (tn h , y n hk 3 )
四 阶 Runge-Kutta 法计算流程图
开始
h 初始条件：t
m 3 x 3 k 3 (x 2 x 3 ) k 4 x 3 0
矩阵表示
m 1 0 0 x 1 k1k2 k2 0x 1 P 0sin t)(
0 m 2 0 x 2 k2
k2 k3
k3 x2
3. Matlab 程序(主程序：ZCX)
t0;Y0;h;N;P0,w; %输入初始值、步长、迭代次数、初始激励力；
for i = 1 : N
t1 = t0 + h
P=[P0*sin(w*t0);0.0;0.0]
%输入t0时刻的外部激励力
K1 = ZCX_sub (t0, Y0, P )
P=
%输入 (t0+h/2) 时刻的外部激励力
例Van der Pol方程
%M function file name: dYdt.m function Yd = f (t, Y) Yd=zeros(size(Y));
Y(1 d)Y(2); Y(d2)(Y(1).2^1)*Y(2)Y(1);
4. 使编写好的ODE函数文件和初值 Y 0 供微分 方程解算指令（solver）调用
初始条件
Y0

y1(0) y2(0)

y10 y20
3. 根据式(2.2)编写计算导数的M函数文件ODE文件
把t,Y作为输入宗量，把 Y 作为输出宗量
%M function file name: dYdt.m function Yd = f (t, Y) Yd = f (t,Y) 的展开式
X (0 ) X 0 ,
X (0 ) X 0
即：
其中：
X X ( (tt) ) M 0 n n 1 K M In n 1 C X X ( (tt) ) M 0 1 n P 1 (t)
ode45——1e-6
3月15日作业: 1.Van der Pol 方程的两种解法:1) 采用ode45命令 2)Runge-Kutta方法 2.Duffing 方程的求解（Runge-Kutta方法，计算步长
h=0.005，计算时间t0=0.0,tN=100）
要求：写出程序体，打印所绘图形，图形标题用个 人的名字。
K2 = ZCX_sub (t0 + h/2, Y0 + hK1/2, P )
K3 = ZCX_sub (t0 + h/2, Y0 + hK2/2, P )
P=
%输入 (t0+h) 时刻的外部激励力
K4 = ZCX_sub (t0 + h, y0 + hK3, P) Y1 = y0 + (h/6) (K1 + 2K2 + 2K3 + K4) t1, Y1 (输出 t1, y1) next i 输出数据或图形
Y (t)A(t) Y P (t)
Matlab 程序(子程序：ZCX_sub.m)
function ydot = f (t, Y,P) M=, K=, C= %输入结构参数 P1=[zeros(3,1);inv（M）*P]； A=[zeros(0,0), eye(n,n) ; -M-1K, -M-1C] ydot =AY+P1
各种solver 解算指令的特点
解法指令 解题类 型
特点
ode45 非刚性 采用4、5阶Runge－Kutta法
适合场合 大多数场合的首选算法
ode23 非刚性 采用Adams算法
较低精度（10－3）场合
wenku.baidu.com
ode113 非刚性
ode23t ode15s
适度刚 性
刚性
ode23s 刚性 ode23tb 刚性
多步法；采用Adams算法；高 低精度均可（10－3～10－6）
采用梯形法则算法
ode45计算时间太长时 取代ode45
适度刚性
多步法；采用2阶Rosenbrock 当ode45失败时使用；
算式，精度中等
或存在质量矩阵时
一步法；采用2阶Rosenbrock 算式，低精度
采用梯形法则－反向数值微分 两阶段算法，低精度
Y0

y1(0)

y
2
(0)
(2.2) (2.3)
例：著名的Van der Pol方程
y (y 2 1 )y y 0
令
y1y,y2y
Y


y y
1 2

降为一阶 Yyy1 2(y12y12 )y2y1
ode113
ode23t
ode15s
ode23s
ode23tb
一.解ODE的基本机理:
1. 列出微分方程
y f(y,y ,t)
(2.1)
y (0 ) y 0 ,y (0 ) y 0 初始条件
令 y1y,y2y
2. 把高阶方程转换成一阶微分方程组
Yyy1 2f(t,Y)ff1 2((tt,,Y Y))
迭代次数：
；y
0
N
0
积分步长：
tn t0
for i = 1 : N
tn1 tn h
子程序计算 K i
Next i yn 1yn6 h(K 12K 22K 3K 4)
tntn1,ynyn1
输出结果
Plot
End
三. Runge-Kutta 法解Van der Pol 方程的 Matlab 程序结构