Solving FPE with boundary condition at t=T rather than t=0

[jdrugo]

First of all, thank you for creating such a nice set of tools for solving FPEs and ODEs!

I ran into an issue when trying to solve an FPE (arising from an HJB equation) that has a boundary condition at t=T rather than at t=0. The FPE has the known solution, given in the code by trueJ, but I wanted to see how well Julia's MOL approximate this solution (as preparation for solving more complex problems). Unfortunately, I don't seem to get the following code to work:

using ModelingToolkit, MethodOfLines
import ModelingToolkit: Interval, infimum, supremum

# define the problem with indep. variables x & t and dependent variable J
@parameters x t
@variables J(..)
Dx = Differential(x)
Dt = Differential(t)
# parameters
T=5.0
α=1.0
# FPE boundaries
X=5.0
trueJ(x, t) = 0.5*x^2 / (T - t + 1/α)
# problem equations
eqs = [ Dt(J(x,t)) ~ 0.5*Dx(J(x,t))^2 ]
bcs = [ J(x,T) ~ trueJ(x, T),
        J(X,t) ~ trueJ(X, t),
        J(-X,t) ~ trueJ(-X, t)]
domains = [ t ∈ Interval(0.0, T),
            x ∈ Interval(-X, X)]
@named pde_system = PDESystem(eqs, bcs, domains, [x, t], [J(x, t)])

# discretized system
N = 4
dx = 2*X/N
order = 2
discretization = MOLFiniteDifference([x=>dx], t, approx_order=order, grid_align=center_align)
@time prob = discretize(pde_system, discretization)

Running the above yields

The system of equations is:
Equation[Differential(t)(J[2](t)) ~ 0.5((0.3J[3](t) - 0.3J[2](t))^2), Differential(t)(J[3](t)) ~ 0.5((0.3J[4](t) - 0.3J[3](t))^2), J[4](t) ~ 12.5 / (6.0 - t), J[1](t) ~ 12.5 / (6.0 - t)][Base.OneTo(4)]

Discretization failed, please post an issue on https://github.com/SciML/MethodOfLines.jl with the failing code and system at low point count.

Term{Real, Base.ImmutableDict{DataType, Any}}[J[1](t), J[2](t), J[3](t), J[4](t)][Base.OneTo(4)]

ArgumentError: Term{Real, Base.ImmutableDict{DataType, Any}}[J[2](t), J[3](t)] are missing from the variable map.

Stacktrace:
  [1] throw_missingvars(vars::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/variables.jl:89
  [2] varmap_to_vars(varmap::Vector{Pair}, varlist::Vector{Term{Real, Base.ImmutableDict{DataType, Any}}}; defaults::Dict{Any, Any}, check::Bool, toterm::Function, promotetoconcrete::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/variables.jl:42
  [3] process_DEProblem(constructor::Type, sys::ODESystem, u0map::Vector{Pair}, parammap::SciMLBase.NullParameters; implicit_dae::Bool, du0map::Nothing, version::Nothing, tgrad::Bool, jac::Bool, checkbounds::Bool, sparse::Bool, simplify::Bool, linenumbers::Bool, parallel::Symbolics.SerialForm, eval_expression::Bool, kwargs::Base.Pairs{Symbol, Bool, Tuple{Symbol}, NamedTuple{(:has_difference,), Tuple{Bool}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:566
  [4] (ODEProblem{true})(sys::ODESystem, u0map::Vector{Pair}, tspan::Tuple{Float64, Float64}, parammap::SciMLBase.NullParameters; callback::Nothing, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:657
  [5] (ODEProblem{true})(sys::ODESystem, u0map::Vector{Pair}, tspan::Tuple{Float64, Float64}, parammap::SciMLBase.NullParameters) (repeats 2 times)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:656
  [6] ODEProblem(::ODESystem, ::Vector{Pair}, ::Vararg{Any}; kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:635
  [7] ODEProblem(::ODESystem, ::Vector{Pair}, ::Vararg{Any})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/a7NhS/src/systems/diffeqs/abstractodesystem.jl:635
  [8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{MethodOfLines.CenterAlignedGrid})
    @ MethodOfLines ~/.julia/packages/MethodOfLines/9MLPx/src/discretization/MOL_discretization.jl:152
  [9] top-level scope
    @ ./timing.jl:220 [inlined]
 [10] top-level scope
    @ ./In[86]:0
 [11] eval
    @ ./boot.jl:373 [inlined]
 [12] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1196

As you can see, I have boundary conditions on some chosen x range, as well as on t=T. That's the natural thing to do for the problem I am interested in. If I instead change the boundary condition to t=0 by using

bcs = [ J(x,0.0) ~ trueJ(x, 0.0),
        J(X,t) ~ trueJ(X, t),
        J(-X,t) ~ trueJ(-X, t)]

the code seems to work without problems. Are boundary conditions on t=T simply not supported? Or is something else going on here?

[xtalax]

Unfortunately boundary conditions on the timespan supremum are not supported, though you can get around this with a change of variables -t => s.

[jdrugo]

Thank you for the reply. I will use a change of variables instead. It would be good to add this as a note to the documentation, maybe at http://methodoflines.sciml.ai/dev/#limitations.