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 的值，还是返回新值。

编辑 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

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

1. 第一个输出值就是z0。我希望首先应用反射，因为我设置了func_start = true。
2. 第二个值也是负数表示两件事：
   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)