使用Julia的Quadgk计算Caputo分数阶导数时遇定义域错误

问题根源 #

当计算t=3时，被积函数g(x)/(γ(1/2)*√(t-x))在x趋近于3时，分母√(3-x)趋近于0，产生可积奇点。QuadGK的自适应积分采样到接近端点的点时，直接计算出Inf，触发定义域错误，但这个奇点本身是可积的，只是默认积分逻辑没适配这类端点奇异性。

解决方案：变量替换消除端点奇点 #

通过变量替换将端点奇点转化为常规积分区间，具体推导：
令x = t - t*s²，则dx = -2t*s ds，当x=0时s=1，x=t时s=0，积分上下限反转后抵消负号：
$$
\int_0^t \frac{g(x)}{\sqrt{t-x}} dx = \int_0^1 \frac{g(t - t s²)}{\sqrt{t s²}} \cdot 2t s ds = 2\sqrt{t} \int_0^1 g(t(1-s²)) ds
$$
这样就消除了分母的奇点，转化为常规的[0,1]区间积分。

修改后的代码 #

using QuadGK, ForwardDiff, SpecialFunctions

f(t) = t^2 * sin(t)
g(t) = ForwardDiff.derivative(f, t)  # 简化写法，直接传入原函数更清晰

# 基于变量替换实现Caputo导数计算
caputo(t) = 2 * sqrt(t) / gamma(1/2) * quadgk(s -> g(t*(1 - s^2)), 0, 1)[1]

# 测试t=3的情况
println(caputo(3))  # 输出接近-3.11522

结果验证 #

运行上述代码，输出结果约为-3.11522，和Mathematica的计算结果一致，且不会触发定义域错误。

内容的提问来源于stack exchange，提问作者Hariz Khaled