How to get PDE solutions when doing foward time steps

Specific Domains Modelling & Simulations
1

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?

2

Yes, the initialising the solution it isn’t yet wrapped, you have to use the old method. Are you sure you don’t want to just solve the problem? Init is called internally in the solve method so you don’t need to do it manually, what are you trying to do?

3

I have been trying to do an animated plot with Makie of the solution in 1-d. I’ve been fiddling with this and eventually found using the “old method” works. I was initially trying to only integrate up to the time step in which a single frame would be drawn.