基于Caputo导数的分数阶非线性振动系统响应计算

， *，

(1.天津大学 机械工程学院 力学系，天津 300072； 2.天津大学 天津市非线性动力学与混沌控制重点实验室，天津 300072)

1 引 言

分数阶微积分是对整数阶微积分的一种推广，将微积分运算的阶次从传统的整数阶扩展到了分数阶和复数阶，至今已有三百多年的发展历史，已成为研究反常扩散、多孔介质力学、非牛顿流体力学、粘弹性力学以及软物质物理等学科领域的有力数学工具。

关于分数阶导数的定义较多，常用的有Riemann-Liouville导数和Caputo导数。Riemann-Liouville导数对应的初始条件没有物理意义，而Caputo导数所描述系统的初始条件有明确的物理意义，且形式上和整数阶微分方程的初始条件类似，更符合工程实际，因而在有明确物理和工程背景的研究工作中应用较多[1-3]。对于Caputo导数下的微分方程，目前导数阶数相同时相关的算法已经成熟[4-10]，阶数不同时，文献[11]考虑了分段线性Chen系统的混沌和多稳态行为，但未给出具体的算法。为计算阶数不同时系统的响应，文献[2，3,17]采用了变通的做法，即用分数阶的Grunwald-Letnikov导数的算法来计算Caputo导数，这种处理是否可行，尚未见报道。

为此，本文研究阶数不同时Caputo导数意义下微分方程的算法。首先，给出一些预备知识；然后，推导含分数阶导数项非线性振动系统状态方程的标准形式及用于计算Caputo导数的迭代格式；最后，通过算例说明本文理论分析的正确性和所得计算方法的可靠性，以及用分数阶的Grunwald-Letnikov导数算法来计算Caputo导数时可能存在的问题。

2 预备知识

为叙述方便起见，先介绍本文用到的三种分数阶导数的定义及其联系。

定义2.1定义在区间[a,b]上的函数x(t)的p阶Caputo导数为[12]

(1)

式中n-1

定义2.2定义在区间[a,b]上的函数x(t)的p阶Riemann-Liouville导数为[13]

(2)

式中n-1

定义2.3函数x(t)的p阶Grunwald-Letnikov导数定义为[14]

Dp[x(t)]=dp[x(t)]/dtp|t=kh=

(3)

设x(t)定义在[a,b]上，若x(t)具有n-1阶连续导函数，且x(n)(t)∈L1(a,b)，则Caputo导数和Riemann-Liouville导数满足关系[15]

(4)

式中n-1

设x(t)定义在[a,b]上，若x(t)具有n阶连续导函数，则Riemann-Liouville导数与Grunwald-Letnikov导数等价[16]。

(5)

式中n-1

3 基于Caputo导数迭代格式的构造

在非线性振动系统的研究中，通常取初始时刻t=0，对应的Caputo导数为

(6)

式中n-1

3.1 含分数阶导数项非线性振动系统状态方程的标准形式

利用Caputo导数算子的叠加关系，可将含m个分数阶导数项的单自由度非线性振动系统化为标准形式的状态方程：

(7)

式中q=[q1,…,qm +2]T

X(t)=[x1(t),…,xm +2(t)]T

f[X(t),t]={f1[x(t),t],…,fm+2[x(t),t]}T

阶数qj满足0

先考虑周期激励下的含单个分数阶导数项的系统：

(8)

式中 0

在x(t)连续且分数阶导数存在的条件下，运用Caputo分数阶算子的叠加关系[17]：

(9)

进行增维处理[13,18-20]，可将系统(8)化为等价系统：

(10)

式中q=[q1,q2,q3]T

X(t)=[x1(t),x2(t),x3(t)]T

f[X(t),t]= {f1[x(t),t],f2[x(t),t],

f3[x(t),t]}T

当0

(11)

当1

(12)

考虑周期激励下含多个分数阶导数项的系统：

(13)

式中 0

同理，可将式(13)化为等价形式：

(14)

式中q1=p1,qj=pj-pj -1(j=2,…,l)

ql +1=1-pl，ql +2=pl +1-1

qj=pj -1-pj -2(j=l+3,…,m+1)

qm +2=2-pm,x(t)=[x1(t),…,xm +2(t)]T

且

(15)

上述分析表明，单自由度含分数阶导数项(单个或多个)的非线性振动系统，均可转化为标准形式(7)。多自由度非线性振动系统也可做同样的转化。因此，只要给出分数阶标准微分方程组(7)的算法，就可以解决含分数阶导数的非线性振动系统的数值求解问题。

3.2 一般迭代格式的构造

由第2节给出的三种导数之间的关系可知

(16)

式中 0

fj[x(tk -1),tk -1]

(17)

式中j=1,…,m+2，从而有

(18)

选取适当的时间步长h，通过以上迭代格式，可以得到x(t)的数值解，即Caputo导数意义下微分方程(7)的解。为方便起见，称该算法为Caputo算法。不要求状态方程中各分数阶导数阶数相等，是该算法的主要优点。

若去掉方程(18)右端第二项，即得到求GL导数意义下标准形式方程(7)解的迭代格式，以下简称GL算法。不难发现，只有在零初值条件下，两种算法的迭代格式才完全一致，对应的解才完全相同。而在非零初始条件下，两种算法的迭代格式形式不同，二者能否得到相同的结果是一个值得讨论的问题，关系到用GL算法计算Caputo分数阶导数的可行性。

4 数值仿真

通过若干算例对本文提出的算法进行验证。首先，通过有解析解的算例1～3说明本文计算方法的正确性，所选算例均来自文献[21]。然后，通过算例4说明用基于Grunwald-Letnikov导数算法计算Caputo导数下的非线性振动系统存在的问题。

算例1分数阶导数阶数小于1的情况：

(19)

其精确解为x(t)=1+t。

由Caputo分数阶导数算子的叠加关系[17]，对方程(19)进行相关处理，可将其转化为等价形式:

(20)

取时间步长h=0.001，用3.2节中的Caputo算法，可得方程(19)的数值解。数值解与解析解对比如图1所示，二者完全重合。

图1 算例1中数值解与解析解对比

Fig.1 Comparison of numerical and analytical solutions

of example 1

算例2分数阶导数阶数大于1的情形：

(21)

同理，可将方程(21)转化为

(22)

从图2可以看出,数值解(时间步长取h=0.001)与解析解同样吻合很好。

算例3含两个分数阶项的情形：

(23)

同上，可将方程(23)转化为

(24)

从图3可以看出，数值解(时间步长取h=0.0001)和解析解吻合很好。

上述算例表明，本文计算方法正确且具有很高的精度。

算例4分数阶Duffing振子系统[2]：

(25)

式中K1>0，0≤p≤2，m，k和c分别为系统的质量、线性刚度系数和线性阻尼系数，K1为分数阶导数项的系数，F1和w分别为外激励的幅值和频率，各参数的取值分别为m=5，k=45，c=0.2，α1=0.5，K1=0.1，p=0.1，F1=100和w=1.0612。

方程(25)可转化为状态方程：

对系统初值分三种情况讨论，分别将Caputo算法及GL算法所得结果进行比较(时间步长均取h=0.002)。

可以看出，初始条件均取(0, 0)时，两种算法等价，计算结果完全重合。

可以看出，两种算法所得结果经过不同过渡过程后，达到相同的稳态。

图2 算例2中数值解与解析解对比

Fig.2 Comparison of numerical and analytical solutions of example 2

图3 算例3中数值解与解析解对比

Fig.3 Comparison of numerical and analytical solutions

of example 3

可以看出，两种算法所得结果经过不同过渡过程后，分别达到两个不同的稳态解。此时GL算法不能得到正确的结果。

从本例可见，用GL算法来计算Caputo导数下的微分方程，除零初始条件下没有问题外，因初始条件不同，计算结果存在如下可能。(1) 仅过渡过程有差异，稳态结果正确。(2) 过渡过程不同，稳态结果也不同。这种现象可能与初始条件是否靠近吸引域的边缘有关，后续工作将对此展开详细讨论。

图4 零初始条件下两种方法的比较

Fig.4 Comparison of the two methods under the zero initial condition

图5 初始条件远离吸引域边界

Fig.5 Initial condition is far from the boundary of the basin of attraction

图6 初始条件靠近吸引域边界

Fig.6 Initial condition is close to the boundary of the basin of attraction

5 结 论

本文基于Caputo导数与Riemann-Liouville导数和Grunwald-Letnikov导数间的关系，建立了基于Caputo导数的数值离散方法，并通过数值算例讨论了所得算法的正确性及基于Grunwald-Letnikov导数和Caputo导数离散方法的差异，比较两种算法所得结果可得，用GL算法替代Caputo算法，不仅可能影响正确结果的计算，也可能影响吸引域的计算。在后续研究中，还要进一步考虑构造该算法的矩阵格式，从而利用Matlab矩阵运算高效的特点，提高其计算效率。