I’ve been trying to get the hang of solving DEs in Julia (i.e. new user). As a simple task, I tried to set up a 1-d Schrödinger equation and perform a time evolution. Here’s my effort so far:
using ModelingToolkit
using MethodOfLines
using DomainSets
@variables t [bounds = (0.0, Inf)]
@variables x [bounds = (-5.0, 5.0)]
@variables ur(..), ui(..)
Dx = Differential(x)
Dt = Differential(t)
equationri = [
Dt(ur(t,x)) ~ -(1/2)*(Dx^2)(ui(t,x)), # real LHS
Dt(ui(t,x)) ~ (1/2)*(Dx^2)(ur(t,x)), # imag LHS
]
domains=[
t ∈ (0.0, Inf),
x ∈ (-5.0,5.0)
]
boundary_conditions=[
ur(0.0,x) ~ cos(x) * exp(-x^2/4) / (2π)^(1/4),
ui(0.0,x) ~ sin(x) * exp(-x^2/4) / (2π)^(1/4),
ur(t,-5.0) ~ 0.0,
ui(t,-5.0) ~ 0.0,
ur(t,5.0) ~ 0.0,
ui(t,5.0) ~ 0.0,
]
@named schro_system = PDESystem(
equationri,
boundary_conditions,
domains,
[t, x],
[ur(t,x), ui(t,x)]
)
dx=0.1
discretization = MOLFiniteDifference([x => dx], t; approx_order = 2)
problem = discretize(schro_system, discretization)
using OrdinaryDiffEq
st=init(problem, Tsit5()) # here we are initialized
# step with `step!(st, 0.05, true)` (don't know why kwarg doesn't work)
# Solution in st.sol
The problem I’m having is that I can not get the solution out of the st.sol
ODESolution object. It seems to suggest that it is in the index u
, but “u” is undefined. Doing an index of ur(t,x)
just gives
ERROR: ArgumentError: ur(t, x) is neither an observed nor a state variable.
which is a somewhat confusing message.
The documentation here seems somewhat lacking. There is some warning about the “old” system for discretising. But the “new” system doesn’t actually seem to work.
So, I guess my questions are: is the above code snippet correct usage? If so, how can I get the solution read out after perfoming forwards integration step?