一、基本定义与性质

1. Riemann-Liouville分数阶积分

定义式：
$$
I^\alpha f(t) = \frac{1}{\Gamma(\alpha)} \int_a^t (t-\alpha)^{\alpha-1} f(\tau) d\tau \quad (\alpha>0)
$$

性质：

  • 半群性：$\displaystyle I^\alpha I^\beta = I^{\alpha+\beta}$
  • 当$\alpha\to n^+$时收敛于n重积分
  • 对常数$C$的积分：$\displaystyle I^\alpha C = \frac{C(t-a)^\alpha}{\Gamma(\alpha+1)}$

2. Caputo分数阶导数

定义式：
$$
^C D^\alpha f(t) = \frac{1}{\Gamma(m-\alpha)} \int_a^t \frac{f^{(m)}(\tau)}{(t-\tau)^{\alpha+1-m}} d\tau
$$
其中$m-1 < \alpha \leq m$, $m\in\mathbb{N}$

优势：

  • 常数导数为零：$^C D^\alpha C = 0$
  • 允许传统初始条件：$\displaystyle f(a), f’(a),…,f^{(m-1)}(a)$

3. Grünwald-Letnikov定义

离散形式：
$$
D^\alpha f(t) = \lim_{h\to 0} h^{-\alpha} \sum_{k=0}^{[(t-a)/h]} (-1)^k \binom{\alpha}{k} f(t-kh)
$$

权重系数递推：
$$
w_0^{(\alpha)} = 1, \quad w_k^{(\alpha)} = \left(1-\frac{\alpha+1}{k}\right)w_{k-1}^{(\alpha)}
$$

本文我们主要关注Riemann-Liouville 分数阶积分与分数阶导数的性质

二、Riemann-Liouville 分数阶积分与分数阶导数的关系推导

1. Riemann-Liouville 分数阶积分的定义

设函数 $f(t)$ 在区间 $[a, t]$ 上定义，阶数 $\alpha > 0 $，则 Riemann-Liouville 左侧分数阶积分 定义为：
$$
I^\alpha f(t) = \frac{1}{\Gamma(\alpha)} \int_a^t (t - \tau)^{\alpha - 1} f(\tau) , d\tau,
$$

其中 $\Gamma(\alpha)$ 是 Gamma 函数。

2. Riemann-Liouville 分数阶导数的定义

设 $n = \lceil \alpha \rceil$ 是不小于 $ \alpha $ 的最小整数，则 Riemann-Liouville 分数阶导数 定义为：

$$
D^\alpha f(t) = \frac{d^n}{dt^n} \left( I^{n - \alpha} f(t) \right) = \frac{1}{\Gamma(n - \alpha)} \frac{d^n}{dt^n} \int_a^t (t - \tau)^{n - \alpha - 1} f(\tau) , d\tau.
$$

即：分数阶导数是对分数阶积分的整数阶导数。

3. 从积分推导到导数

设有：

$$
I^\alpha f(t) = \frac{1}{\Gamma(\alpha)} \int_a^t (t - \tau)^{\alpha - 1} f(\tau) , d\tau,
$$

我们希望找到某种操作 $D^\alpha$ 使得：

$$
D^\alpha \left( I^\alpha f(t) \right) = f(t).
$$

使用分数阶导数的定义，我们有：

$$
D^\alpha f(t) = \frac{d^n}{dt^n} \left( I_a^{n - \alpha} f(t) \right),
$$

那么对于 $f(t) = I^\alpha g(t)$，有：

$$
D^\alpha f(t) = D^\alpha \left( I^\alpha g(t) \right) = g(t),
$$

因为：

$$
D^\alpha \left( I^\alpha \right) = \text{Identity operator}.
$$

4. 总结关系

  • 分数阶积分是分数阶导数的逆操作：
    $$
    D^\alpha \left( I^\alpha f(t) \right) = f(t).
    $$

  • 分数阶导数可以表示为先做积分，再做整数阶导数：
    $$
    D^\alpha f(t) = \frac{d^n}{dt^n} \left( I^{n - \alpha} f(t) \right),
    $$
    其中 $n = \lceil \alpha \rceil$。

5. 特殊情况举例 $( 0 < \alpha < 1 )$

当 $0 < \alpha < 1$ 时，有：

  • 分数阶积分：
    $$
    I^\alpha f(t) = \frac{1}{\Gamma(\alpha)} \int_a^t (t - \tau)^{\alpha - 1} f(\tau) , d\tau.
    $$

  • 分数阶导数：
    $$
    D^\alpha f(t) = \frac{d}{dt} \left( \frac{1}{\Gamma(1 - \alpha)} \int_a^t (t - \tau)^{-\alpha} f(\tau) , d\tau \right).
    $$

6. 结论

Riemann-Liouville 分数阶导数与积分之间的本质关系是：

$$
D^\alpha = \frac{d^n}{dt^n} \circ I^{n - \alpha},
$$

也即：

  • 分数阶导数 = 整数阶导数 ∘ 分数阶积分；
  • 分数阶积分 = 分数阶导数的逆操作。

三、重要数学特性

1. 非局部性证明

对于任意$t>a$：
$$
D^\alpha f(t) = \frac{d^n}{dt^n}\left[\frac{1}{\Gamma(n-\alpha)} \int_a^t \frac{f(\tau)}{(t-\tau)^{\alpha+1-n}} d\tau\right]
$$
其中$n=\lceil\alpha\rceil$

2. Leibniz法则推广

$$
D^\alpha(fg) = \sum_{k=0}^\infty \binom{\alpha}{k} (D^{\alpha-k}f) D^k g
$$

3. Laplace变换关系

$$
\mathcal{L}{^C D^\alpha f(t)} = s^\alpha F(s) - \sum_{k=0}^{m-1} s^{\alpha-k-1} f^{(k)}(0^+)
$$

四、数值实现方法

1. 分数阶微积分的离散化方案

Grunwald-Letnikov改进格式：
$$
D^\alpha f(t_n) \approx \frac{1}{h^\alpha} \sum_{k=0}^{n} w_k^{(\alpha)} f(t_{n-k}) + \mathcal{O}(h^p)
$$
其中$p$为精度阶数，权重系数满足：
$$
w_0^{(\alpha)} = 1, \quad w_k^{(\alpha)} = \left(1-\frac{\alpha+1}{k}\right)w_{k-1}^{(\alpha)}
$$

2. 短记忆原理实现

MATLAB代码示例：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
function [dy] = fgl_deriv(a, y, h, L)
    % 参数说明：
    % a: 分数阶次 (0 < a < 1)
    % y: 输入信号序列
    % h: 采样步长
    % L: 记忆窗口长度
    
    n = length(y);
    dy = zeros(size(y));
    weights = zeros(n,1);
    weights(1) = 1;
    
    for k = 1:n-1
        weights(k+1) = (1 - (a+1)/k)*weights(k);
    end
    
    M = min(floor(L/h), n-1);
    for i = M+1:n
        dy(i) = sum(weights(1:M+1).*y(i:-1:i-M)')/h^a;
    end
end

3. 分数阶微分方程求解器

3.1 Predictor-Corrector算法完整推导

对于分数阶初值问题：
$$
\begin{cases}
^C D^\alpha_t y(t) = f(t,y(t)), \quad 0 < \alpha < 1 \
y(0) = y_0
\end{cases}
$$

离散化步骤：

  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
    $$

  2. 预测项计算：
    $$
    y^{P}{n+1} = 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
    $$

  3. 校正项计算：
    $$
    y_{n+1} = y_0 + \frac{h^\alpha}{\Gamma(\alpha+2)} \left[ f(t_{n+1},y^{P}{n+1}) + \sum{j=0}^n a_{j,n+1} f(t_j,y_j) \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 \leq j \leq n \
    1, & j=n+1
    \end{cases}
    $$

3.2 Python实现代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
import numpy as np
from scipy.special import gamma

def frac_pc_solver(f, alpha, y0, t_span, h):
    """
    Predictor-Corrector方法求解分数阶微分方程
    
    参数：
        f: 右端函数 f(t,y)
        alpha: 分数阶次 (0 < alpha < 1)
        y0: 初始条件
        t_span: 时间区间 [t0, tf]
        h: 步长
        
    返回：
        t: 时间节点数组
        y: 解向量
    """
    t0, tf = t_span
    t = np.arange(t0, tf + h, h)
    n = len(t) - 1
    y = np.zeros(n + 1)
    y[0] = y0
    
    # 预计算权重
    b_weights = [(k+1)**alpha - k**alpha for k in range(n+1)]
    a_weights = np.zeros((n+1, n+1))
    for j in range(n+1):
        if j == 0:
            a_weights[:,j] = np.arange(n+1)**(alpha+1) - \
                            (np.arange(n+1)-alpha)*(np.arange(n+1)+1)**alpha
        else:
            a_weights[j:,j] = (np.arange(n+1-j)+2)**(alpha+1) + \
                             (np.arange(n+1-j))**(alpha+1) - \
                             2*(np.arange(n+1-j)+1)**(alpha+1)
    
    for n_step in range(n):
        # Predictor
        y_p = y0 + (h**alpha/gamma(alpha+1)) * \
              sum(b_weights[n_step - k] * f(t[k], y[k]) for k in range(n_step+1))
        
        # Corrector
        sum_c = sum(a_weights[n_step+1, k] * f(t[k], y[k]) for k in range(n_step+1))
        y[n_step+1] = y0 + (h**alpha/gamma(alpha+2)) * \
                      (f(t[n_step+1], y_p) + sum_c)
    
    return t, y

3.3 收敛性分析

理论误差界

对于分数阶微分方程：

$$
{}^C D^\alpha y(t) = f(t, y(t)), \quad y(0) = y_0
$$

Predictor-Corrector 方法的全局误差为：

$$
\max_{0 \leq t_n \leq T} |y(t_n) - y_n| \leq C \cdot h^{\min{2,1+\alpha}}
$$

其中：

  • $C$ 依赖于 $\alpha$、$f$ 的 Lipschitz 常数和终态时间 $T$。
  • 当 $\alpha \to 1^-$，收敛阶趋近于 2（与经典 ODE 一致）。

数值验证实验

测试方程：

$$
{}^C D^{0.6} y(t) = \frac{\Gamma(4)}{\Gamma(3.4)} t^{2.4} - y(t), \quad y(0) = 0
$$

解析解：

$$
y(t) = t^3 - \frac{\Gamma(4)}{\Gamma(3.4)} \int_0^t (t - \tau)^{2.4} \tau^3 , d\tau
$$
python误差计算代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
import numpy as np
from scipy.special import gamma
from scipy.integrate import quad

def exact_solution(t, alpha=0.6):
    # 解析解计算（需数值积分）
    integral = quad(lambda tau: (t - tau)**(alpha - 1) * tau**3, 0, t)[0]
    return t**3 - gamma(4)/gamma(3.4) * integral / gamma(alpha)

def compute_errors(alpha=0.6, T=1, steps=[0.1, 0.05, 0.025, 0.0125]):
    errors = {}
    for h in steps:
        t, y = frac_pc_solver(lambda t, y: gamma(4)/gamma(3.4) * t**2.4 - y,
                             alpha, 0, [0, T], h)
        y_exact = np.array([exact_solution(ti) for ti in t])
        errors[h] = np.max(np.abs(y - y_exact))
    return errors

# 输出误差表
errors = compute_errors()
print("| 步长 h | 最大误差 E_inf | 收敛阶 p |")
print("|--------|----------------|----------|")
h_list = sorted(errors.keys())
for i in range(len(h_list)):
    h = h_list[i]
    if i == 0:
        p = "--"
    else:
        p = np.log(errors[h_list[i-1]] / errors[h]) / np.log(h_list[i-1] / h)
        p = f"{p:.2f}"
    print(f"| {h:.4f} | {errors[h]:.4e} | {p} |")

高阶方法扩展（Adams 型）

对于更高精度需求，可采用 分数阶 Adams-Bashforth-Moulton (ABM) 方法：

Predictor:

$$
y_{n+1}^P = y_0 + \frac{h^\alpha}{\Gamma(\alpha)} \sum_{j=0}^{n} b_{j,n+1} f(t_j, y_j)
$$

Corrector:

$$
y_{n+1} = y_0 + \frac{h^\alpha}{\Gamma(\alpha + 2)} \left[ f(t_{n+1}, y_{n+1}^P) + \sum_{j=0}^{n} a_{j,n+1} f(t_j, y_j) \right]
$$

其中权重系数：

$$
b_{j,n+1} = \frac{(n+1 - j)^\alpha - (n - j)^\alpha}{\alpha}
$$

$$
a_{j,n+1} = \frac{(n - j + 2)^{\alpha + 1} - 2(n - j + 1)^{\alpha + 1} + (n - j)^{\alpha + 1}}{\alpha + 1}
$$
python实现代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def frac_abm_solver(f, alpha, y0, t_span, h, corrector_steps=2):
    t = np.arange(t_span[0], t_span[1] + h, h)
    n = len(t) - 1
    y = np.zeros(n + 1)
    y[0] = y0
    
    # 计算权重
    def b_weight(k, alpha):
        return ((k+1)**alpha - k**alpha) / alpha
    
    def a_weight(k, alpha):
        return ((k+2)**(alpha+1) - 2*(k+1)**(alpha+1) + k**(alpha+1)) / (alpha+1)
    
    for i in range(n):
        # Predictor
        y_p = y0 + (h**alpha / gamma(alpha)) * sum(
            b_weight(i - j, alpha) * f(t[j], y[j]) for j in range(i+1))
        
        # 多步校正
        for _ in range(corrector_steps):
            sum_c = sum(a_weight(i - j, alpha) * f(t[j], y[j]) for j in range(i+1))
            y[i+1] = y0 + (h**alpha / gamma(alpha+2)) * (
                f(t[i+1], y_p) + sum_c)
            y_p = y[i+1]  # 迭代校正
    
    return t, y

3.4 改进算法变体

自适应步长控制

通过局部误差估计：

$$
\epsilon_{n+1} \approx \frac{\Gamma(\alpha + 2)}{h^\alpha} \left| y_{n+1} - y_{n+1}^P \right|
$$

步长调整策略：

$$
h_{\text{new}} = 0.9 h_{\text{old}} \left( \frac{\tau}{\epsilon_{n+1}} \right)^{\frac{1}{1 + \alpha}}
$$

高阶修正格式

采用加权 Grünwald 近似：

$$
D^\alpha y(t_n) \approx \frac{1}{h^\alpha} \sum_{k=0}^n w_k^{(\alpha)} y_{n-k} + \frac{1}{h^\alpha} \sum_{k=1}^p c_k \Delta^k y_n
$$

其中 $c_k$ 为修正系数，$\Delta^k$ 为阶差分算子。

差分算子 $\Delta^k$
$$
\Delta^k y_n = \sum_{i = 0}^{k} (-1)^i \binom{k}{i} y_{n - i}
$$

修正系数 $c_k$
$$
c_k = \frac{(-1)^k}{k!} \prod_{j = 1}^{k} (\alpha - j + 1)
$$

MATLAB高阶实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
function [t, y] = high_order_frac_solver(f, alpha, y0, t_span, h, p)
    t = t_span(1):h:t_span(2);
    n = length(t);
    y = zeros(1, n);
    y(1) = y0;
    
    % 计算 Grünwald 权重
    w = zeros(1, n);
    w(1) = 1;
    for k = 1:n-1
        w(k+1) = (1 - (alpha + 1)/k) * w(k);
    end
    
    % 计算差分修正系数
    c = zeros(1, p);
    for k = 1:p
        c(k) = (-1)^k / factorial(k) * prod(alpha - (0:k-1) + 1);
    end
    
    for i = 2:n
        % 基础 Grünwald 近似
        D_alpha = sum(w(1:i) .* y(i:-1:1)) / h^alpha;
        
        % 高阶修正项
        correction = 0;
        for k = 1:p
            if i > k
                delta_k = sum((-1).^(0:k) .* nchoosek(k, 0:k) .* y(i:-1:i-k));
                correction = correction + c(k) * delta_k / h^alpha;
            end
        end
        
        % 更新解
        y(i) = y(i-1) + h * f(t(i-1), D_alpha + correction);
    end
end

3.5 刚性系统处理

隐式格式的牛顿迭代法

  1. 隐式方程
    $$
    y_{n + 1} = y_0 + \frac{h^{\alpha}}{\Gamma(\alpha + 2)} \left[ f(t_{n + 1}, y_{n + 1}) + \sum_{j = 0}^{n} a_{j, n + 1} f(t_j, y_j) \right]
    $$
  2. 牛顿迭代步骤
    $$
    y_{n + 1}^{(k + 1)} = y_{n + 1}^{(k)} - \left[ I - \frac{h^{\alpha}}{\Gamma(\alpha + 2)} J_f(t_{n + 1}, y_{n + 1}^{(k)}) \right]^{-1} \cdot \text{Residual}
    $$
    其中 $J_f$ 为 Jacobi 矩阵。

python刚性系统求解

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def implicit_frac_solver(f, jacobian, alpha, y0, t_span, h, max_iter=10, tol=1e-8):
    t = np.arange(t_span[0], t_span[1] + h, h)
    n = len(t)
    y = np.zeros((n, len(y0)))
    y[0] = y0
    
    for i in range(1, n):
        # 显式预测
        y_pred = y[0] + (h**alpha / gamma(alpha + 1)) * sum(
            ((i - j + 1)**alpha - (i - j)**alpha) * f(t[j], y[j]) 
            for j in range(i))
        
        # 隐式校正（牛顿迭代）
        y_current = y_pred.copy()
        for _ in range(max_iter):
            residual = y_current - y[0] - (h**alpha / gamma(alpha + 2)) * (
                f(t[i], y_current) + 
                sum(a_weight(i - j - 1, alpha) * f(t[j], y[j]) for j in range(i)))
            
            J = np.eye(len(y0)) - (h**alpha / gamma(alpha + 2)) * jacobian(t[i], y_current)
            delta = np.linalg.solve(J, -residual)
            y_current += delta
            
            if np.linalg.norm(delta) < tol:
                break
        
        y[i] = y_current
    
    return t, y

# 辅助函数：计算隐式权重
def a_weight(k, alpha):
    return ((k + 2)**(alpha + 1) - 2 * (k + 1)**(alpha + 1) + k**(alpha + 1)) / (alpha + 1)

理论分析

1. 稳定性
  • 隐式格式对刚性系统的稳定性优于显式格式，稳定区域覆盖整个负半平面。
2. 计算复杂度
方法 每步计算量 适用场景
自适应步长 $O(n)$ 非光滑解
高阶修正 $O(n + p^2)$ 高精度需求
隐式迭代 $O(n \cdot d^3)$ 刚性系统 ($d$ 为维数)
3. 收敛性

$$
| y(t_n) - y_n | \leq C \cdot \left( h^{\min{2, 1 + \alpha}} + \tau \right)
$$
其中 $\tau$ 为牛顿迭代容差。

五、典型应用案例

1. 异常扩散建模

时间分数阶扩散方程：
$$
\frac{\partial^\gamma u}{\partial t^\gamma} = K_\gamma \frac{\partial^2 u}{\partial x^2}, \quad 0<\gamma<1
$$

解析解：
$$
u(x,t) = \frac{1}{2\sqrt{K_\gamma t^\gamma}} M_{1-\gamma/2}\left(\frac{|x|}{\sqrt{K_\gamma t^\gamma}}\right)
$$
其中$M_\beta(z)$是Wright函数。

2. 粘弹性材料本构关系

分数阶Kelvin-Voigt模型：
$$
\sigma(t) = E_\infty \epsilon(t) + E_\alpha D^\alpha \epsilon(t)
$$

3. 锂电池状态估计

分数阶等效电路模型：
$$
D^\alpha U_{oc}(t) = \frac{I(t)}{C} + R D^\beta I(t)
$$

参数辨识算法：

1
2
3
4
5
6
7
8
9
def fode_identify(data, alpha_range, beta_range):
    from scipy.optimize import minimize
    def cost(params):
        alpha, beta, C, R = params
        # 实现分数阶数值求解
        return error
    res = minimize(cost, x0=[0.5,0.5,1,1], 
                  bounds=[alpha_range, beta_range, (0,None), (0,None)])
    return res.x