一种求解非良态高阶ODE初值问题的类解析法

陈 毅,伍儒康,王 伟,凌家胜

(南京工程学院机械工程学院, 江苏 南京 211167)

1 常微分方程的研究现状

微分方程的解是存在且唯一的,但由于理论方法的局限性,很多高阶微分方程无法求出解析解.从实际应用上讲,人们往往不需要解析解,而是关心某个定义范围内求出的对应近似值,即数值解[1-2].

常微分方程(简称ODE)的数值解法包括欧拉方法、龙格-库塔方法和差分法等,其中高阶(包括三阶)常微分方程一般可以化为一阶方程组的初值问题进行求解或者采用状态空间法运用经典的4级四阶龙格-库塔方法进行求解[3-5].但是差分法和龙格-库塔方法都存在求解的局限性,当微分方程三阶以上时,差分后的离散方程难以保证对角占优、绝对收敛;使用龙格-库塔方法求解非良态高阶常微分方程时,存在计算周期较长的问题[6-8].本文针对上述问题,结合微分方程解析解法与数值解法,研究微分方程的类解析法.

2 三阶ODE初值问题

定义三阶ODE的初值问题:

(1)

式中,y为输出,需要求解.

旋转机械发电装置导出的数学模型为:

(2)

式中:g(x)为输入;p、q、s′和d分别为三阶、二阶、一阶和y的系数;p≠0且g(x)≠0,即对应的三阶微分方程是非齐次的.

3 类解析法的介绍

3.1 输入函数的分解

周期为T的输入函数g(x)按三角形式的傅里叶级数展开为:

(3)

式中,a0、ai、bi为傅里叶系数.

解方程可得:

py‴+qy″+s′y′+dy=

(4)

展开有:

(5)

解出每一个方程对应的解析解,然后将上述2n+1项线性叠加.

3.2 特征方程的求解

式(4)对应的特征方程可以通过Matlab的solve命令求得,也可以通过迭代得到.特征根的情况分为四种:

情况1,px3+qx2+s′x+d=p(x+r1)(x+r2)×(x+r3),其中r1、r2、r3为特征方程的3个全不相等实根;

情况2,px3+qx2+s′x+d=p(x+r1)[(x+α)2+β2],1个实根、2个共轭复根;

情况3,3个相等实根;

情况4,3个实根,包含2个相等实根.

本文主要讨论情况1和情况2,情况3和情况4与情况1求解类似,且求解更简单,不再赘述.

3.3 拉氏变换与拉氏反变换

3.3.1 情况1

1) 输入为e1=a0/2,利用拉氏变换中微分定理[9]可得:

py‴+qy″+s′y′+dy=e1

(6)

变换为:

(7)

式中:e1为傅里叶多项式中的常数项a0/2;s为复数σ+jω;s′为微分方程一阶系数,即特征方程一次系数;r1、r2、r3为情况1的3个全不相等实根;k1、k2、k3、k4为分母连乘真分式化成部分真分式的系数.

将式(7)进行通分得到s的三次方程,s的三次、二次和一次系数应为0,常数项应等于1,4个方程可确定4个系数:

(8)

下面系数求解同理,不再赘述.通过拉氏反变换得到常数项输入对应的输出为:

(9)

(10)

系数可以通过Matlab矩阵求解得到:

(11)

通过拉氏反变换可得:

k3e-r2t+k4e-r3t)

(12)

(13)

式中:ei+n+1为第i个余弦输入的振幅,即bi;k11、k12、k2、k3、k4为部分真分式的系数.

系数可以通过Matlab矩阵求解得到:

(14)

通过拉氏反变换可得:

k3e-r2t+k4e-r3t)

(15)

式(15)与式(12)的区别在于求解的系数不同.

3.3.2 情况2

(16)

系数的求解为:

(17)

对应的拉氏反变换为:

(18)

(19)

系数通过Matlab矩阵求解得到：

(20)

同理,通过拉氏反变换可得:

(21)

3) 输入为ei+n+1coswit(i=1,2,…,n),其中

(22)

系数通过Matlab矩阵求解得到：

(23)

拉氏反变换为:

(24)

3.4 最终输出

将2n+1项线性叠加,得到最终的输出函数方程为:

(25)

当n较小时,这是一个有意义的解析解;当n较大时,可以根据步长得出数值解.

4 类解析法的运用

求解

式中:右边为系统输入;左边时变函数y(t)为系统输出.

已知系统输入为周期函数且周期为2.5 s,如图1所示,图1中仅展示了系统输入的5个周期.

图1 系统输入随时间变化图

很明显,上述三阶ODE各系数数量级相差较大,采用本文类解析法进行求解，步骤为:

1) 将系统输入进行傅里叶级数展开,展开阶数n=4,从图1中不难看出,9项傅里叶展开式可以很好地代替输入函数;

2) 运用Matlab进行求解特征根,求解结果为1个实根和2个共轭虚根;

3) 利用循环指令和矩阵运算求出常数输入的4个系数、正弦输入的20个系数和余弦输入的20个系数,共计44个常数,再回代到每个输出分量的表达式中,得出每项输入对应的输出;

4) 将输出分量线性叠加,得出系统输出,如图2所示,可知本方法收敛;

图2 系统输出随时间变化图

5) 利用得到的输出数据进行一阶差分、二阶差分和三阶差分,代入例题中与图1的系统输入进行比较,明显可见误差较小,验证了本方法求解结果的正确性.

5 结语

本文对三阶常系数ODE的数值求解方法进行探讨,给出了一个行之有效的求解方法——类解析法,即将三阶常系数ODE的输入按傅里叶级数展开成若干项,将一个输入的微分方程分解为若干输入分量的微分方程,分别求解并线性叠加,减少了单个方程的求解难度.通过计算可以证明,该方法适应性强、构造简单、可复制性强、求解准确且计算周期短,可以应用于机械、液压和电气等三阶系统中,是一个有效求解类解析解的数值解法.