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

Question

我尝试了解如何以数值方式求解随机微分方程 (SDE)（我没有任何语言的经验，但出于某些原因我选择了 Julia）。作为初始模型，我决定使用Lotka-Volterra equations。我为DifferentialEquations.jl 阅读了手册和tutorial。虽然我的模型是一个简单的 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);

要结束这个话题，我还有两个问题：

1. 它是正确的代码吗？ （恐怕在这种情况下我的噪音不是对角线）

2. 我可以为两个方程中的每一个设置不同的初始噪声 (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

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

您可以通过更聪明地处理上述随机数、如何定义误差估计以及使用假设“Ito 可微性”（即 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 上会好得多。