第四章连续系统的离散化方法

相应的输出为 Yk 1 CX k 1 即可求得所需时刻
按上式，取 k 0,1, 2,, n
t1 , t1 ,, tn 的状态变量 X (tk )
和输出值 Y (tk ) 德国学者Felhberg对传统的龙格－库塔法提出了改进，在每一个计算步长内 对f( )函数进行了六次求值，以保证更高的精度和数值稳定性，假设当前的 步长为 ，定义下面 h 6个 Ki
Ki hf ( xk ij K j , tk i h), i 1, 2,,6;
j 1 i 1
j 1, 2,i 1
6
下一步的状态变量可由下式求出：
xk 1 xk i Ki
i 1
四阶/五阶龙格－库塔法系数表
i* ij i
0 16/135 25/216
xi ”的固定顺序输入i 1, 2, n 表示微分方程的序数。
，
4、输入参数“Tspan”规定了常微分方程的自变量取值范围， 它以矩阵[t0,tf]的形式输入，表示自变量 t t 0, tf
(t 0) x0 x(t 0) x x(t 0) 5、输入参数x0表示初始条件向量，
x x2 例 4－1 求解常微分方程
t2
(0) 1; t 0, 10 x(0) 1; x

0，初始条件为：
解：方法1 把二阶微分方程化成两个一阶微分方程组： 令 1 x2 则： x
2 x1 2 x t2
x1 x, x2 x

%定义输入、输出变量和函数文件名 %明确xdot的维数 % 第一个微分方程表示形式 % 第二个微分方程表示形式
以上几种递推公式均称为二阶龙格－库塔公式。是较典型的几 个常用算法。其中的方法3又称为预估－校正法，或梯形法。
x1 x0 hK1 其意义如下：用欧拉法以斜率先求取一点， 再由此点求得另一斜率 K2 f (t1 , x1 ) f (t0 h, x0 hK1 )
然后，从
1
x0 点开始，既不按该点斜率 K 变化，也不按预估点斜率
f (t , X ) AX BU 计算公式为：
h X k 1 X k ( K1 2 K 2 2 K 3 K 4 ) 6
其中， K1 f (tk , X k ) AX k BU (tk )
h h h h K 2 f (tk , X k K1 ) A( X k K1 ) BU (tk ) 2 2 2 2 h h h h K3 f (tk , X k K 2 ) A( X k K 2 ) BU (tk ) 2 2 2 2 K 4 f (tk h, X k hK3 ) A( X k hK3 ) BU (tk h)
K1 K 2 2
K 2 变化，而是取两者平均值 K h x1 x0 hK x0 ( K1 K 2 ) 2 h x1 x0 ( f 0 f1 ) 2
f
f1 f0
求得校正点，即：
。
0
t0
t1
t
四阶龙格－库塔法的计算公式为：
K1 f (tk , xk )
h xk 1 xk ( K1 2 K 2 2 K 3 K 4 ) 6
编制M文件，并且函数名和M文件名相同。 function xdot=vdpl(t,x) xdot=zeros(2,1); ％赋初值，并规定向量的维数。 xdot(1)=x(2); ％对第一个微分方程进行描述 xdot(2)=2*(1-x(1)^2)*x(2)-x(1); ％对第二个微分方程进行描 述 在命令窗口键入[t,x]=ode45('vdpl',[0,20],[1;1])，可得微分方程的 数值解，其前10组数据如下： x= t= 1.0000 1.0000 0 1.0489 0.9437 0.0502 1.0946 0.8762 0.1005 1.1368 0.7995 0.1507 1.1748 0.7158 0.2010 1.2441 0.5157 0.3136 1.2908 0.3156 0.4262 1.3159 0.1311 0.5389 1.3215 -0.0278 0.6515 1.3131 -0.1431 0.7484
将 K1 K 2 代入式
f f x1 x0 a1hf (t0 , x0 ) a2 h[ f (t0 , x0 ) b1h b2 hK1 ] t t t0 x x x0
a1 a2 1, a2b1 1 1 , a2b2 2 2
比较各项系数得
待定系数个数超过方程个数，必须先设定一个系数，然后即可求得其 参数。一般有以下几种取法： 1、 a1 0, a2 1, b1 b2
h h K 2 f (tk , xk K1 ) 2 2 h h K 3 f (tk , xk K 2 ) 2 2 K 4 f (tk h, xk hK 3 )
AX BU X
对于用状态方程表示的高阶线性系统 Y CX
X (0) X 0
其中状态变量为 X x1, x2 , xn 用四阶龙格－库塔法时，有
首先编制M文件，并且函数名和M文件名相同。 function xdot=wffc_1(t,x) %定义t,x,xdot和文件名 xdot=[0 1;-1 0]*x+[0;1]*(2+t^2/pi); %状态方程的表示形式 在命令窗口键入[t,x]=ode45(@wffc_1,[0,10],[-1;1])，可得微分方程的数值解， 其前10组数据如下： x= t= 0 0.0167 0.0335 0.0502 0.0670 0.1507 0.2344 0.3182 0.4019 0.5964 ------
0.8333333
0.819
0.7692307
0.5719
0.5
0.4627810
4.1.3 龙格－库塔法
K1 f (t0 , x0 )
h2 (2) (t0 ) x (t0 ) x(t0 h) x(t0 ) hx 2!
K2 f (t0 b1h, x0 b2 K1h) x1 x0 h(a1K1 a2 K2 )
x2 0.9 0.1 0.92 0.819 x3 0.819 0.1 0.8192 0.7519 x10 0.4627810
1 x 上式的精确解是 1 t
t 0 0.1
数值计算与精确解的比较见表
0.2 0.3 1.0
精确解x(t)
1
1
0.9090909
0.9
x x2 பைடு நூலகம்x (0) 1
解：因欧拉法的递推公式为
2 f (tk , xk ) xk
xk 1 xk hf (tk , xk )
所以
则递推公式为： 由t=0开始计算可得
2 xk 1 xk hxk
若取步长 h 0.1
x1 1 0.112 0.9
-1.0000 1.0000
-0.9828 1.0501 -0.9648 1.0999 -0.9460 1.1494 -0.9263 1.1986 -0.8158 1.4395 -0.6856 1.6709 -0.5363 1.8917 -0.3691 2.1007 0.0830 2.5349 ------
一般形式
3、a1 a2 , b1 b2 1 则 2
h x1 x0 ( K1 K 2 ) 2 K1 f (t0 , x0 )
一般形式
1
K 2 f (t0 h, x0 hK1 )
h xk 1 xk ( K1 K 2 ) 2 K1 f (tk , xk ) K 2 f (tk h, xk hK1 )
微分方程组中的方程个数必须等于初始条件数，这是求微分方程 特解所必须的条件。 6、输入参数options表示选项参数（包括tol，trace），可缺省， 3 即取默认值，tol是控制结果精度的选项对ode23( )函数取 10 ， 6 对ode45( )函数取10 。trace为输出形式控制变量，如果trace 不为0，则会将仿真中间结果逐步地由频幕显示出来，否则将不 显示中间结果 7、输出参数[t,x]为微分方程组解函数的列表（t和x都是列矩阵），它 包含向量t各节点 ti 和与ti 对应向量x的第j个分量值 x j g (ti ) （即第j个方程解），i表示节点序列数。 8、输出参数[t,x]缺省时，输出解函数的曲线，即函数 x g (t ) 及其各 ( j) 阶导数x j g (t ) 的曲线。 求解微分方程的指令还有ode113（多步解法器），ode15s（基于数字 微分公式的解法器），ode23s（单步解法器），ode23T（梯形规则的 一种自由插值实现），ode23TB（二阶隐式龙格－库塔公式）等。
首先编制M文件，并且函数名和M文件名相同。 function xdot=wffc_1(t,x) xdot=zeros(2,1); xdot (1)=x(1); xdot(2)=-x(1)+2+t^2/pi; 方法2 写出系统的状态方程
0 1 0 t2 x x (2 ) 1 0 1
1 3 2 , a2 , b1 b2 4 4 3
则
h xk 1 xk ( K1 3K 2 ) 4 K1 f (tk , xk ) 2 2 K 2 f (tk h, xk hK1 ) 3 3
h x1 x0 ( K1 3K 2 ) 4 K1 f (t0 , x0 ) 2 2 K 2 f (t0 h, x0 hK1 ) 3 3
' Fun '为定义微分方程组 2、输入参数
i fi (t , xi ) , x i 1 , 2,n
M－函数文件名，可以在文件名加写@，或用英文格式单引号界定文件名。
i 3、在编辑调试窗口中编写一阶常微分方程组 x fi (t , xi ) , i 1, 2,n
i fi (t , xi )一致，即等号 的M－函数文件时，每个微分方程的格式必须与 x 左边为带求函数的一阶导数，右边函数的变量 t , xi 严格以“先自变量t，后函数
x1 x0 hK 2 K1 f (t0 , x0 ) h h K 2 f (t0 , x0 K1 ) 2 2
1 则 2
一般形式
xk 1 xk hK 2 K1 f (tk , xk ) h h K 2 f (tk , xk K1 ) 2 2
2、 a1
。
例 4－2 求Var der Pol微分方程 (0) 下的数值 x(0) 1; x 1 解 t 0,
则 解：将二阶微分方程变换成一阶微分方程组x1 x, x2 x 1 x2 x
2 2*(1 x12 ) x2 x1 x

20
x ，在初始条件 x 2(1 x2 ) x 0
1/4
1/4
0
0
3/8
3/32
9/32
6656/12825
1408/2565
12/13
1932/2197
-7200/2197
7296/2197
28561/56430
2197/4104
1
439/216
-8
3680/513
-845/4104
-9/50
-1/5
1/2
-8/27
2
3544/2565
1859/4104
其一般公式为
xk 1 xk hf (tk , xk )
f1
f
f0
c
0
t0
t1
t
h2 h3 (2) h k ( k 1) Rn f (t0 , x0 ) f (t0 , x0 ) f (t0 , x0 ) 称为截断误差 2! 3! k! 例4－1 用欧拉法求下述微分方程的数值解。
-11/40
2/55
0
这一方法又称为四阶/五阶龙格－库塔法。
4.1.4 微分方程数值解的MATLAB实现
1、
t, x ode23(' Fun ',Tspan, x0, options) t, x ode45(' Fun ',Tspan, x0, options)
i fi (t , xi ), i 1, 2,n 该指令适用于一阶常微分方程组 x 如遇到高阶常微分方程，必须先将他们转换成一阶微分方程组，即状态方 程方可使用。
h2 f f x1 x0 hf (t0 , x0 ) ( f ) 2 t x t t0 , x x0
将 K 2在点 K1 f (t0 , x0 ) 处作台劳级数展开，并取线性部分可得
f f K 2 f (t0 , x0 ) b1 h b2 K1 h t t t0 x x x0