Julia 中带回调的随机微分方程

Question

我正在尝试使用 DifferentialEquations.jl 中的各种 SDE 积分器来解决具有反射边界的扩散问题。我想我可以使用 FunctionCallingCallback 来处理边界，通过在每个积分器步骤之后反映关于域边界的解决方案。

这是我的代码

using DifferentialEquations

K0 = 1e-3
K1 = 5e-3 
alpha = 0.5

K(z)    = K0 + z*K1*exp(-z*alpha)
dKdz(z) = K1*exp(-alpha*z) - K1*alpha*z*exp(-alpha*z)

a(z,p,t) = dKdz(z)
b(z,p,t) = sqrt(2*K(z))

dt = 0.1
tspan = (0.0,1.0)

z0 = 1.0

prob  = SDEProblem(a,b,z0,tspan)

function reflect(z, p, t, integrator)
    bottom = 2.0
    if z < 0
        # Reflect from surface
        z = -z
    elseif z > bottom
        # Reflect from bottom
        z = 2*bottom - z
    end
    return z
end

cb = FunctionCallingCallback(reflect;
                 func_everystep = true,
                 func_start = true,
                 tdir=1)

sol = solve(prob, EM(), dt = dt, callback = cb)

编辑： 由于 Chris Rackauckas 的评论解决了我最初的问题后，我修改了我的反射函数。现在代码运行了，但是解决方案包含负值，应该在每一步之后通过大约 0 的反射来防止它。

任何关于这里出了什么问题的想法都将不胜感激。

请注意，FunctionCallingCallback 示例发现 here 包含两个不同的回调函数签名，但我遇到了同样的问题。我也不清楚回调是否应该修改 z 的值，或者 return 新值。

编辑 2： 根据 Chris Rackauckas 的回答，并查看 this example，我通过反射函数进行了修改：

function reflect(z, t, integrator)
    bottom = 2.0
    if integrator.u < 0
        # Reflect from surface
        integrator.u = -integrator.u
    elseif integrator.u > bottom
        # Reflect from bottom
        integrator.u = 2*bottom - integrator.u
    end
    # Not sure if the return statement is required
    return integrator.u
end

运行这与初始条件 z0 = -0.1 产生以下输出：

retcode: Success
Interpolation: 1st order linear
t: 11-element Array{Float64,1}:
 0.0                
 0.1                
 0.2                
 0.30000000000000004
 0.4                
 0.5                
 0.6                
 0.7                
 0.7999999999999999 
 0.8999999999999999 
 1.0                
u: 11-element Array{Float64,1}:
 -0.1                 
 -0.08855333388147717 
  0.09862543518953905 
  0.09412012313587219 
  0.11409372573454478 
  0.10316400521980074 
  0.06491042188420941 
  0.045042097789392624
  0.040565317051189105
  0.06787136817395374 
  0.055880083559589955

在我看来，这里发生的事情是：

第一个输出值刚好z0。我希望首先应用反射，因为我设置了 func_start = true.
第二个值也是负值表示两件事：
1. 在第一次集成商调用之前未调用回调。
2. 在第一次集成器调用后存储输出之前未调用回调。

我本以为输出中的所有值都是正数（即，在存储输出之前对它们应用回调）。我是做错了什么，还是应该简单地调整一下我的期望？

Answer 1

FunctionCallingCallback 是一个函数 (u,t,integrator)，所以我不确定您的代码如何没有出错。应该是：

using DifferentialEquations

K0 = 1e-3
K1 = 5e-3
alpha = 0.5

K(z)    = K0 + z*K1*exp(-z*alpha)
dKdz(z) = K1*exp(-alpha*z) - K1*alpha*z*exp(-alpha*z)

a(z,p,t) = dKdz(z)
b(z,p,t) = sqrt(2*K(z))

dt = 0.1
tspan = (0.0,1.0)

z0 = 1.0

prob  = SDEProblem(a,b,z0,tspan)

function reflect(z, t, integrator)
    bottom = 2.0
    if z < 0
        # Reflect from surface
        z = -z
    elseif z > bottom
        # Reflect from bottom
        z = 2*bottom - z
    end
    return z
end

cb = FunctionCallingCallback(reflect;
                 func_everystep = true,
                 func_start = true,
                 tdir=1)

sol = solve(prob, EM(), dt = dt, callback = cb)

编辑

您不需要函数调用回调。只需使用正常回调：

using DifferentialEquations

K0 = 1e-3
K1 = 5e-3
alpha = 0.5

K(z)    = K0 + z*K1*exp(-z*alpha)
dKdz(z) = K1*exp(-alpha*z) - K1*alpha*z*exp(-alpha*z)

a(z,p,t) = dKdz(z)
b(z,p,t) = sqrt(2*K(z))

dt = 0.1
tspan = (0.0,1.0)

z0 = 1.0

prob  = SDEProblem(a,b,z0,tspan)

condition(u,t,integrator) = true
function affect!(integrator)
    bottom = 2.0
    if integrator.u < 0
        # Reflect from surface
        integrator.u = -integrator.u
    elseif integrator.u > bottom
        # Reflect from bottom
        integrator.u = 2*bottom - integrator.u
    end
end

cb = DiscreteCallback(condition,affect!;save_positions=(false,false))

sol = solve(prob, EM(), dt = dt, callback = cb)

Julia 中带回调的随机微分方程

Stochastic differential equation with callback in Julia

callback

julia

differentialequations.jl

编辑