As a basis for more complex models, I’m trying to get a simplified version of the Von Foerster equation - Wikipedia discretized and running in MethodOfLines.jl.
Here’s my model:
using ModelingToolkit
using MethodOfLines
using OrdinaryDiffEq
using DomainSets
@parameters t a μ
@variables v(..) u(..)
Dt = Differential(t)
Da = Differential(a)
p = [μ => 0.04]
tmax = 100.0
amax = 100.0
da = 0.1
dt = 0.1
pde_eqns = [Dt(u(t,a)) + Da(u(t,a)) ~ -μ*u(t,a)]
domains = [t ∈ Interval(0,tmax), a ∈ Interval(0,amax)]
function a0(a)
if a > 0
return 0.0
else
return 1.0
end
end
@register a0(a)
bcs = [u(t,0) ~ 0.0, u(0,a) ~ a0(a)]
@named pde_sys = PDESystem(pde_eqns, bcs, domains, [t,a], [u(t,a)], p)
discretization = MOLFiniteDifference([a=>da], t, approx_order = 2, grid_align = center_align)
pde_prob = discretize(pde_sys, discretization)
pde_sol = solve(pde_prob, Tsit5(), saveat = dt)
The boundary conditions for u(0,a) are meant to be 1 when a=0 and 0 when a>0, which is what I’m trying to do with the registered function. However, when I look at the discretized problem, the initial conditions are 0 for all values of a:
julia> pde_prob.u0
1000-element Vector{Float64}:
0.0
0.0
0.0
0.0
0.0
⋮
0.0
0.0
0.0
0.0
0.0
I think I understand what’s going on here - while the above function a0 is correct in the continuous sense, it doesn’t work in the discretized system. If I change function:
function a0(a)
if a > da
return 0.0
else
return 1.0
end
end
This now does something more sensible:
julia> pde_prob.u0
1000-element Vector{Float64}:
1.0
0.0
0.0
0.0
0.0
⋮
0.0
0.0
0.0
0.0
0.0
However, now the solution gives nonsensical results - the integral over a should be an exponential decay, but it shows exponential growth, which again points to another boundary issue problem. Any tips on how to fix this? Full code, along with the ‘correct’ exponential decay as an ODE is given below.
using ModelingToolkit
using MethodOfLines
using OrdinaryDiffEq
using DomainSets
@parameters t a μ
@variables v(..) u(..)
Dt = Differential(t)
Da = Differential(a)
p = [μ => 0.04]
tmax = 100.0
amax = 100.0
da = 0.1
dt = 0.1
# ODE version
ode_eqns = [Dt(v(t)) ~ -μ*v(t)]
v0 = [v(t) => 1.0]
@named ode_sys = ODESystem(ode_eqns)
ode_prob = ODEProblem(ode_sys, v0, (0,tmax), p)
ode_sol = solve(ode_prob, Tsit5(), saveat = dt)
# PDE version
pde_eqns = [Dt(u(t,a)) + Da(u(t,a)) ~ -μ*u(t,a)]
domains = [t ∈ Interval(0,tmax), a ∈ Interval(0,amax)]
function a0(a)
if a > da
return 0.0
else
return 1.0
end
end
@register a0(a)
bcs = [u(t,0) ~ 0.0, u(0,a) ~ a0(a)]
@named pde_sys = PDESystem(pde_eqns, bcs, domains, [t,a], [u(t,a)], p)
discretization = MOLFiniteDifference([a=>da], t, approx_order = 2, grid_align = center_align)
pde_prob = discretize(pde_sys, discretization)
pde_sol = solve(pde_prob, Tsit5(), saveat = dt)
