Finite difference and Runge-Kutta for PDEs

[EmpaDoc]

I made a small program to simulate the time development of a 1D wavepacket obeying the Schrodinger equation, mostly in order to learn a new programming language - so in order to not have to invoke big numerical methods packages, I opted for the simplest solution: The standard three-point finite-difference discretization for the Laplacian (kinetic term), and the 4-order Runge-Kutta method to propagate the wavepacket forward.

I notice that sooner or later, my wavepacket will start to develop oscillations on the smallest wavelengths, followed by a blowup. If I include a step potential to reflect against, this happens once the wavepacket encounters the potential. Otherwise it seems to happen when the wavepacket scatters against itself across the periodic boundary (as it expands).

Changing my three-point finite-difference to a five-point stencil only makes the problem worse.

I know that there are much better methods for PDE solving (e.g. Crank-Nicolson). But that's not really the point here. Now I am curious: What is the cause of this behavior? And-can I tweak the FD+RK so I can get it to work?

 

[AlephZero]

The first thing I would try is reducing your time step (e.g. halve it). Possibly (by happenstance) your time step is close to the limit for stability, so it appears to be working OK until the instability has grown from 10^-6 or (10^-16) to be bigger than the true solution.

Using a solution scheme that has a higher order of accuracy than the order of the differential equation almost always results in conditional stability. That's why a 4th order method like RK is rarely "better" than a well-chosen lower order method like Crank Nicholson - there are very few 4th order differential equations in physics!

 

[EmpaDoc]

Forget it! I just hadn't debugged the program properly...
Sorry. :(
At least I know now that FD+RK4 works fine with the Schrodinger equation.

 

[Sumert]

Can you explain what was the problem in your code ? I am also trying such a way but norm of wave function isn't constant and differs from analytical solution (Im trying to propagate free gaussian wave packet).