宇哥说的分数阶求导到底是个啥？

0.前言

相信各位研友在上张宇的课时，会听到宇神提到的分数阶求导。还记得当时在复习导数，一个『分数阶求导』的词不经意间跑过，但却把我的兴趣给抓牢了。

现在，我就来系统地叨叨分数阶求导的知识，也当是给各位消遣时的学习了。如有缺漏，还请指正，谢谢！

1. 什么是分数阶导数？

分数阶导数的通俗解释是，它是一种介于微分和积分之间的运算，它可以用来描述一些既不是局部的也不是全局的，而是介于两者之间的现象。

例如，如果一个物体的位移与时间的关系是 $x(t)=t^{\alpha }$，那么它的速度就是 $v(t)=\alpha t^{\alpha -1}$，它的加速度就是 $a(t)=\alpha (\alpha -1)t^{\alpha -2}$。如果 $\alpha$ 是一个整数，那么这些都是标准的导数。但是如果 $\alpha$ 是一个分数，那么这些可以被理解为分数阶导数，但并没有那么简单。

分数阶导数可以用来描述一些具有非线性或非平稳的特征的运动，例如分数阶谐振子，分数阶阻尼振荡器，分数阶布朗运动等。

分数阶导数有多种不同的定义，其中最常用的是 $Riemann-Liouville$ 分数阶导数和 $Caputo$ 分数阶导数¹：

$$D^{\alpha }f(x)=\frac{1}{\Gamma (n-\alpha )}\frac{d^n}{d{x}^n}{\int }_0^x\frac{f(t)}{(x-t{)}^{\alpha -n+1}}dt,n-1<\alpha <n$$

$${}_a{D}_x^{\alpha }f(x)=\frac{1}{\Gamma (n-\alpha )}{\int }_a^x\frac{f^{(n)}(t)}{(x-t{)}^{\alpha -n+1}}dt,n-1<\alpha <n$$

其中 $\alpha$ 是任意正实数，$\Gamma$ 是 $Gamma$ 函数，$n$ 是大于 $\alpha$ 的最小整数，$f^{(n)}$ 是 $f$ 的 $n$ 阶导数。这两种定义的区别在于 $Riemann-Liouville$ 分数阶导数是对原函数 $f$ 进行分数阶积分，然后再求整数阶导数，而 $Caputo$ 分数阶导数是是对 $f$ 的最高整数阶导数积分，即，先对原函数 $f$ 求整数阶导数，然后再进行分数阶积分。$Caputo$ 分数阶导数的优点是它可以保持原函数的初值条件，而 $Riemann-Liouville$ 分数阶导数则需要引入分数阶初值条件²。

2. 怎样理解分数阶导数？

相信看到这一大串公式以后，大家都会比较茫然，不知道该怎样理解。我们可以从它的最基本的物理意义和几何意义着手进行理解。

2.1 物理意义

所以，请让我们回顾一下整数阶导数的定义，它是通过极限的概念来定义的，例如一阶导数的定义是：

$$f^{\prime }(x)=\underset{h\to 0}{\lim }\frac{f(x+h)-f(x)}{h},$$

它表示函数 $f$ 在点 $x$ 处的变化率，也就是切线的斜率。如果 $f$ 表示位移，那么 $f^{\prime }$ 就表示速度。二阶导数的定义是：

$$f^{\prime \prime }(x)=\underset{h\to 0}{\lim }\frac{f^{\prime }(x+h)-f^{\prime }(x)}{h},$$

它表示函数 $f^{\prime }$ 在点 $x$ 处的变化率。如果 $f^{\prime }$ 表示速度，那么 $f^{\prime \prime }$ 就表示加速度。这些都是整数阶导数的情况，它们都是基于极限的概念，也就是说，它们只考虑了函数在某一点附近的局部性质，而忽略了函数在其他地方的全局性质³。这样的导数可以用来描述一些简单的运动，例如匀速或匀加速运动，但是对于一些更复杂的运动，例如非线性或非平稳的运动，这样的导数就不够用了。

分数阶导数的定义是通过积分的概念来定义的，它是通过对函数进行多次积分，然后再求导的方式来得到的，例如 $Riemann-Liouville$ 分数阶导数的定义是：

$$D^{\alpha }f(x)=\frac{1}{\Gamma (n-\alpha )}\frac{d^n}{d{x}^n}{\int }_0^x\frac{f(t)}{(x-t{)}^{\alpha -n+1}}dt,n-1<\alpha <n,$$

它表示函数 $f$ 在点 $x$ 处的 $\alpha$ 阶导数，其中 $\alpha$ 是任意正实数，$n$ 是大于 $\alpha$ 的最小整数，$\Gamma$ 是 $Gamma$ 函数。

这个定义的含义是，先对函数 $f$ 在区间 $[0,x]$ 上进行 $n-\alpha$ 次积分，得到一个新的函数 $F$，然后再对 $F$ 在点 $x$ 处求 $n$ 阶导数，得到 $D^{\alpha }f(x)$。这样的导数可以用来描述一些复杂的运动，例如具有记忆效应或历史依赖性的运动，因为它不仅考虑了函数在某一点的局部性质，还考虑了函数在其他地方的全局性质。这是因为积分的过程相当于对函数的历史进行了累积，而导数的过程相当于对函数的当前进行了微分。所以分数阶导数可以看作是一种结合了积分和导数的运算，它可以捕捉到函数的更多信息，而不是只关注函数的一阶或二阶变化。

2.2 几何意义

比如，对于整数阶导数的图形，它们都是与曲线的切线或曲率有关的。

例如一阶导数表示曲线 $y=f(x)$ 在点 $(x,f(x))$ 处的切线的斜率，也就是曲线在该点的方向。如果曲线是直线，那么它的一阶导数就是常数，表示曲线的斜率不变。如果曲线是曲线，那么它的一阶导数就是变化的，表示曲线的斜率随着位置的变化而变化。

二阶导数的表示曲线 $y=f(x)$ 在点 $(x,f(x))$ 处的曲率，也就是曲线在该点的弯曲程度。如果曲线是直线，那么它的二阶导数就是零，表示曲线的曲率不存在。如果曲线是圆，那么它的二阶导数就是常数，表示曲线的曲率不变。如果曲线是曲线，那么它的二阶导数就是变化的，表示曲线的曲率随着位置的变化而变化。这些都是整数阶导数的情况，它们都是基于切线或曲率的概念，也就是说，它们只反映了曲线在某一点的局部性质，而忽略了曲线的整体形状。这样的导数可以用来描述一些简单的曲线，例如直线或圆，但是对于一些更复杂的曲线，例如分形或混沌，这样的导数就不够用了。

分数阶导数的图形是与曲线的整体形状有关的，它们不是切线或曲率，而是一种更抽象的概念，例如 $Riemann-Liouville$ 分数阶导数的图形是：

它表示曲线 $y=f(x)$ 在点 $(x,f(x))$ 处的 $\alpha$ 阶导数，其中 $\alpha$ 是任意正实数。

这个图形的含义是，先对曲线 $y=f(x)$ 在区间 $[0,x]$ 上进行 $n-\alpha$ 次积分，得到一个新的曲线 $y=F(x)$，然后再对 $F$ 在点 $(x,F(x))$ 处求 $n$ 阶导数，得到 $D^{\alpha }f(x)$。这样的导数可以用来描述一些复杂的曲线，例如具有非整数维度或非马尔可夫性质的曲线，因为它不仅反映了曲线在某一点的局部性质，还反映了曲线的整体形状。这是因为积分的过程相当于对曲线的形状进行了平滑，而导数的过程相当于对曲线的形状进行了放大⁴。所以分数阶导数可以看作是一种结合了平滑和放大的运算，它可以反映出曲线的更多特征，而不是只关注曲线的一阶或二阶变化。

3. 计算方法

3.1 近似为有限差分

要计算分数阶导数，有多种数值方法，其中一种是基于 $Grunwald-Letnikov$ 公式的方法，它可以将分数阶导数近似为有限差分的形式⁵：

$$D^{\alpha }f(x)\approx \frac{1}{h^{\alpha }}\sum _{k=0}^N(-1{)}^k(\genfrac{}{}{0ex}{}{\alpha }{k})f(x-kh)$$

其中 $h$ 是步长，$N$ 是取决于 $x$ 和 $h$ 的整数，$\left(\genfrac{}{}{0ex}{}{\alpha }{k}\right)$ 是广义二项式系数，定义为：

$$\left(\genfrac{}{}{0ex}{}{\alpha }{k}\right)=\frac{\alpha (\alpha -1)\cdots (\alpha -k+1)}{k!},k=0,1,2,\dots$$

我用了Matlab和C++分别对这个公式进行了代码实现。以下是一个用 MATLAB 编写的函数，它可以计算任意函数 $f$ 在任意点 $x$ 处的分数阶导数：

function y = fgl_deriv(f, x, alpha)
% f: a function handle
% x: a scalar value
% alpha: a positive real number
% y: the fractional derivative of f at x of order alpha

h = sqrt(eps); % choose a small step size
N = floor(x / h); % choose the number of terms
y = 0; % initialize the output
for k = 0:N % loop over the summation
    y = y + (-1)^k * nchoosek(alpha, k) * f(x - k * h); % add each term
end
y = y / h^alpha; % divide by the power of h
end

或者我们用C++来编写：

#include <cmath>
#include <functional>

double fgl_deriv(std::function<double(double)> f, double x, double alpha)
{
    // f: a function object
    // x: a scalar value
    // alpha: a positive real number
    // return: the fractional derivative of f at x of order alpha

    double h = std::sqrt(std::numeric_limits<double>::epsilon()); // choose a small step size
    int N = std::floor(x / h); // choose the number of terms
    double y = 0; // initialize the output
    for (int k = 0; k <= N; k++) // loop over the summation
    {
        y += std::pow(-1, k) * std::tgamma(alpha + 1) / (std::tgamma(k + 1) * std::tgamma(alpha - k + 1)) * f(x - k * h); // add each term
    }
    y /= std::pow(h, alpha); // divide by the power of h
    return y;
}

3.2 基于具体问题的求解

举一个分数阶微分方程的例子：

$${}_{\mathrm{C}}{D}^{0.5}y(t)+y(t)=1,y(0)=0.$$

这是一个分数阶阻尼振荡器的模型，它描述了一个具有分数阶粘弹性的弹簧的运动。它的解析解是：

$$y(t)=1-E_{0.5}(-t^{0.5}),$$

其中 $E_{0.5}$ 是 $Mittag-Leffler$ 函数，它是分数阶微分方程的常见解。这个方程的数值解可以用 MATLAB 的 ode45 函数求得：（要用到前面定义的fderiv函数哦）：

% define the function handle for the fractional derivative
f = @(t,y) 1 - y - 0.5 * fderiv(@(s) y, t, 0.5, 0.01);

% solve the fractional differential equation using ode45
[t,y] = ode45(f, [0,10], 0);

% plot the solution
plot(t,y)
xlabel('t')
ylabel('y')
title('Solution of the fractional differential equation')

3.3 几个特殊的阶数

对于特殊阶数的分数阶求导， $\frac{1}{2}$ 阶导数是再合适不过的例子了，在这里给各位站友们介绍如何求解它们的简便方法。

首先，分数阶导数的定义可以通过多种方式给出，最常见的是$Riemann-Liouville$分数阶导数和$Caputo$分数阶导数（老生常谈的话题了）。对于函数$f(t)$，其$Riemann-Liouville$分数阶导数定义为：

$$D_{a+}^{\alpha }f(t)=\frac{1}{\Gamma (n-\alpha )}\frac{d^n}{d{t}^n}{\int }_a^t\frac{f(\tau )}{(t-\tau {)}^{\alpha -n+1}}d\tau ,n-1<\alpha <n$$

其中，$\Gamma$是伽马函数，$n$ 是大于$\alpha$的最小整数。

对于$1/2$阶导数，我们可以使用上述定义，其中$\alpha =1/2$，并且$n=1$。

于是：

$$D_{a+}^{1/2}f(t)=\frac{1}{\Gamma (1/2)}\frac{d}{dt}{\int }_a^t\frac{f(\tau )}{\sqrt{\pi (t-\tau )}}d\tau$$

求解分数阶导数的方法通常涉及到变换方法，如拉普拉斯变换和傅里叶变换，以及特殊函数，如$Mittag-Leffler$函数。

以下是一些特殊函数的$1/2$阶导数的表格：

函数	$1/2$阶导数
$e^t$	$\frac{e^t}{\sqrt{\pi t}}$
$t^n$	$\frac{t^{n-\frac{1}{2}}}{\Gamma (n+\frac{1}{2})}$
$\sin (t)$	$\frac{\cos (t)}{2\sqrt{\pi t}}-\frac{\sin (t)}{2\sqrt{\pi }}\Gamma (\frac{1}{2},t)$
$\cos (t)$	$\frac{\sin (t)}{2\sqrt{\pi t}}-\frac{\cos (t)}{2\sqrt{\pi }}\Gamma (\frac{1}{2},t)$

同样，对于$\frac{1}{3}$阶导数，直接把值代入就可以了。

$$D_{a+}^{1/3}f(t)=\frac{1}{\Gamma (2/3)}\frac{d}{dt}{\int }_a^t\frac{f(\tau )}{(t-\tau {)}^{2/3}}d\tau$$
给出一些特殊函数的$1/3$阶导数的表格：

函数	$1/3$阶导数
$e^t$	$\frac{3e^t}{2\Gamma (\frac{2}{3})t^{1/3}}$
$t^n$	$\frac{t^{n-\frac{1}{3}}}{\Gamma (n+\frac{2}{3})}$
$\sin (t)$	$\frac{3\cos (t)}{2\Gamma (\frac{2}{3})t^{1/3}}-\frac{\sin (t)}{2\Gamma (\frac{2}{3})}\Gamma (\frac{1}{3},t)$
$\cos (t)$	$\frac{3\sin (t)}{2\Gamma (\frac{2}{3})t^{1/3}}-\frac{\cos (t)}{2\Gamma (\frac{2}{3})}\Gamma (\frac{1}{3},t)$

4. 分数阶导的应用

4.1 初析

分数阶导数的求解应用涉及到分数阶微分方程，即包含分数阶导数的方程。分数阶微分方程可以用来建模许多自然科学和工程领域的问题，例如流体力学，电路理论，粘弹性，生物物理，信号处理，图像处理，金融数学等。

以下是一些具体的例子：

1. 在电路和网络中，分数阶导数可以用来模拟一些具有分数阶阻抗或电容的元件，例如分数阶电感器、分数阶电容器、分数阶电阻器等⁶。这些元件可以提高电路的灵敏度、带宽、稳定性等性能。例如，以下是一个含有分数阶电感器的电路的微分方程：

$$L^{\alpha }\frac{d^{\alpha }i}{d{t}^{\alpha }}+Ri=V$$

• 其中 $L^{\alpha }$ 是分数阶电感器的分数阶电感，$R$ 是电阻，$V$ 是电压，$i$ 是电流，$\alpha$ 是分数阶参数。这个方程可以用数值方法求解，例如 $Grunwald-Letnikov$ 公式⁷。

2. 在非局部效应中，分数阶导数可以用来描述一些具有长程相互作用或历史依赖的系统，例如粘弹性材料、变粘性流体、热传导、波动等⁸。这些系统的动力学行为不能用整数阶微分方程来刻画，而需要用分数阶微分方程来描述。例如，以下是一个描述粘弹性材料的分数阶 $Kelvin-Voigt$ 模型的微分方程：

$$E_0\frac{d^{\alpha }u}{d{t}^{\alpha }}+E_{1}u=F$$

• 其中 $E_0$ 和 $E_1$ 是材料的弹性模量，$u$ 是位移，$F$ 是外力，$\alpha$ 是分数阶参数。这个方程可以用数值方法求解，例如 $Adams-Bashforth-Moulton$ 方法⁹。

3. 在概率和统计中，分数阶导数可以用来建立一些具有非高斯分布或非马尔可夫性的随机过程，例如分数阶布朗运动、分数阶扩散过程、分数阶随机微分方程等¹⁰。这些过程可以用来模拟一些具有异常扩散或长记忆的现象，例如生物细胞的运动、金融市场的波动、地震的发生等¹¹。例如，以下是一个描述分数阶布朗运动的随机微分方程：

$$d{X}_t=\mu dt+\sigma d{W}_t^{\alpha }$$

• 其中 $X_t$ 是分数阶布朗运动的位置，$\mu$ 是漂移系数，$\sigma$ 是扩散系数，$W_t^{\alpha }$ 是分数阶维纳过程，$\alpha$ 是分数阶参数。这个方程可以用数值方法求解，例如 $Euler-Maruyama$ 方法¹²。

4. 在分形和随机介质中，分数阶导数可以用来刻画一些具有自相似性或多尺度性的结构，例如分形维数、分形维度谱、分形布朗运动、分形随机场等 。这些结构可以用来描述一些具有复杂几何形状或不均匀性质的物理系统，例如多孔介质、裂缝岩石、分形图像、分形信号等 。例如，以下是一个描述分形维数的公式：

$$D_f=\underset{ϵ\to 0}{\lim }\frac{\log N(ϵ)}{\log (1/ϵ)}$$

• 其中 $D_f$ 是分形维数，$N(ϵ)$ 是用边长为 $ϵ$ 的小方块覆盖分形所需的最小数量。这个公式可以用数值方法求解，例如盒计数法 。

4.2 方法介绍

所以，光看别人是怎么在实际工程问题应用的，还是感觉有点拿捏不住。那这里就具体地来介绍一种方法吧。

分数阶微分方程的求解方法其实有 $1mol$ 多,其中一种常用的方法是基于 $Gr\ddot{u}nwald-Letnikov$ 公式的数值方法。$Gr\ddot{u}nwald-Letnikov$ 公式是一种用差分近似分数阶导数的公式：

$${}_a{D}_t^{\alpha }f(t)\approx \frac{1}{h^{\alpha }}\sum _{j=0}^k(-1{)}^j(\genfrac{}{}{0ex}{}{\alpha }{j})f(t-jh),h=\frac{t-a}{k}$$

$p.s.$ 其实也就是上面说过的$Grunwald-Letnikov$ 公式，就不再赘述啦~

关键是下面的——

基于 $Gr\ddot{u}nwald-Letnikov$ 公式的数值方法可以将分数阶微分方程转化为一阶常微分方程组，然后用常用的数值方法如 $Euler$ 法，$Runge-Kutta$ 法等求解。具体的算法步骤如下：

1. 给定分数阶微分方程 ${}_a{D}_t^{\alpha }y(t)=f(t,y(t))$，初始条件 $y(a)=y_0$，分数阶导数的阶数 $\alpha \in (0,1)$，求解区间 $[a,b]$，步长 $h$，步数 $k$。
2. 用 $Gr\ddot{u}nwald-Letnikov$ 公式将分数阶微分方程转化为一阶常微分方程组¹³：

$$\{\begin{array}{l}{y}_0^{\prime }(t)=\frac{1}{h^{\alpha }}\left[f(t,y_0(t))-\left(\genfrac{}{}{0ex}{}{\alpha }{1}\right)y_0(t)\right]\\ y_j^{\prime }(t)=\frac{1}{h^{\alpha }}\left[y_{j-1}(t)-\left(\genfrac{}{}{0ex}{}{\alpha }{j+1}\right)y_j(t)\right],j=1,2,\dots ,k-1\\ y_k^{\prime }(t)=\frac{1}{h^{\alpha }}{y}_{k-1}(t)\end{array}$$

• 其中 $y_0(t)=y(t)$ 是原方程的解，$y_j(t)$ 是辅助变量，初始条件为 $y_j(a)=0$，$j=1,2,\dots ,k$。

3. 用数值方法如 $Euler$ 法，$Runge-Kutta$ 法等求解一阶常微分方程组，得到 $y_0(t)$ 的近似解，即 $y(t)$ 的近似解。

4.3 一起来实战吧

以一个例题为例，假设要求解如下的分数阶微分方程¹⁴：

$${}_0{D}_t^{0.5}y(t)=-y(t)+t,t\in [0,10],y(0)=1$$

其中 $\alpha =0.5$，$f(t,y(t))=-y(t)+t$，$y_0=1$，$a=0$，$b=10$。我们取步长 $h=0.01$，步数 $k=100$，用四阶 Runge-Kutta 法求解一阶常微分方程组。

数理过程我就不写啦，就说说具体是怎么代码实现的吧~

MATLAB 代码如下：

% Define the parameters
alpha = 0.5; % fractional order
f = @(t,y) -y + t; % right hand side function
y0 = 1; % initial condition
a = 0; % left endpoint
b = 10; % right endpoint
h = 0.01; % step size
k = 100; % step number

% Define the binomial coefficients
binom = zeros(1,k+1);
binom(1) = 1;
for j = 1:k
    binom(j+1) = binom(j) * (alpha - j + 1) / j;
end

% Define the Runge-Kutta function
function [t,y] = rk4(f,a,b,h,y0)
    % Input: f is the function handle
    % a and b are the left and right endpoints
    % h is the step size
    % y0 is the initial condition
    % Output: t is the vector of time points
    % y is the vector of solution values
    t = a:h:b; % create the time vector
    n = length(t); % get the number of steps
    y = zeros(1,n); % initialize the solution vector
    y(1) = y0; % set the initial condition
    for i = 1:n-1 % loop over the steps
        % calculate the slopes
        k1 = f(t(i),y(i));
        k2 = f(t(i)+h/2,y(i)+h*k1/2);
        k3 = f(t(i)+h/2,y(i)+h*k2/2);
        k4 = f(t(i)+h,y(i)+h*k3);
        % update the solution
        y(i+1) = y(i) + h/6 * (k1 + 2*k2 + 2*k3 + k4);
    end
end

% Define the auxiliary variables
y = zeros(k+1,1); % initialize the solution vector
y(1) = y0; % set the initial condition
z = zeros(k,1); % initialize the auxiliary vector
z(1) = y0; % set the initial condition

% Solve the fractional differential equation
for i = 1:k % loop over the steps
    % solve the first equation using Runge-Kutta method
    [t,y] = rk4(@(t,y) 1/h^alpha * (f(t,y) - binom(2) * y),a,a+h,y(1));
    % update the solution vector
    y(i+1) = y(end);
    % solve the other equations using Runge-Kutta method
    for j = 2:k
        [t,z] = rk4(@(t,z) 1/h^alpha * (y(j-1) - binom(j+1) * z(j)),a,a+h,z(j));
        % update the auxiliary vector
        z(j) = z(end);
    end
end

% Plot the solution
plot(0:h:b,y,'b')
xlabel('t')
ylabel('y(t)')
title('Solution of the fractional differential equation')
grid on

然后是C++ 代码：

// Include the necessary libraries
#include <iostream>
#include <cmath>
#include <vector>
#include <limits>
using namespace std;

// Define the parameters
const double alpha = 0.5; // fractional order
const double y0 = 1.0; // initial condition
const double a = 0.0; // left endpoint
const double b = 10.0; // right endpoint
const double h = 0.01; // step size
const int k = 100; // step number

// Define the right hand side function
double f(double t, double y) {
    return -y + t;
}

// Define the Runge-Kutta function
void rk4(double (*f)(double, double), double a, double b, double h, double y0, vector<double>& t, vector<double>& y) {
    // Input: f is the function pointer
    // a and b are the left and right endpoints
    // h is the step size
    // y0 is the initial condition
    // Output: t is the vector of time points
    // y is the vector of solution values
    int n = (b - a) / h; // get the number of steps
    t.resize(n+1); // resize the time vector
    y.resize(n+1); // resize the solution vector
    t[0] = a; // set the initial time
    y[0] = y0; // set the initial condition
    for (int i = 0; i < n; i++) { // loop over the steps
        // calculate the slopes
        double k1 = f(t[i], y[i]);
        double k2 = f(t[i] + h/2, y[i] + h*k1/2);
        double k3 = f(t[i] + h/2, y[i] + h*k2/2);
        double k4 = f(t[i] + h, y[i] + h*k3);
        // update the solution
        y[i+1] = y[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4);
        // update the time
        t[i+1] = t[i] + h;
    }
}

// Define the main function
int main() {
    // Define the binomial coefficients
    vector<double> binom(k+1);
    binom[0] = 1;
    for (int j = 1; j <= k; j++) {
        binom[j] = binom[j-1] * (alpha - j + 1) / j;
    }

    // Define the auxiliary variables
    vector<double> y(k+1); // initialize the solution vector
    y[0] = y0; // set the initial condition
    vector<double> z(k); // initialize the auxiliary vector
    z[0] = y0; // set the initial condition

    // Solve the fractional differential equation
    for (int i = 0; i < k; i++) { // loop over the steps
        // solve the first equation using Runge-Kutta method
        vector<double> t, y_temp;
        rk4([](double t, double y) {return 1/pow(h,alpha) * (f(t,y) - binom[1] * y);}, a, a+h, h, y[i], t, y_temp);
        // update the solution vector
        y[i+1] = y_temp.back();
        // solve the other equations using Runge-Kutta method
        for (int j = 1; j < k; j++) {
            vector<double> t, z_temp;
            rk4([=](double t, double z) {return 1/pow(h,alpha) * (y[j-1] - binom[j+1] * z);}, a, a+h, h, z[j], t, z_temp);
            // update the auxiliary vector
            z[j] = z_temp.back();
        }
    }

    // Print the solution
    for (int i = 0; i <= k; i++) {
        cout << "y(" << i*h << ") = " << y[i] << endl;
    }

    return 0;
}

尾声

本文讲了好多好多干货，可能一时间消化不了。建议抽空的时候反复看，多消化。

当你看懂之后，也许，会浮现出一年前，在那个汗水被冷风凝固的冬夜，你聚精会神地听着张宇的闭关修炼时，期盼现在的你一个月后折桂的一路生花吧~

祝君成功上岸，功不唐捐！

参考

---

1. ^Fractional calculus.” Wikipedia, 6 Jan. 2024. ↩︎

2. ^Weisstein, Eric W. “Fractional Derivative.” MathWorld, Wolfram Research. ↩︎

3. ^Flammable Maths. “The Fractional Derivative, what is it? | Introduction to Fractional Calculus.” YouTube, 20 Jan. 2018. ↩︎

4. ^Almeida, Ricardo, et al. “Fractional derivatives and the fundamental theorem of fractional calculus.” Fractional Calculus and Applied Analysis, vol. 23, no. 6, 2020, pp. 1634-1652. ↩︎

5. ^Li, Changpin, et al. “FdeSolver: A Julia Package for Solving Fractional Differential Equations.” arXiv, Cornell University, 28 Dec. 2021. ↩︎

6. ^user1709. “Applications of Fractional Calculus.” Mathematics Stack Exchange, 30 Sept. 2010. ↩︎

7. ^Almeida, Ricardo, et al. “A Generalized Definition of the Fractional Derivative with Applications.” Mathematical Problems in Engineering, vol. 2021, 2021, ↩︎

8. ^The Organic Chemistry Tutor. “Limit Definition of Derivative Square Root, Fractions, 1/sqrt(x), Examples - Calculus.” YouTube, 17 May 2017. ↩︎

9. ^Erhan, Ozan. “Fractional Derivatives. And How to Calculate Them.” Medium, 10 Feb. 2019. ↩︎

10. ^Tenreiro Machado, J. A., et al. Fractional Derivative Modeling in Mechanics and Engineering. Springer, 2021. ↩︎

11. ^Sen, Mihir. “Introduction to Fractional-Order Operators and Their Engineering Applications.” University of Notre Dame. ↩︎

12. ^Podlubny, Igor, et al. “Selected Engineering Applications of Fractional-Order Calculus.” In Fractional Calculus: Theory, Applications, and Numerics, edited by Igor Podlubny and J. A. Tenreiro Machado, 323-358. Springer, 2021. ↩︎

13. ^Li, Changpin, et al. “FdeSolver: A Julia Package for Solving Fractional Differential Equations.” arXiv, Cornell University, 28 Dec. 2021. ↩︎

14. ^user1559695. “Implementing the derivative in C/C++.” Stack Overflow, 11 Oct. 2009, ↩︎