1

我正在尝试使用来自DifferentialEquations.jl. 我想我可以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找到的示例包含回调函数的两个不同的函数签名，但我对两者都有相同的问题。我也不清楚回调是应该修改in 的值，还是返回新值。z

编辑 2： 基于 Chris Rackauckas 的回答，并查看此示例，我已通过反射函数进行了修改：

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. 在第一次积分器调用后存储输出之前未调用回调。

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

4

1 回答 1

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)
于 2019-03-05T21:09:59.817 回答