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. 第二个值也是负数表示两件事： 在第一次积分器调用之前未调用回调。 在第一次积分器调用之后存储输出之前未调用回调。

我希望输出中的所有值都是正的（即，在存储输出之前应用回调）。我做错了什么，还是应该调整我的期望？

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)