constant norm condition

[berceanu]

I am working on a simple quantum mechanics problem for benchmarking the new ode solver. It involves solving the discretized Schrodinger equation on a 2D grid with a harmonic potential and looking for the steady state. Since the analytical solution is known, I thought it would be a nice test, and might also prove instructive as an example. I am solving the time propagation in imaginary time, and the condition I must impose is for the number of particles to stay constant, meaning the norm of my (discretized) wavefunction should not change. Is there a way to do that using the current API? Should I implement something? What do you guys think?

[mauro3]

The norm a conserved quantity of the equations, right? Or do you actually need norm==const as an additional constraint to solve the system?

For the former, there are special classes of integrators which conserve certain quantities https://en.wikipedia.org/wiki/Geometric_integrator. But I don't know much about them. I don't think ODE.jl offers any at the moment but I might be wrong. For the latter, I don't know.

Anyway, it would be good to have that example either as an example in an example/ folder or add it to https://github.com/mauro3/IVPTestSuite.jl.

[acroy]

For methods which do not conserve the norm (like RK) you probably want to (re)normalize your wave function after each time-step. Currently, we don't support this, but I can imagine that the event system (see #11) could be used. Of course you can just use small time-steps and points=:specified to renormalize after each call to ode45 (or whatever).

[acroy]

BTW, you might want to have a look at Expmv.jl and/or Expokit.jl, which could be more efficient for such problems (the latter package also has a mutating method).

[mauro3]

I just stumbled across this http://radio.feld.cvut.cz/matlab/techdoc/math_anal/ch_8_od8.html#670396 bottom example. There they replace one of the equations of the ODE with a algebraic relation of the conserved quantity. Maybe this could work for this problem too? (However, I have no personal experience with using that trick). To solve the DAE you could either try DASSL.jl or #72.

[berceanu]

Well, so what I need to do is basically divide by sum(abs2(psi[i])) after every time step. Perhaps an iterator version of ode45 (#27) would make this easier.

[berceanu]

What I ended up doing for the time being is

h = 0.1
t₀ = 0.
q = 50

aout = Array[ones(Float64, N^2), ones(Float64, N^2)]
for k in 0:q
    tk = t₀ + k*h
    tout, aout = ode45((t,a)->F(t,a, N, m₀,n₀, 25.), aout[end]/norm(aout[end]), [tk, tk+h], points=:specified)
end

[ChrisRackauckas]

Closing this because @mauro3 is right that the right way to do this is via a DAE solver or a geometric integrator. That's a whole separate issue, and packages like Sundials.jl and DASSL.jl exist for this.