#10 Create a propagator

Produce a propagator that, given a Hamiltonian operator that defines the problem and a wave function at some time t0, propagates the wave function to a later time t1.

One part is the design of the interface. The other part is the implementation of a propagator that uses Runge-Kutta and similar methods to solve the differential equation. At this point, it would be useful to make use of the GNU Scientific Library; it provides several special functions that we can also use later.

[Anonymous]

• Propgator interface that takes wavefunction and propagates to a later time

  • setStepsize(double dt) sets step width and reinitializes internal variables
  • CTensor propagate(CTensor psi0, int t0, int t1) propagates wave function psi0 from dtt0 to dtt1

• implementation RungeKuttaPropagator using GSL
  - Operator in Constructor; IllegalArgument if null
  - solve differential equation with gsl_odeiv2_step_rk4 (Runge-Kutta 4th order)
  - Problem: How to test? Wrapper around gsl and use/mock the wrapper? No testing for now; to be discussed next time.

[Anonymous]

Chat summary today:

• We need three classes: a wrapper around the GSL, a wrapper around the differential equation, and a propagator to propagate a wave function.

• Definition of the differential equation:

  • We create an interface ODESystem
  • it overloads the operator() such that it gets t, y as input, and returns dy/dt as output.

• Definition of the wrapper

  • A class ODESolver
  • Constructor takes an ODESystem, a relative error (default 1e-6), an absolute error (default 1e-10)
  • there is a function y1 = solve(t0, t1, y0), which takes an initial vector y0 and propagates it from t0 to t1 using the ODESystem. It returns the vector at time t1.
  • Propagation is done with Runge-Kutta 4th order or something like that; choice of the method comes later.
  • The step size is auto-determined by the wrapper, unless this turns out to be cumbersome; it should start with something like t1-t0, and cache the result for future application.

• To test that the wrapper works correctly, there are two tests:

  1. Use a classical harmonic oscillator d^2x/dt^2 = -x and propagate it with the default settings. This should give a solution close to the analytic solution.
  2. Use an exponential function dy/dt = y to check for correct passing of errors. Three test cases:
     • Use the default errors. The numerical solution should be close to the exponential function exp(t).
     • Set the absolute error to 0, and the relative error to something large (1e-2 or more). The solution should deviate significantly.
     • Set the relative error to 0, but the absolute error to something like 1e-6 or so. The solution should be recovered well.
     • Note that the parameters and such should be tweaked that the whole test is run in <= 2 seconds.

[Piotr Bienkiewicz]

The first version of the propagator is ready. A wrapper around the GSL can be found in src/gsl. Both a propagator class and a wrapper around the differential equation are in src/ode subfolder.

Although some unit tests are already implemented and the propagator is functioning, there are a couple of points to catch up with:

1. The propagator uses arrays instead of tensors to solve a ode system. This is a temporary solution. For my convenience it was easier for me to handle with C-style arrays since the GNU GSL-library also do so. Nevertheless, the next step should be an adoption of the propagator to be able to work with tensor-arguments. That should not be a big deal. We'll see :)

2. The remaining test case have to be implemented. More precisely - the second test case, mentioned in the note from 2013-05-05: "Use an exponential function dy/dt = y to check for correct passing of errors. Three test cases: ..."

[Piotr Bienkiewicz]

During code review some shortcomings shows up regarding to the usage of GSL, particular by the fact that an GSL ode system can be treated only using double values.

The odeint V2-library(http://www.codeproject.com/Articles/268589/odeint-v2-Solving-ordinary-differential-equations) follows a different design approach. Since odeint V2 takes advantage of template metaprogramming there are no such limitations.

In the next step i would try out the odeint V2-library.