Стохастическое дифференциальное уравнение с обратным вызовом по Юлии - PullRequest
Новая
       71

Стохастическое дифференциальное уравнение с обратным вызовом по Юлии

1 голос
/ 05 марта 2019

Я пытаюсь решить диффузионную проблему с отражением границ, используя различные интеграторы SDE из 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)

Редактировать: После решения моей первоначальной проблемы благодаря комментарию Криса Ракауцкаса я изменил свою функцию отражения. Теперь код выполняется, но решение содержит отрицательные значения, которые должны были предотвращаться отражением около 0 после каждого шага.

Будем весьма благодарны за любые идеи относительно того, что здесь происходит не так.

Обратите внимание на то, что найденный пример FunctionCallingCallback здесь содержит две разные сигнатуры функций для функции обратного вызова, но у меня одна и та же проблема с обеими. Мне также не ясно, должен ли обратный вызов изменить значение z на месте или вернуть новое значение.

Редактировать 2: Основываясь на ответе Криса Ракауцкаса и глядя на этот пример , я изменил функцию рефлекса следующим образом: 

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. Обратный вызов не был вызван до сохранения результата после первого вызова интегратора.

Я ожидал бы, что все значения в выходных данных будут положительными (то есть обратный вызов будет применен к ним перед сохранением выходных данных). Я делаю что-то не так или мне просто нужно скорректировать свои ожидания?

1 Ответ

1 голос
/ 06 марта 2019

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)
...