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

Question

我尝试了解如何用数值方法求解随机微分方程（SDE）（我没有任何语言的经验，但出于某些原因，我选择了Julia）。作为初始模型，我决定使用Lotka-Volterra方程。我阅读了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文档，但是如果没有很好的例子，我将无法理解如何处理它。此外，我的数学背景很差，似乎我无法正确理解表示法。如何使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小于1/2时，它都是可微的。因此，当直接将ODE方法用于SDE时，错误估计和时间步长选择将大相径庭。

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

您可以通过更聪明地掌握所说的随机数，如何定义误差估计以及使用假设“ Ito可微性”（即，根据SDE组件的“可微性”）的高阶方法来解决此问题。本文对此进行了描述，这是当前自适应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 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上会更好。