
引言
这篇文章整理自我在 2021 年完成的学位论文《分数阶系统的短记忆方法研究》,主要记录我当时围绕这个题目展开研究的过程,以及最后形成的一些理解。
我最早真正对分数阶系统产生兴趣,不是因为它"定义新奇",而是因为它有一种很特别的建模方式:系统当前的状态不只由此刻决定,还会持续受到整段历史过程的影响。也正因为这种记忆特性,分数阶模型在黏弹性材料、异常扩散、控制系统这些问题里,往往比普通整数阶模型更贴近实际。
这篇论文想解决的核心问题,其实很明确:分数阶系统虽然更能描述"带记忆"的动力学过程,但它也因此更难计算。像分数阶积分项 $\int_0^t (t-\tau)^{-\alpha}f(\tau)d\tau$ 这样的历史项,在数值计算里需要把从 0 到 $t$ 的全部历史离散后逐项累加。时间一长,历史项会越来越多,每往前走一步,计算量都会比前一步更大,最后整体算法的复杂度迅速上升。
围绕这个问题,我在论文里主要采用了两条方法线。第一条是预估校正法(Predictor-Corrector Method),用它来构造常分数阶与变分数阶系统的数值求解格式;第二条是短记忆原理(Short Memory Principle),用它把原本需要保留的完整历史积分截断为最近长度 $T$ 的一段记忆,从而把计算复杂度从 $O(N^2)$ 压到 $O(N \cdot T/h)$。论文的重点,不只是把这两种方法放在一起使用,而是分析它们结合之后,能在多大程度上同时兼顾精度与效率。
具体来说,这篇论文做了两件事。第一件事,是把预估校正法系统地应用到常分数阶系统和变分数阶系统中,并通过数值实验比较不同算法在精度上的表现。第二件事,是把短记忆策略引入这些算法,考察在截断历史积分之后,误差会如何变化、计算时间能缩短多少,以及记忆长度 $T$ 应该怎样选取才更合适。
最后真正得到的结论,也正是我后来觉得最值得单独写出来的部分:短记忆方法并不只是一个"为了提速而做的工程技巧",它背后还有一个更关键的定量问题,即分数阶阶次 $\alpha$ 与最优记忆长度 $T$ 之间到底存在什么关系。论文通过常分数阶和变分数阶两组实验说明,短记忆策略确实能够有效降低计算成本,而 $\alpha$ 的变化会直接影响 $T$ 的选取,这也是全文最核心的认识之一。
所以这篇文章会尽量保留我当时的思考路径,但引言先把综述框架交代清楚。为了让全文更适合博客阅读,下面按"概念基础 → 常分数阶算法 → 变分数阶推广 → 结论启发"的顺序展开。你可以把它理解为在回答三个问题:
- 分数阶系统为什么比整数阶系统更难算?
- 短记忆原理为什么能显著提速?
- 这种截断策略在常分数阶与变分数阶场景下分别能做到什么程度?
1. 分数阶系统的基础概念
本章先介绍后文会反复用到的几个核心概念。读到这里,最需要抓住的其实只有三件事:分数阶导数到底和普通导数有什么不同、预估校正法为什么会成为这类问题里最常见的数值工具、短记忆原理为什么有机会成立。
1.1 分数阶微积分的几种常见定义
在正式进入算法之前,先区分三种最常见的分数阶定义。它们并不是彼此竞争的三套理论,而更像是同一个问题在不同语境下的表达方式:有的更适合离散计算,有的更适合理论推导,有的更适合处理带初值的微分方程。后面论文真正用得最多的,其实是 Caputo 导数的建模形式,以及 GL 思路背后的离散化直觉。
1. 格伦瓦尔德-莱特尼科夫(GL)导数
格伦瓦尔德-莱特尼科夫导数是最容易和"离散历史求和"联系起来的一种定义。在数值计算里它尤其重要,因为它几乎直接告诉我们:分数阶导数不是只看当前点附近的局部变化,而是在把过去一串历史点按权重重新累加。$\alpha$ 阶 GL 导数写作:
$$ D_{GL}^{\alpha}f(t) = \lim_{h \to 0} h^{-\alpha} \sum_{j=0}^{\lfloor t/h \rfloor} (-1)^j \binom{\alpha}{j} f(t-jh) $$
其中二项式系数 $\displaystyle \binom{\alpha}{j} = \frac{\alpha(\alpha-1)(\alpha-2)\cdots(\alpha-j+1)}{j!}$。从这个形式里,已经能直观看到分数阶系统的"记忆"特征:当前时刻的导数值会吸收一整串过去状态,而不是像整数阶导数那样只由局部邻域决定。
2. Riemann-Liouville 积分
Riemann-Liouville 分数阶积分定义如下:
$$ I_{RL}^{\alpha}f(t) = \frac{1}{\Gamma(\alpha)} \int_0^t (t-\tau)^{\alpha-1} f(\tau) d\tau $$
其中 $\Gamma$ 为伽马函数。RL 积分的特点是可以对任意阶次 $\alpha>0$ 定义,且满足类微分的性质,常用于理论分析和解析求解。
3. Caputo 导数
Caputo 导数是另一种常用的分数阶定义:
$$ {}^C D_t^{\alpha}f(t) = \frac{1}{\Gamma(1-\alpha)} \int_0^t (t-\tau)^{-\alpha} f’(\tau) d\tau $$
在处理初值问题时,Caputo 导数通常更顺手,因为它的初值条件和整数阶微分方程保持一致,不需要额外引入 RL 型初始条件。这也是为什么后面在建立数值格式时,我会更自然地以 Caputo 形式作为出发点。
4. 三种定义的适用差异
| 定义 | 适用场景 | 初值处理 |
|---|---|---|
| GL 导数 | 数值离散计算 | 需要额外处理 |
| RL 积分/导数 | 理论分析、解析求解 | 需 RL 初始条件 |
| Caputo 导数 | 工程应用、初值问题 | 与整数阶相同 |
1.2 预估校正法
在真正讨论"历史能不能截断"之前,先要回答另一个更基础的问题:原方程本身该怎么稳定地数值求解。预估校正法(Predictor-Corrector Method)的思路并不复杂:先用一个较便宜的公式做预估步(P),再用更精确的公式做校正步(C)。它兼顾了实现难度、计算效率和精度,因此在分数阶数值计算里非常常见。
对于 Caputo 分数阶微分方程 $D^{\alpha}y(t) = f(t,y(t)),\ y(0)=y_0$,其积分形式为:
$$ y(t) = y_0 + \frac{1}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1} f(\tau,y(\tau)) d\tau $$
把积分形式离散化以后,就会得到后文反复出现的预估系数和校正系数。也就是说,后面所有看起来稍微复杂的递推公式,本质上都是从这里一步步展开的。后文会先在常分数阶系统中给出标准形式,再推广到变分数阶场景。
1.3 短记忆原理
有了基本的数值求解框架之后,接下来才会进入这篇论文真正关心的核心问题:历史积分到底能不能截短。分数阶积分 $\int_0^t (t-\tau)^{-\alpha}f(\tau)d\tau$ 会把整段历史输入做加权累积,但当 $\alpha \in (0,1)$ 时,越久远的历史影响通常越弱,因此就有理由去问:能不能只保留"最近一段"历史,而把更早的部分截掉?
短记忆原理将积分历史截断至最近长度 $T$ 的区间:
$$ I_{short}^{\alpha}f(t) \approx \frac{1}{\Gamma(\alpha)} \int_{t-T}^{t} (t-\tau)^{\alpha-1} f(\tau) d\tau $$
真正困难的不是"要不要截断",而是:$T$ 该取多大,才能在提速的同时把误差控制住? 这个问题看起来只是参数选择,实际上正是整篇论文最核心的研究入口。
2. 常分数阶系统的数值计算
2.1 预估校正法的基本思路
在研究变分数阶系统之前,论文先从较为简单的常分数阶系统入手。考虑最一般的常分数阶微分方程初值问题:
$$ {}^C D_t^\alpha y(t)=f(t,y(t)), \quad 0<\alpha<1,\quad y(0)=y_0 $$
从相关文献可知,这个方程等价于如下 Volterra 积分方程:
$$ y(t)=y_0+\frac{1}{\Gamma(\alpha)}\int_0^t (t-\tau)^{\alpha-1}f(\tau,y(\tau)),d\tau $$
论文在这里的思路,并不是直接从分数阶情形硬推,而是先回顾整数阶一阶微分方程里熟悉的预估校正方法,再把这种求积思想推广到常分数阶方程中。
先考虑本科阶段最熟悉的一阶微分方程初值问题。对于这类问题,如果假设方程在区间内存在唯一解,就可以用 Adams 预估校正方法求解。设区间长度为 $T$,总步数为 $N$,步长 $h=T/N$,其核心思想是:先通过一个较简单的预测公式得到预估值,再把这个预估值代回更精细的校正公式,得到下一步更精确的近似值。
有了这个铺垫之后,再回到常分数阶方程。由于公式中的积分区间从 $0$ 一直延伸到当前时刻,这体现了分数阶系统的非局部性,但并不妨碍继续使用预估校正的求积思想。下面取等距网格
$$ t_k=kh,\quad k=0,1,\dots,n+1 $$
并记 $y_k\approx y(t_k)$。在 $t_{n+1}$ 处,上式变成
$$ y(t_{n+1})=y_0+\frac{1}{\Gamma(\alpha)}\int_0^{t_{n+1}} (t_{n+1}-\tau)^{\alpha-1}f(\tau,y(\tau)),d\tau $$
接下来将积分区间划分为若干小区间,并在每个小区间上对被积函数做近似处理:
$$ \int_0^{t_{n+1}} = \sum_{j=0}^{n}\int_{t_j}^{t_{j+1}} $$
论文在这里采用的近似方式,是对函数 $g$ 在节点处做线性插值,并据此得到校正公式和预测公式。整理后,常分数阶微分方程的预估校正格式可以写为:
$$ y_{n+1}^{P} =y_0+\frac{h^\alpha}{\Gamma(\alpha+1)} \sum_{j=0}^{n} b_{j,n+1} f(t_j,y_j) $$
其中预测系数为
$$ b_{j,n+1}=(n+1-j)^\alpha-(n-j)^\alpha,\qquad j=0,1,\dots,n $$
校正公式为
$$ y_{n+1} =y_0+\frac{h^\alpha}{\Gamma(\alpha+2)} \left[ \sum_{j=0}^{n} a_{j,n+1}f(t_j,y_j) + f(t_{n+1},y_{n+1}^{P}) \right] $$
其中校正系数为
$$ a_{j,n+1}= \begin{cases} n^{\alpha+1}-(n-\alpha)(n+1)^\alpha, & j=0,\\ (n-j+2)^{\alpha+1}+(n-j)^{\alpha+1}-2(n-j+1)^{\alpha+1}, & 1\le j\le n,\\ 1, & j=n+1. \end{cases} $$
到此,常分数阶微分方程的预估校正公式已经建立完成。后续数值实验,正是围绕这组公式展开。
2.2 短记忆原理在常分数阶系统中的引入
在原始的预估校正公式中,不管是预估公式还是校正公式,除了最后一项以外,都只与当前时刻附近的状态有关;而最后的求和项由于积分区间从 $0$ 开始,一直到当前步为止,因此始终与全部历史函数值有关。这正是短记忆原理要处理的问题。
截断思想
论文采用的是固定积分长度的思路。设保留的积分长度为 $T$,当求和长度超过这个区间时,就从起点一侧截断更早的历史,只保留最近长度为 $T$ 的部分。
记
$$ M=\left\lfloor \frac{T}{h}\right\rfloor $$
如果当前步数 $n\le M$,说明当前积分区间还没有超过需要保留的长度 $T$,这时不需要截断,求和项仍然从 $0$ 开始累计。
如果 $n>M$,就截去更早的历史,只保留
$$ j=n-M,; n-M+1,;\dots,;n $$
这时预估步改写为
$$ y_{n+1}^{P} =y_0+\frac{h^\alpha}{\Gamma(\alpha+1)} \sum_{j=n-M}^{n} b_{j,n+1} f(t_j,y_j) $$
校正步改写为
$$ y_{n+1} =y_0+\frac{h^\alpha}{\Gamma(\alpha+2)} \left[ \sum_{j=n-M}^{n} a_{j,n+1}f(t_j,y_j) + f(t_{n+1},y_{n+1}^{P}) \right] $$
也就是说,当积分区间长度超过 $T$ 之后,求和下限不再固定为 $0$,而是改成从 $n-T/h$ 附近开始。
关于这样做为什么可行,可以从两个角度理解。
更具体地说,后面这两种分析方法,讨论的重点都不是短记忆原理本身是否成立,而是在实际计算中应该怎样选择保留的积分长度 $T$。
1. 解析方法
在计算截断误差时,若给定允许精度 $E$,理论上可以由相应误差公式反推出积分长度 $T$。论文据此指出:当参数值位于 $[0,1]$ 时,可以通过固定积分长度的方法减少计算量,而且所选积分长度与原始积分总长度无关;而当参数值超过这个范围时,固定积分长度的方法不再容易发挥作用,$T$ 的选取也会变得更困难。
2. 权函数图像分析法
除了理论误差分析之外,论文还给出了一种更直观的判断方式:重新观察预测系数 $b$ 和校正系数 $a$ 的图像分布。图像显示,当步数靠近当前状态时,权函数值会明显增大,而较早历史节点对应的权值较小;同时,参数值越小,前期权值越小。这说明参数越小时,系统对更早历史的依赖性越低,权函数图像也因此可以为积分长度 $T$ 的选取提供参考。
2.3 算法思路
在完成上述分析之后,论文将引入短记忆原理后的预估校正算法整理为一组实现步骤:
- 输入初始条件、步长和完整积分区间等已知参数。
- 结合权函数图像,确定保留积分长度 $T$。
- 按预测公式计算下一步的预估值。
- 将预估值代入校正公式,得到下一步的精确值。
- 当步数小于 $T/h$ 时,求和从 $0$ 开始。
- 当步数大于 $T/h$ 时,求和从 $n-T/h$ 开始。
至此,引入短记忆原理后的预估校正算法思路就建立完成,论文随后进入数值算例验证。
2.4 常分数阶算例与误差表现
2.4.1 算例一
$$D_*^\alpha x(t)=\frac{\Gamma(9)}{\Gamma(9-\alpha)}t^{8-\alpha}+\frac{9}{4}\Gamma(\alpha+1)+t^8+\frac{9}{4}t^\alpha-x(t)$$
初始条件: $x(0)=0,\ x’(0)=0$
精确解: $x(t)=t^8+\frac{9}{4}t^\alpha$
本文考虑 α=0.7 和 α=1.7 两种情况,对应 α 的两个重要区间 (0,1) 和 (0,2)。
区间 $(0,1)$,α=0.7,不同步长比较:
$h=1/160$

| $t$ | 精确值 $x(t)$ | 数值解 | 绝对误差 | 相对误差% |
|---|---|---|---|---|
| 0.1 | 0.4489 | 0.4563 | 0.0074 | 1.64 |
| 0.3 | 0.9687 | 0.9824 | 0.0136 | 1.41 |
| 0.5 | 1.3889 | 1.4064 | 0.0175 | 1.26 |
| 0.7 | 1.8105 | 1.8316 | 0.0211 | 1.16 |
| 0.9 | 2.5205 | 2.5494 | 0.0289 | 1.15 |
$h=1/320$

| $t$ | 精确值 $x(t)$ | 数值解 | 绝对误差 | 相对误差% |
|---|---|---|---|---|
| 0.1 | 0.4489 | 0.45343 | 0.00453 | 1.01 |
| 0.3 | 0.9687 | 0.97707 | 0.00837 | 0.86 |
| 0.5 | 1.3889 | 1.3997 | 0.01080 | 0.78 |
| 0.7 | 1.8105 | 1.8234 | 0.01290 | 0.71 |
| 0.9 | 2.5205 | 2.5382 | 0.01770 | 0.70 |
$h=1/640$

| $t$ | 精确值 $x(t)$ | 数值解 | 绝对误差 | 相对误差% |
|---|---|---|---|---|
| 0.1 | 0.4489 | 0.45169 | 0.00279 | 0.62 |
| 0.3 | 0.9687 | 0.97384 | 0.00514 | 0.53 |
| 0.5 | 1.3889 | 1.3955 | 0.00660 | 0.48 |
| 0.7 | 1.8105 | 1.8184 | 0.00790 | 0.44 |
| 0.9 | 2.5205 | 2.5314 | 0.01090 | 0.43 |
从以上的表格数据可以看到,在这个算例下,预估校正法得到的数值解精度较高,而且随着步长取值越来越小,计算精度也随之提高。
区间 $(0,2)$,α=1.7,$h=1/160$:

| $t$ | 精确值 $x(t)$ | 数值解 | 绝对误差 | 相对误差% |
|---|---|---|---|---|
| 0.1 | 0.0448934 | 0.044895 | 0.0000016 | 0.004 |
| 0.3 | 0.2906610 | 0.290670 | 0.0000090 | 0.003 |
| 0.5 | 0.6964250 | 0.696460 | 0.0000350 | 0.005 |
| 0.7 | 1.2846611 | 1.284700 | 0.0000389 | 0.003 |
| 0.9 | 2.3114931 | 2.311700 | 0.0002069 | 0.009 |
从图表中可以看到,当 α=1.7 时,误差几乎可以忽略不计,从该算例可以看到,本文提出的预估校正公式在 α∈(0,2) 时的精确程度较高。
2.4.2 算例二
$$D_*^\alpha x(t)=x^3\sin t-tx^2+t^2x-t^3,\ \alpha\in(0,1),\ x\in(0,1)$$
初始条件: $x(0)=0,\ x’(0)=0$
由于阿贝尔分数阶微分方程难以直接写出解析解,论文这里通过与文献中的其他算法数值结果对比,来观察预估校正算法的精确情况。
我们先通过权函数分析法来初步确定最少要保留的 T 值。选取步长 $h=1/5000$,则 $n=1/h=5000$,观察预测系数和校正系数的权函数图像:

从图中可以看到,当步数大于 4000 后,预测系数和校正系数的权函数值显著增加,因此可以大致确定保留的积分区间的步数为 $5000-4000 = 1000$,也就是保留的积分区间 $T$ 至少为 0.2。
我们先考察未引入短记忆原理的原始预估校正算法的精度情况,考察定义域为 $(0,1)$,以 Propsed method 的数值解结果作为基准值进行对比。
预估校正法与其他方法在各点的误差表:
| $t$ | Propsed method | 预估校正法 | 绝对误差 | 相对误差% |
|---|---|---|---|---|
| 0.2 | -0.00074434 | -0.00074432 | 0.00000002 | -0.00269 |
| 0.4 | -0.010522 | -0.010521 | 0.00000100 | -0.00950 |
| 0.6 | -0.05108 | -0.05107 | 0.00001000 | -0.01958 |
| 0.8 | -0.16592 | -0.16583 | 0.00009000 | -0.05424 |
| 1.0 | -0.46934 | -0.46864 | 0.00070000 | -0.14915 |
从图表中可以看到,未引入短记忆原理的预估校正法相较于其他方法依然能够保持优秀的精度。下面我们在预估校正法的基础上引入短记忆原理,设保留的积分区间长度为 $T$,可以得到不同 $T$ 值情况下的数值解。
短记忆方法结果:$T=0.2$ 至 $T=1.0$,给出 $x(1)$、绝对误差、相对误差%、运行时间/s
| $T$ | $x(1)$ | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|
| 0.2 | -0.29587 | 0.17347 | 36.96 | 11.46s |
| 0.3 | -0.37486 | 0.09448 | 20.13 | 24.04s |
| 0.4 | -0.42228 | 0.04706 | 10.03 | 24.78s |
| 0.5 | -0.44837 | 0.02097 | 4.47 | 27.73s |
| 0.6 | -0.46139 | 0.00795 | 1.69 | 37.51s |
| 0.7 | -0.467 | 0.00234 | 0.50 | 39.22s |
| 0.8 | -0.46891 | 0.00043 | 0.09 | 41.00s |
| 0.9 | -0.46927 | 0.00007 | 0.01 | 43.39s |
| 1.0(不截断) | -0.46864 | 0.00070 | 0.15 | 45.21s |

从图表可以看到,$T<0.3$ 时解的情况不佳,这与刚才通过预测校正权函数的分析结果相一致;当 $T=0.8$ 时误差已经足够小,因此我们重点研究 $T$ 取值 0.4~0.8 这个区间段的解的情况。
此外在实际测试中发现,对于预测方程来说,如果将保留的区间仅限为一个步长,最终得到的结果与原来的结果同样精确,因此我们可以通过这种方法大大降低计算量。定义此步骤为改进后的预估校正法。
改进方法结果:$T=0.4$ 至 $T=1.0$
| $T$ | $x(1)=$ | 绝对误差 | 相对误差 | 运行时间/s |
|---|---|---|---|---|
| 0.4 | -0.42169 | 0.04765 | 10.15 | 11.76s |
| 0.45 | -0.43678 | 0.03256 | 6.94 | 13.01s |
| 0.5 | -0.44771 | 0.02163 | 4.61 | 13.07s |
| 0.55 | -0.45543 | 0.01391 | 2.96 | 13.36s |
| 0.6 | -0.4607 | 0.00864 | 1.84 | 14.01s |
| 0.65 | -0.46416 | 0.00518 | 1.10 | 15.03s |
| 0.7 | -0.46632 | 0.00302 | 0.64 | 18.18s |
| 0.75 | -0.46757 | 0.00177 | 0.38 | 19.15s |
| 0.8 | -0.46822 | 0.00112 | 0.24 | 19.16s |
| 1.0(不截断) | -0.46864 | 0.00070 | 0.15 | 21.05s |

最优 $T$ 值选取:给定允许误差下的最优 $T$ 选择
| 允许相对误差/% | 最佳 $T$ | 相对误差/% | 运行时间/s | 保留所有积分区间下的运行时间/s |
|---|---|---|---|---|
| 10 | 0.4 | 10.15 | 11.76s | 45.21s |
| 5 | 0.5 | 4.61 | 13.07s | 45.21s |
| 2 | 0.6 | 1.84 | 14.01s | 45.21s |
| 1 | 0.65 | 1.10 | 15.03s | 45.21s |
从图表可以看到,通过这个算例可以看到,在运用了短记忆原理并不断改进后,在保证可允许范围内的精度下,大大减少了计算量,计算时间缩短了一倍有余。
2.4.3 算例三
$$D_*^\alpha x(t)=\dfrac{2}{\Gamma(3-\alpha)}t^{2-\alpha}-x(t)+t^2-t$$
初始条件:$x(0)=0,\ x’(0)=-1$
精确解:$x(t)=t^2-t$
由权函数分析法,先观察 α=1.5 时预测系数和校正系数的图像。可以发现,两个系数随自变量增大而逐渐减小,因此无法通过系数值图像来确定保留的积分区间长度 T。
于是我们采用实际截断测试的方式逐次确定 T。首先验证未引入短记忆原理时预估校正方法的解函数精度。
$t=10$ 至 $t=50$,精确值与数值解比较:
| $t$ | 精确解 | 数值解 | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|---|
| 10 | 90 | 90.1896 | 0.1896 | 0.2107 | 21.3s |
| 20 | 380 | 380.1303 | 0.1303 | 0.0343 | - |
| 30 | 870 | 870.1082 | 0.1082 | 0.0124 | - |
| 40 | 1560 | 1560.0952 | 0.0952 | 0.0061 | - |
| 50 | 2450 | 2450.0865 | 0.0865 | 0.0035 | - |
可以看到预估校正法得到的数值解与精确解近似度很高,因此可以认为预估校正法适用于本算例。
再考察保留的积分区间长度 $T$,设时刻 $t=50$,精确值 $x(50)=2450$,分数阶 $\alpha=1.5$,步长 $h=1/80$,令 $T$ 分别取 10,20,30,40,50(当 $T=50$ 即为不截取积分区间的情况)我们可以得到解函数图像。
短记忆方法:$T=10$ 至 $T=50$,$x(50)$ 及误差
| $T$ | $x(50)$ | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|
| 10 | 2411.74 | 38.26 | 1.562 | 3.96s |
| 20 | 2419.30 | 30.70 | 1.253 | 6.91s |
| 30 | 2434.18 | 15.82 | 0.646 | 9.04s |
| 40 | 2434.43 | 15.57 | 0.636 | 10.05s |
| 50 | 2450.09 | 0.09 | 0.004 | 10.83s |

从图表的数据可知当 $T=10$ 时,解函数图像仍有一部分震荡的情况,当 $T=20$ 及以上时,解函数图像基本稳定且误差情况足够优秀。因此再考察从 $T=10$ 开始的解函数图像如何,可以看到

从图表中可以看到当 $T=16$ 之后震荡的情况变优。下面再重点观察 $T=16$ 到 $T=20$ 的解函数情况。

$T=16$ 至 $T=20$ 详细结果
| $T$ | $x(50)$ | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|
| 16 | 2364.49 | 85.51 | 3.490 | 5.82s |
| 17 | 2418.09 | 31.91 | 1.302 | 5.85s |
| 18 | 2416.85 | 33.15 | 1.353 | 6.70s |
| 19 | 2415.32 | 34.68 | 1.416 | 9.64s |
| 20 | 2419.30 | 30.70 | 1.253 | 13.69s |
从表格中可以看到 T=17 的积分保留长度,计算得到的解精确度尚可,因此拟采用该长度的积分保留。
对比:非截断长记忆 $T=17$ vs $T=50$
| $T$ | $x(50)$ | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|
| 17(截断) | 2418.09 | 31.91 | 1.302 | 3.8s |
| 50(未截断) | 2450.0865 | 0.0865 | 0.004 | 21.3s |
通过这个算例同样可以发现,当 $T=17$ 时,同样可以运用短记忆原理大大减少计算时间。
从图表可以看到,预估校正法得到的数值解与精确解近似度很高。再看短记忆方法的各组结果,当 $T=10$ 时,解函数图像仍然存在较明显的震荡;当 $T=20$ 及以上时,图像基本稳定且误差已经较小。进一步细看 $T=16$ 到 $T=20$ 的情形后发现,当 $T=17$ 之后,函数的震荡情况明显改善,因此在兼顾函数收敛性的前提下,可以采用 $T=17$ 及以上的积分保留长度。
3. 变分数阶系统的算法推广
3.1 算法思路
变分数阶微分方程将常数阶次 $\alpha$ 推广为时间函数 $\alpha(t)$,即
$$ D_t^{\alpha(t)}x(t)=f(t,x(t)),\quad x(0)=x_0 $$
此时每个历史节点 $t_j$ 对应的阶次为 $\alpha(t_j)$,权重系数在每一步都需要重新计算。算法思路与常分数阶完全相同,只需在循环中将 $\alpha$ 替换为时变函数 $\alpha(t_j)$。引入短记忆原理后,保留最近 $T$ 时间区间内各节点对应的时变权重即可。
3.2 数值算例
3.2.1 预估校正法验证
首先不加短记忆,验证预估校正法对变分数阶系统是否依然有效。
算例I($\alpha(t)=t\cos t$):
$$ D_t^{\alpha(t)}x(t)+x(t)+\sqrt{t}x^2(t)=f(t),\quad x(0)=x’(0)=0 $$
其中 $f(t)=t^a\left(1+t^{\frac12+a}+\dfrac{\Gamma(1+a)t^{-\alpha(t)}}{\Gamma(1+a-\alpha(t))}\right),\ a=1.2$,精确解 $x(t)=t^a$。
取步长 $h=1/80,\ 1/160,\ 1/320$,数值解与精确解对比如下:
| $t$ | 解析解 | $h=1/80$ 数值解 | $h=1/80$ 相对误差% | $h=1/160$ 数值解 | $h=1/160$ 相对误差% | $h=1/320$ 数值解 | $h=1/320$ 相对误差% |
|---|---|---|---|---|---|---|---|
| 0.1 | 0.0631 | 0.07352 | 16.513 | 0.06644 | 5.293 | 0.06433 | 1.946 |
| 0.2 | 0.1450 | 0.14774 | 1.918 | 0.14585 | 0.614 | 0.14526 | 0.207 |
| 0.4 | 0.3330 | 0.33322 | 0.060 | 0.33300 | -0.006 | 0.33297 | 0.015 |
| 0.6 | 0.5417 | 0.54150 | -0.042 | 0.54157 | -0.030 | 0.54164 | 0.017 |
| 0.8 | 0.7651 | 0.76495 | -0.017 | 0.76498 | -0.013 | 0.76502 | 0.008 |
| 1.0 | 1.0000 | 1.00040 | 0.040 | 1.00010 | 0.010 | 1.00000 | 0.000 |
步长越小精度越高,可见 $\alpha\in(0,1)$ 时预估校正法对变分数阶系统同样有效。
算例II($\alpha(t)=t/2$):
$$ D_t^{\alpha(t)}x(t)=\frac{3t^{1-\alpha(t)}}{\Gamma(2-\alpha(t))}+\frac{2t^{2-\alpha(t)}}{\Gamma(3-\alpha(t))},\quad x(0)=x’(0)=0 $$
精确解 $x(t)=t^2+3t$。
| $t$ | 解析解 | $h=1/80$ 数值解 | $h=1/80$ 相对误差% | $h=1/160$ 数值解 | $h=1/160$ 相对误差% | $h=1/320$ 数值解 | $h=1/320$ 相对误差% |
|---|---|---|---|---|---|---|---|
| 0.2 | 0.6400 | 0.63993 | 0.011 | 0.63998 | 0.003 | 0.63999 | 0.002 |
| 0.4 | 1.3600 | 1.35970 | 0.022 | 1.35990 | 0.007 | 1.36000 | 0.000 |
| 0.6 | 2.1600 | 2.15940 | 0.028 | 2.15980 | 0.009 | 2.15990 | 0.005 |
| 0.8 | 3.0400 | 3.03910 | 0.030 | 3.03970 | 0.010 | 3.03990 | 0.003 |
| 1.0 | 4.0000 | 3.99920 | 0.020 | 3.99970 | 0.008 | 3.99990 | 0.003 |
结论与算例I一致,步长减少精度单调提升。
3.2.2 短记忆原理验证
在常分数阶分析中我们已经知道:当 $\alpha<1$ 时,预测和校正权函数在尾部几乎直线上升,$\alpha$ 越接近零则前部权重越小——这意味着短记忆截断的效果与 $\alpha$ 直接相关。下面通过变分数阶算例验证:$\alpha$ 越小,所需保留的积分长度越短。
算例I($\alpha(t)=t\cos t$),取 $h=1/5000$,逐步截断积分区间:
| $T$ | $x(1)$ | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|
| 0.2 | 0.7930 | 0.2071 | 20.71 | 9.09 |
| 0.3 | 0.8597 | 0.1403 | 14.03 | 12.37 |
| 0.4 | 0.9006 | 0.0994 | 9.94 | 15.29 |
| 0.5 | 0.9498 | 0.0502 | 5.02 | 17.87 |
| 0.6 | 0.9554 | 0.0446 | 4.46 | 20.30 |
| 0.7 | 0.9673 | 0.0327 | 3.27 | 22.45 |
| 0.8 | 0.9814 | 0.0186 | 1.86 | 23.70 |
| 0.9 | 0.9915 | 0.0085 | 0.85 | 23.89 |
| 1.0(不截断) | 1.0051 | 0.0051 | 0.51 | 24.03 |
对预测方程仅保留一个步长,精度几乎不受影响,据此改进后的结果:
| $T$ | $x(1)$ | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|
| 0.2 | 0.7984 | 0.2016 | 20.16 | 4.84 |
| 0.3 | 0.8653 | 0.1347 | 13.47 | 6.65 |
| 0.4 | 0.9064 | 0.0936 | 9.36 | 8.46 |
| 0.5 | 0.9340 | 0.0661 | 6.61 | 9.97 |
| 0.6 | 0.9554 | 0.0446 | 4.46 | 11.36 |
| 0.7 | 0.9732 | 0.0268 | 2.68 | 12.55 |
| 0.8 | 0.9878 | 0.0123 | 1.23 | 13.18 |
| 0.9 | 0.9989 | 0.0011 | 0.11 | 13.29 |
| 1.0(不截断) | 1.0053 | 0.0053 | 0.53 | 13.48 |
在允许误差范围内,改进方法在 $T=0.4\sim0.8$ 区间内均可保持较高精度。
| 允许相对误差/% | 最佳 $T$ | 实际相对误差/% | 运行时间/s | 不截断运行时间/s |
|---|---|---|---|---|
| 10 | 0.4 | 9.36 | 8.46 | 29.38 |
| 5 | 0.6 | 4.46 | 11.36 | 29.38 |
| 2 | 0.8 | 1.23 | 13.18 | 29.38 |
本算例中 $\alpha(t)=t\cos t\in[0,0.5611]$。若将阶次范围缩小一半至 $[0,0.5101]$(即 $\alpha(t)=\dfrac{t\cos t}{1.1}$),观察积分长度是否可以进一步缩短:
| $T$ | $x(1)$ | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|
| 0.2 | 0.8301 | 0.1700 | 17.00 | 5.59 |
| 0.3 | 0.8888 | 0.1112 | 11.12 | 13.76 |
| 0.4 | 0.9248 | 0.0752 | 7.52 | 16.15 |
| 0.5 | 0.9490 | 0.0511 | 5.11 | 18.17 |
| 0.6 | 0.9676 | 0.0324 | 3.24 | 22.73 |
| 0.7 | 0.9827 | 0.0173 | 1.73 | 23.31 |
| 0.8 | 0.9946 | 0.0054 | 0.54 | 26.87 |
| 0.9 | 1.0035 | 0.0035 | 0.35 | 27.85 |
| 1.0(不截断) | 1.0082 | 0.0082 | 0.82 | 28.11 |

相同 $T$ 值下相对误差明显减小,验证了 $\alpha$ 越小所需积分长度越短的猜想。
算例II($\alpha(t)=t/2$),取 $h=1/5000$,精确解 $x(t)=t^2+3t$:
| $T$ | $x(1)$ | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|
| 0.2 | 2.3370 | 1.6630 | 41.58 | 9.71 |
| 0.3 | 2.7806 | 1.2194 | 30.49 | 24.09 |
| 0.4 | 3.1170 | 0.8830 | 22.08 | 30.12 |
| 0.5 | 3.3802 | 0.6198 | 15.50 | 30.85 |
| 0.6 | 3.5877 | 0.4123 | 10.31 | 35.59 |
| 0.7 | 3.7495 | 0.2505 | 6.26 | 35.72 |
| 0.8 | 3.8722 | 0.1278 | 3.20 | 38.89 |
| 0.9 | 3.9576 | 0.0424 | 1.06 | 38.98 |
| 1.0(不截断) | 4.0000 | 0.0000 | 0.00 | 39.82 |
采用改进预测公式后精度不变,计算量大幅降低:
| $T$ | $x(1)$ | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|
| 0.2 | 2.3370 | 1.6630 | 41.58 | 4.36 |
| 0.3 | 2.7806 | 1.2194 | 30.49 | 10.00 |
| 0.4 | 3.1170 | 0.8830 | 22.08 | 10.09 |
| 0.5 | 3.3802 | 0.6198 | 15.50 | 10.97 |
| 0.6 | 3.5877 | 0.4123 | 10.31 | 19.20 |
| 0.7 | 3.7495 | 0.2505 | 6.26 | 21.53 |
| 0.8 | 3.8722 | 0.1278 | 3.20 | 23.38 |
| 0.9 | 3.9576 | 0.0424 | 1.06 | 23.48 |
| 1.0(不截断) | 4.0000 | 0.0000 | 0.00 | 25.08 |
| 允许相对误差/% | 最佳 $T$ | 实际相对误差/% | 运行时间/s | 不截断运行时间/s |
|---|---|---|---|---|
| 10 | 0.6 | 10.31 | 19.20 | 54.06 |
| 5 | 0.8 | 3.20 | 23.38 | 54.06 |
| 2 | 0.9 | 1.06 | 23.48 | 54.06 |
改进后计算量减少一倍以上,短记忆原理具有实用价值。
进一步将阶次范围从 $\alpha\in(0,0.5)$ 缩小至 $\alpha\in(0,\dfrac13)$(即 $\alpha(t)=t/3$),验证相同规律:
| $T$ | $x(1)$ | 绝对误差 | 相对误差% | 运行时间/s |
|---|---|---|---|---|
| 0.2 | 2.9046 | 1.0954 | 27.39 | 4.62 |
| 0.3 | 3.2430 | 0.7570 | 18.93 | 10.64 |
| 0.4 | 3.4803 | 0.5197 | 12.99 | 12.86 |
| 0.5 | 3.6539 | 0.3461 | 8.65 | 13.15 |
| 0.6 | 3.7822 | 0.2178 | 5.45 | 15.49 |
| 0.7 | 3.8758 | 0.1242 | 3.11 | 16.52 |
| 0.8 | 3.9416 | 0.0584 | 1.46 | 17.79 |
| 0.9 | 3.9829 | 0.0171 | 0.43 | 19.09 |
| 1.0(不截断) | 4.0000 | 0.0000 | 0.00 | 20.51 |

结果再次印证:$\alpha$ 越小,相同 $T$ 值下相对误差越小,短记忆效果越好。
4. 核心结论与启发
本文将短记忆原理与预估校正法结合,用于求解变分数阶系统。主要结论:
-
预估校正法在常、变分数阶系统中均有效:当 $\alpha\in(0,1)$ 时,无论是常数阶还是时变阶次,预估校正法均能保持较高精度。
-
短记忆截断能显著降低计算量:引入短记忆原理后,在保证精度的前提下可大幅减少计算时间——改进后的预估校正法计算量减少一倍以上。
-
$\alpha$ 与 $T$ 的定性关系:这是全文最重要的发现——$\alpha$ 越小,能够保留的积分长度越短。这意味着 $\alpha$ 较大时只需极短的记忆区间即可达到精度要求,而 $\alpha$ 较小时虽然所需记忆区间更长,但依然远小于全历史积分。该结论对不同 $\alpha$ 取值下的 $T$ 值选取具有实际参考意义。
-
变分数阶系统的短记忆推广可行:在每一步根据 $\alpha(t_j)$ 重新计算权重系数,即可将常分数阶的预估校正格式推广至变分数阶系统,短记忆策略同样有效。
这些结果表明,短记忆方法并非普适截断技巧,而是一种依赖阶次 $\alpha$、误差容忍度与问题结构的数值策略。其核心价值在于它提供了一个可分析、可实验验证、可推广至变分数阶系统的近似框架。
5. 参考文献
[1] Miller, K. S., & Ross, B. (1993). An Introduction to the Fractional Calculus and Fractional Differential Equations. Wiley.
[2] Diethelm, K. (2002). A Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations. Nonlinear Dynamics, 29(1-4), 3-22.
[3] Xu, Y., & He, Z. (2011). The short memory principle for solving Abel differential equation of fractional order. Journal of Computational and Applied Mathematics, 5-6.
[4] Parand, K., & Nikarya, M. (2015). New numerical method based on Generalized Bessel function to solve nonlinear Abel fractional differential equation of the first kind. Nonlinear Engineering, 5-7.
[5] Yousefi, F. S., Ordokhani, Y., & Yousefi, S. (2020). Numerical solution of variable order fractional differential equations by using shifted Legendre cardinal functions and Riez method. Engineering with Computers, 6.
[6] Patnaik, S., & Semperlotti, F. (2020). Application of variable and distributed order fractional operators to the dynamic analysis of nonlinear oscillators. Nonlinear Dynamics, 1-2.
[7] Sun, H. G., Chen, W., Wei, H., & Chen, Y. Q. (2011). A comparative study of constant-order and variable-order fractional models in characterizing memory property of systems. The European Physical Journal Special Topics, 2011.
[8] Ma, C. Y., Shiri, B., Wu, G. C., & Baleanu, D. (2018). New fractional signal smoothing equations with short memory and variable order. Signal Processing, 2018.
[9] Wu, F., Gao, R., Liu, J., & Li, C. (2020). New fractional variable-order creep model with short memory. Applied Mathematical Modelling, 2020.