微分方程的数值解法matlab(四阶龙格—库塔法)

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)];
f
2
(t
,
Y
)
Y0
y1(0)
y2
(0)
(2.2) (2.3)
例：著名的Van der Pol方程
y ( y2 1) y y 0
令 y1 y, y2 y
Y
y1 y2
降为一阶
Y
y1
y2
( y12
y2 1) y2
y1
初始条件
Y0
y1(0)
y2
(0)
y10
y20
tn t0
for i = 1 : N
tn1 tn h
子程序计算 Ki
Next i
yn1
yn
h 6
(K1
2K2
2K3
K4)
tn tn1, yn yn1
输出结果
Plot
End
三. Runge-Kutta 法解Van der Pol 方程的 Matlab 程序结构
主程序：RK_vanderpol.m 子程序：RK_sub.m（函数文件）
Y (t) AY (t) 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
四. Matlab对应命令：ode23,ode45
说明：
t0：初始时刻；tN：终点时刻
y0：初值； tol：计算精度
调用格式： [t, y]=ode23 (‘函数文件名’, t0, tN, y0, tol) [t, y]=ode45 (‘函数文件名’, t0, tN, y0, tol)
默认精度： ode23——1e-3
初值问题的数值 解法分为两大类
单步法-Runge-Kutta 方法 多步法-Admas方法
计算 y(tn1) 的近似值 yn1 时只用到 tn , yn ，是自开始
方法
Runge-Kutta法是常微分方程的一种经典解法 MATLAB 对应命令：ode45
四阶Runge-Kutta公式
yn1
yn
各种solver 解算指令的特点
解法指令 解题类 型
特点
ode45 非刚性 采用4、5阶Runge－Kutta法
适合场合 大多数场合的首选算法
ode23 非刚性 采用Adams算法
较低精度（10－3）场合
ode113 非刚性
ode23t ode15s
适度刚 性
刚性
ode23s 刚性 ode23tb 刚性
低精度时，比ode15s有 效；或存在质量矩阵时
低精度时，比ode15s有 效；或存在质量矩阵时
二. 四 阶 Runge-Kutta 法
dy f (t, y) a t b dt
y(t0 ) y0
对 I=[a,b]作分割
a t0 t1 tN1 b
步长
hi hi1 hi ,i 0,1,, N 1
多步法；采用Adams算法；高 低精度均可（10－3～10－6）
采用梯形法则算法
ode45计算时间太长时 取代ode45
适度刚性
多步法；采用2阶Rosenbrock 当ode45失败时使用；
算式，精度中等
或存在质量矩阵时
一步法；采用2阶Rosenbrock 算式，低精度
采用梯形法则－反向数值微分 两阶段算法，低精度
例题2：三自由 度质量弹簧系统
P0sin(wt) x1
k1
k2
m1
x2
k3 m2
x3
k4 m3
m1x1 k1x1 k2 (x1 x2 ) P0 sin(t) 0
m2x2 k2 (x1 x2 ) k3(x2 x3) 0
m3x3 k3(x2 x3) k4x3 0
矩阵表示
m1 0
0
m2
0 0
0 x1 k1 k2
0
x2
k2
m3 x3 0
k2 k2 k3
k3
0 x1 P0 sin(t)
k3
x2
0
k3 k4 x3 0
其中：
x1
X
x2
x3
m1 0 0
M
0
m2
0
0 0 m3
k1 k2
K
k2
0
k2 k2 k3
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) 时刻的外部激励力
K2 = ZCX_sub (t0 + h/2, Y0 + hK1/2, P )
微分方程的数值解法
四阶龙格—库塔法 （The Fourth-Order Runge－Kutta
Method）
常微分方程（Ordinary differential equations, ODE)
初值问题---给出初始值 边值问题---给出边界条件
与初值常微分方程解算有关的指令
ode23
ode45
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 输出数据或图形
4. 使编写好的ODE函数文件和初值 Y0 供微分 方程解算指令（solver）调用
Solver解算指令的使用格式
说明：
t0：初始时刻；tN：终点时刻 Y0：初值； tol：计算精度
[t, Y]=solver (‘ODE函数文件名’, t0, tN, Y0, tol);
ode45
输出宗量形式
y1(t0 )
Y
y1
(t1
)
y1
(t2
)
y2 (t0 )
y2
(t1
)
y2
(t2
)
Y (:,1) y1(t) Y (:,2) y2 (t)
例题1：著名的Van der Pol方程
y ( y2 1) y y 0
解法1：采用 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))
X (t) 令： Y (t) X (t),
Y (t) R2n,
X (0)
Y
(0)
X
(0)
Y (t) AY (t) P (t)
(2)
MX (t) CX (t) KX (t) P(t), X (t) Rn
X (0) X 0 ,
X (0) X 0
即：
其中：
X X
(t (t
) )
h 6
(k1
2k2
2k3
k4 )
k1 f (tn , yn )
1
h
k2 f (tn 2 h, yn 2 k1)
1
h
k3 f (tn 2 h, yn 2 k2 )
k4 f (tn h, yn hk3 )
四 阶 Runge-Kutta 法计算流程图
开始
h 初迭始代条次件数： ：t0N；y0 积分步长：
k3
0
k3
k3 k4
P0 sin(t)
P(t)
0
0
动力学方程：
例如：
MX (t) KX (t) P(t)
已知参数：m1=m2=m3=1, k1=2, k2=2, K3=1, K4=2, P0=1, 3
要求：采用四阶龙格－库塔法编程计算三个质量的响应时程.计算时 间 0 ~ 50
解法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
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
0 -0.1 -0.2 -0.3 -0.4
0
0
-0.1
-0.2
-0.3
5 10 15 20 25 30 35 40 45 50
-0.4 0
5 10 15 20 25 30 35 40 45 50
结果完全一致
MATLAB程序 （1）4阶RK方法： （2）采用ode45： m_chap2_ex2_1.m，m_chap2_ex2_1_sub.m
解析解：
x1 x2
(t) (t)
1 0.0882
P0 k
x3 (t )
1
3
sin(wt) 2.63
0
3
P0 k
sin(wt) 0.21
2
2
2
P0 k
sin(wt)
第一个质量的位移响应时程
4阶龙格－库塔法的结果
0.4
ode45 的结果
0.4
0.3
0.3
0.2
0.2
0.1
0.1
3. 根据式(2.2)编写计算导数的M函数文件ODE文件
把t,Y作为输入宗量，把 Y 作为输出宗量
%M function file name: dYdt.m function Yd = f (t, Y) Yd = f (t,Y) 的展开式
例Van der Pol方程
%M function file name: dYdt.m function Yd = f (t, Y) Yd=zeros(size(Y)); Yd (1) Y (2); Yd (2) (Y (1).^2 1)*Y (2) Y (1);
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）
要求：写出程序体，打印所绘图形，图形标题用个 人的名字。
Duffing 方程
y 0.3y y y3 0.28 cos(1.2t) y(0) 0.01, y(0) 0.0
五. 动力学系统的求解
1. 动力学方程
MX (t) CX (t) KX (t) P(t),
X (0) X 0 ,
X (0) X 0
X (t) Rn
(1)
2. 二阶方程转成一阶方程
ode113
ode23t
ode15s
ode23s
ode23tb
一.解ODE的基本机理:
1. 列出微分方程
y f ( y, y,t)
(2.1)
y(0) y0 , y(0) y0
令 y1 y, y2 y
初始条件
2. 把高阶方程转换成一阶微分方程组
Y
y1
y2
f
(t,Y )
f1(t,Y )
ቤተ መጻሕፍቲ ባይዱ
子程序：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)];
0nn M 1K
I nn M 1C
X X
(t (t
) )
0n1 M 1P
(t
)
A
0nn M 1K
Inn M 1C
R2n2n
,
P
(t)
0n1 M 1P(t
)
R
2n1
Y (t) AY (t) P (t)
(2)
Y (t) AY (t) P (t)
3. Matlab 程序(主程序：ZCX)