g
h x(t0 )

h2 2!
(
ftk

h2 2
x(2) (t0
' fxk ' fk )
)
f 'tk f 'xk 等各阶导数不易计算，用下式中 ki的线性组合代替
r
xk1 xk h aiki
线性组合
i 1
r为精度阶次，ai为待定系数，由精度确定；ki用下 式表示
i 1
ki f (tk b1h, xk hb2 k j ) , i 2, 3L , r
2

A(xk

h 2
K1) Bu(tk

h) 2

K
3

A(xk

h 2
K2)

Bu(tk

h) 2
K4 A(xk hK3) Bu(tk h)
yk+1 Cxk+1
RK法的特点：
1 需要存储的数据少，占用的存储空间少； 2 只需知道初值，即可启动递推公式进行计
算，可自启动； 3 容易实现变步长运算。 4 每积分一步需要计算多次右函数，计算量
求出微分方程在这些时刻的近似值 x0 x1 x2 L xN
x(t0

h)

x(t0
)

h
g
x(t0 )

h2 2!
x(2) (t0
)
L

hk k!
x(k)
(t0 )
1. 欧拉法
矩形面积
取前两项近似： xk1 hf (tk , xk ) xk
取k=0,1,2,…N,从t0开始，逐点递推求解t1时的y1， t2时的 y2…,直至tn时的yn，称之为欧拉递推公式。
(1) 自启动 r阶多步算式无法自启动，需要用单步法求出前r步值
(2) 预估迭代 迭代方法要求收敛性良好，否则在大步长时会造成数 值发散。
(3) 变步长 初始阶段采用小步长，随后可逐步放大步长。
对不同精度要求的系统仿真，要考虑变阶次问题， 即为减小每一步计算的截断误差，以提高精度，应选用 较高的阶次，而当精度较低时，为减少工作量，则应选 取较低的阶次。仿真时应根据估计误差 与给定的误差精 度相比较改变步长或阶次来重新计算。
2 病态系统中绝对值最小的特征值对应于系统动态性能 解中瞬态分量衰减最慢的部分，它决定了整个系统的动 态过渡过程时间的长短。一般与系统中具有最小时间常 数Tmax的环节有关，要求计算步长h取得很大。
3 对于病态问题的仿真需要寻求更加合理的算法，以解 决病态系统带来的选取计算步长与计算精度，计算时间 之间的矛盾。
第四章 连续系统的离散化方法
4.1 常微分方程的数值解法
一. 数值求解的基本概念
设微分方程为
dx

dt

f
(t,
x)
x(t0 ) x0
则求解方程中函数x(t)问题的常微分方程初值问题
所谓数值求解就是要在时间区间[a, b]中取若干离散点
a t0 t1 L tN b hn tn1 tn

y1
y1( 0 ) 2, y2( 0 ) 0
1、建立M-文件vdp.m如下：
function dy=vdp(t,y) dy=zeros(2,1); dy(1)=y(2);
dy(2)=2*(1-y(1)^2)*y(2)-y(1);
2、取t0=0，tf=20，输入命令：
[T,Y]=ode45(‘vdp’,[0 10],[1;1]); plot(T,Y(:,1),'-‘, T,Y(:,2))
例 设系统方程为
g
x x2, x0 1
用Euler法求其数值解（取步长 h0.1, 0 t 1 )
递推公式为
xn1 xn hf tn , xn xn 1 0.1xn
则
t0 0, y0 1
t1 0.1, y1 y010.1y0 0.9 t2 0.2, y2 y110.1y1 0.90.91 0.819 t3 0.3, y3 y210.1y2 0.81910.10.819 0.7519
3、结果
练习 解微分方程组.

y1' y2 y3

y2 ' y1 y3 y3' 0.51y1 y2
y1(0) 0, y2 (0) 1, y3(0) 1
解 1、建立m-文件rigid.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);

t10 1.0, y10 y910.1y9 0.4628
已知方程的解析解为
y 1 1 t
精确解和解析解作比较：
t
0
精确解 1
数值解 1
0.1 0.2 0.3 0.4 0.5 …… 1.0
0.909 0.833 0.769 0.666 0.625 …… 0.5
0.9 0.819 0.752 0.659 0.647
大。
3.常微分方程Matlab求解
基于龙格－库塔法，MATLAB提供了求常微分方程数值解 的函数，一般调用格式为：
2/3阶龙格-库塔算法
[t, x]=ode23(‘xfun’, [t0, tf] ,x0)
4/5阶龙格-库塔算法 [t, x]=ode45(‘xfun’, [t0, tf] ,x0)
输入 输出
将 f tk b1h,xk hb2k1 在点 tk ,xk 展成Taylor级数
k2 fk b1hftk b2hk1 fxk
与台劳公式的二阶展开近似公式相比，可得以下关系：
ba12
a2
a2
1 12
a2b1 1 2
三个方程，四个未知数，解不唯一
0.463
误差在 10 2 数量级， 精度较差。
2. 龙格库塔法※
基本思想：取Taylor级数展开式前三项
近似求解，并利用线性组合代替导数的 求解。
既可避免计算高阶导数，又可提高
数值积分的精度，这就是Runge-Kutta
法的基本思想。
2.
龙格库塔法※
x(t0 h) x(t0 )
xk1 xk h fk
X&(t) AX (t) BU (t), X (t0) X 0
的系数矩阵A的特征值具有如下特征：
Re(i ) 0
m1iaxn Re(i ) ?
min
1i n
Re(i )
则称为“病态”方程。
2 控制系统仿真中的“病态”问题
1 病态系统中绝对值最大的特征值对应于系统动态性能 解中瞬态分量衰减最快的部分，它反映了系统的动态响 应和系统的反应灵敏度。一般与系统中具有最小时间常 数Tmin的环节有关，要求计算步长h取得很小。
j 1
1)当r=1时：
xk1 xk ha1k1,
k1 f (tk , xk ),
与Taloy展开式相比较，可得a1=1,则上式成为
xk1 xk hk1 xk hf (tk , xk ),
欧拉递推公式
2)当r=2时：
k2

k1 f tk ,xk f tk b1h, xk hb2k1
4.4 连续系统状态方程的离散化
上章所述的连续系统数学模型的离散化， 是通过数值积分法实现的，尽管面向结 构图的仿真方法是按环节给定参数，但 是在计算时还是按整个系统进行离散化， 这就不便于引进非线性环节以进行非线 性系统的仿真。在本节，将介绍连续系 统离散模型的建立和仿真。
常微分方 程函数名
起始 时间
终止 时间
初始状 态向量
例 1 经典的微分方程
解: 令 y1=x，y2=x’

d2x dt 2

2(1 x2 ) dx dt
x( 0 ) 1; x'( 0 ) 1

x

0
则微分方程变为一阶微分方程组：

y2
'

y1' y2 2(1 y12 )y2
(i 1, 2,L , n)
ki1 ki 2
ai1x1k ai2 x2k L ain xnk biu(tk )

ai1 ( x1k

h 2
ki1 )
L

ain (xnk

h 2
ki1 )
biu(tk

h) 2
ki3

ai1 ( x1k

h 2
ki2 )
L
仿真中遇到的大多数工程实际问题，四阶龙格库塔法以能 满足精度要求，其截断误差o(h5) 与h5同数量级。该法可以 自启动。
4)、状态空间四阶龙格-库塔递推式 若单输入单输出系统的状态空间表达式为：
X& AX BU Y CX DU
在仿真中，对于n阶系统，状态方程可以写成一阶微分方程
各个系数的几种取法——见书上。
3) r=4时，四阶龙格库塔公式-最常用：
h
xk 1

xk

( 6
K1

2K2

2K3

ቤተ መጻሕፍቲ ባይዱK4
)

K1 f tk ,xk

K2

K3

f f

tk

tk

h 2
,
xk
h 2 , xk