如何在 Julia 中求解随机微分方程？

Question

我试图了解如何用数值方法求解随机微分方程 (SDE)（我没有使用任何语言的经验，但出于某些原因我选择了 Julia）。作为起始模型，我决定使用 Lotka-Volterra equations. I read manual and tutorial for DifferentialEquations.jl。虽然我的模型是一个简单的 ODE 系统：

一切正常：

using DifferentialEquations
using Plots
function lotka(du,u,p,t);
    , , ,  = p; 
    du[1] = *u[1]-*u[1]*u[2];
    du[2] = *u[1]*u[2]-*u[2];
end
u0 = [15.0,1.0,0.0,0.0];
p = (0.3,0.05,0.7,0.1);
tspan = (0.0,50.0);
prob = ODEProblem(lotka,u0,tspan,p);
sol = solve(prob);
plot(sol,vars = (1,2))

现在我想添加 Ornstein-Uhlenbeck 噪声：

愚蠢直接的解决方案是：

using DifferentialEquations
using Plots
function lotka(du,u,p,t);
    , , , , ,  = p; 
    du[1] = *u[1]-*u[1]*u[2]+u[3];
    du[2] = *u[1]*u[2]-*u[2]+u[4];
    du[3] = -u[3]/+sqrt((2.0*^2.0/))*randn();
    du[4] =  -u[4]/+sqrt((2.0*^2.0/))*randn();
end
u0 = [15.0,1.0,0.0,0.0];
p = (0.3,0.05,0.7,0.1,0.1,0.1);
tspan = (0.0,50.0);
prob = ODEProblem(lotka,u0,tspan,p);
sol = solve(prob);

但是，正如预期的那样，它没有工作，因为求解器不适用于此类 SDE 问题。

┌ Warning: Interrupted. Larger maxiters is needed.
└ @ DiffEqBase /home/jrun/.julia/packages/DiffEqBase/ujEgN/src/integrator_interface.jl:116

我试着阅读 SDE Julia documentation，但没有很好的例子，我无法理解如何处理它。而且，我的数学基础很差，似乎我没有正确理解符号。我怎样才能使它适用于 SDE？

更新：

最后，我有以下代码：

using DifferentialEquations,Plots;

function lotka(du,u,p,t);
    , , , , ,  = p; 
    du[1] = *u[1]-*u[1]*u[2];
    du[2] = *u[1]*u[2]-*u[2];
end
function stoch(du,u,p,t)
  , , , , ,  = p; 
  du[1] = 1 
  du[2] = 1
end
u0 = [15.0,1.0];
p = (0.3,0.05,0.7,0.1,0.9,0.1);
, , , , ,  = p; 
OU = OrnsteinUhlenbeckProcess(1/, 0.0, , 0.0, 0.0);
tspan = (0.0,100.0);
prob = SDEProblem(lotka,stoch,u0,tspan,p,noise = OU);
sol = solve(prob,EM(),dt = 1e-3, adaptive=false);

为了结束这个话题，我还有最后两个问题：

代码正确吗？（恐怕，在这种情况下我的噪声不是对角线的）
我可以为两个方程中的每一个设置不同的初始噪声 (W0) 吗？

Answer 1

您似乎有一个 SDE，其中前两项由随机的后两项驱动？在这种情况下，您可以使用 lotka 确定性术语：

function lotka(du,u,p,t);
    , , , , ,  = p; 
    du[1] = *u[1]-*u[1]*u[2]+u[3];
    du[2] = *u[1]*u[2]-*u[2]+u[4];
    du[3] = -u[3]/
    du[4] = -u[4]/
end

然后是stoch随机部分：

function stoch(du,u,p,t)
  , , , , ,  = p; 
  du[1] = 0 
  du[2] = 0
  du[3] = sqrt((2.0*^2.0/))
  du[4] = sqrt((2.0*^2.0/))
end

现在写成du = f(u,p,t)dt + g(u,p,t)dW的形式。请注意，就像您不写 dt 一样，您也不会写 randn()，因为它已在求解器中处理（相当复杂）。使用这些定义 SDEProblem 并求解：

u0 = [15.0,1.0,0.0,0.0];
p = (0.3,0.05,0.7,0.1,0.1,0.1);
tspan = (0.0,50.0);
prob = SDEProblem(lotka,stoch,u0,tspan,p);
sol = solve(prob);

（你确实需要小心这个模型，因为它不能保证保持积极，所以如果有太多噪音它可能会变得不稳定。不仅仅是积分器，还有解决方案本身）

为什么使用 ODE 求解器不起作用的快速解释

为清楚起见，您在上面所做的操作不起作用的原因有两个。一方面，自适应性根据解决方案属性做出假设。例如，标准的 5 阶积分器假设解是 5 次可微分的，并在其误差估计中使用它来选择步长。 SDE 是 0 倍可微的：对于每个 epsilon<1/2，它只有 epsilon 可微。因此，当直接对 SDE 使用 ODE 方法时，误差估计和时间步长的选择将会大大偏离。

但其次，ODE 求解器使用拒绝采样进行调整。他们将选择一个 dt，试一试，如果失败，则减少 dt。在这里，你在你的导数中取了一个随机数，5 阶积分器将调用你的导数函数 7 次，得到它认为导数的 7 个不同的值，计算误差估计，意识到它非常偏离（因为它甚至不是可微所以整个算法做出了错误的假设），减少时间步长，然后采用完全不同的随机数，使得较小的 dt 最终成为完全不同的轨迹。如您所知，这整个事情非常不稳定，实际上并没有解决我们最初拥有的真正 SDE。

您可以通过更聪明地处理上述随机数、定义误差估计以及使用假定 "Ito differentiability"（即 "differentiability" 在SDE 组件的条款）。这被描述为 in the paper which is the foundation of the current crop of adaptive SDE solvers。

回答更新

对于您拥有的代码

using DifferentialEquations,Plots;

function lotka(du,u,p,t);
    , , , , ,  = p; 
    du[1] = *u[1]-*u[1]*u[2];
    du[2] = *u[1]*u[2]-*u[2];
end
function stoch(du,u,p,t)
  , , , , ,  = p; 
  du[1] = 1 
  du[2] = 1
end
u0 = [15.0,1.0];
p = (0.3,0.05,0.7,0.1,0.9,0.1);
, , , , ,  = p; 
OU = OrnsteinUhlenbeckProcess(1/, 0.0, , 0.0, 0.0);
tspan = (0.0,100.0);
prob = SDEProblem(lotka,stoch,u0,tspan,p,noise = OU);
sol = solve(prob,EM(),dt = 1e-3, adaptive=false);

这两个变量都添加了 1 个 OU 进程。如果您希望它们是对角线的，即两个独立的 OU 进程，您使用：

OU = OrnsteinUhlenbeckProcess(1/, 0.0, , 0.0, [0.0,0.0]);

更一般地说，[0.0,0.0] 是起点，您更改这些以解决 (2)。对于小的性能优化，您可以选择：

OU = OrnsteinUhlenbeckProcess!(1/, 0.0, , 0.0, [0.0,0.0]);

这个公式和使用布朗 SDE 的公式都有效。 SDE 公式适用于自适应方法，但系统较大，而 OUProcess 公式是较小的系统，但仅适用于 EM() 和固定时间步长。哪个更好将取决于问题，但在大多数情况下我可能更喜欢 SDE。不过，OUProcess 表单在 RODE 上会好得多。

如何在 Julia 中求解随机微分方程？

How do I solve stochastic differential equations in Julia?

julia

numerical-methods

differential-equations

stochastic

更新：

为什么使用 ODE 求解器不起作用的快速解释

回答更新