Using Chebychev solvers

This web page can be downloaded as notebook: chebychev_solvers.ipynb (Jupyter) or chebychev_solvers.md (Markdown)

Note

This tutorial focuses on the usage of the wavepacket.solver.ChebychevSolver. If you want to know more about the theory, see Background theory: Polynomial solvers or the original paper by Tal-Ezer and Kosloff [1]

Before deciding on the solver to use, be aware of the different tradeoffs between ODE solvers and the Chebychev solver:

Advantages of the Chebychev solver:

much faster, something like a factor of 4-5.
highly accurate with errors on the order of the numerical precision.

Drawbacks of the Chebychev solver:

does not work for time-dependent systems
does not work with complex eigenvalues (open systems, absorbing boundary conditions)
setup is more complex

Setting up a Chebychev solver

To set up the solver, you need three data items:

The expression that describes the equation of motion.
The spectrum of the Hamiltonian / Liouvillian, or rather upper and lower bounds for the maximum and minimum eigenvalue, respectively.
The length of one elementary time step.

The expression is trivial, you always need to define the differential equation to solve.

The tricky part is the determination of the spectral bounds. If these are too generous, we lose efficiency proportional to the quality of the guess; for example, if the Hamiltonian’s spectrum fits twice into our guessed interval, computation takes twice as long as for a perfect guess. If, however, we choose them too tight, that is if the Liouvillian or Hamiltonian has an eigenvalue outside the supplied bounds, the solution usually diverges. We will discuss this in the next section.

The length of the time step is finally given through the alpha value:

$\alpha = \frac{\Delta t \ \delta E}{2},$

where $\delta E$ is the difference between our guessed upper and lower bound of the spectrum (the “spectral range” of the Hamiltonian/Liouvillian). The alpha value should be larger than 40, otherwise the efficiency goes down. If it is much larger than say 100, you might run into numerical problems.

Note

The Chebychev solver is only efficient for rather large time steps. If you need the wave function at small timesteps, it is possible to implement (exact) interpolation, but that has not been done yet. Open a ticket or send a mail to the list if this is an issue for your application.

Determining the spectrum

Let us come to the tricky part: How do we obtain good bounds for the spectrum of the Hamiltonian? As an example, let us set up a simple harmonic oscillator.

import wavepacket as wp

grid = wp.grid.Grid(wp.grid.PlaneWaveDof(-10, 10, 128))

kinetic = wp.operator.CartesianKineticEnergy(grid, 0, mass=1)
potential = wp.operator.Potential1D(grid, 0, lambda x: 0.5 * x ** 2)
hamiltonian = kinetic + potential

equation = wp.expression.SchroedingerEquation(hamiltonian)
liouvillian = wp.expression.CommutatorLiouvillian(hamiltonian)

For didactic reasons, we will discuss three items in the following:

How is the spectrum of the Liouvillian related to that of the Hamiltonian?
How can I estimate the bounds of a Hamiltonian’s spectrum?
Are there better ways to increase the efficiency of the solver?

How can I estimate the bounds of a Hamiltonian’s spectrum?

First, we suggest not to spend too much effort on the exact spectrum, but take reasonable bounds instead. Your time is usually more valuable than the computing time penalty from guessing generous bounds to the spectrum.

A simple lower bound of the Hamiltonian is given by the smallest value of the potential, that is, zero in our harmonic oscillator example.

For an estimate of the largest eigenvalue, we can use the power iteration: Repeatedly applying the Hamiltonian to an almost arbitrary state converges exponentially towards the eigenvector with the largest (absolute) eigenvalue. The implementation is straight forward:

import math

psi = wp.builder.unit_wave_function(grid)
for iteration in range(10):
    psi = equation.apply(psi, t=0)
    psi = wp.normalize(psi)

energy_guess = wp.expectation_value(hamiltonian, psi).real
energy_guess = 1.2 * energy_guess
print(f"Energy guess = {energy_guess:.4}")

Energy guess = 280.1

The initial state can be almost arbitrary; it only needs to contain the highest-energy eigenstate as a non-negligible component. The factor of 1.2 is a guessed safety margin because the result may not have been converged yet. Thus, we arrive at an estimate of the spectrum of about [0, 280]. For an alpha value of 40 or more, our time step must be at least $2 \cdot 40 / 280 = 2/7$ . Let us evolve a Gaussian wave packet with such values:

psi0 = wp.builder.product_wave_function(grid, wp.special.Gaussian(-5, 0, rms=1))
solver = wp.solver.ChebychevSolver(equation, math.pi/10, (0, energy_guess))

for t, psi in solver.propagate(psi0, t0=0.0, num_steps=10):
    wp.log(t, psi)

-----------------------------------------------
t = 0.0,     trace = 1.0

<x_0> = -5.0  +/- 0.707107

-----------------------------------------------
t = 0.314159,     trace = 1.0

<x_0> = -4.75528  +/- 0.707107

-----------------------------------------------
t = 0.628319,     trace = 1.0

<x_0> = -4.04508  +/- 0.707107

-----------------------------------------------
t = 0.942478,     trace = 1.0

<x_0> = -2.93893  +/- 0.707107

-----------------------------------------------
t = 1.25664,     trace = 1.0

<x_0> = -1.54508  +/- 0.707107

-----------------------------------------------
t = 1.5708,     trace = 1.0

<x_0> = -2.25583e-12  +/- 0.707107

-----------------------------------------------
t = 1.88496,     trace = 1.0

<x_0> = 1.54508  +/- 0.707107

-----------------------------------------------
t = 2.19911,     trace = 1.0

<x_0> = 2.93893  +/- 0.707107

-----------------------------------------------
t = 2.51327,     trace = 1.0

<x_0> = 4.04508  +/- 0.707107

-----------------------------------------------
t = 2.82743,     trace = 1.0

<x_0> = 4.75528  +/- 0.707107

-----------------------------------------------
t = 3.14159,     trace = 1.0

<x_0> = 5.0  +/- 0.707107

We can also check what happens if we get the spectrum wrong. Let us say we cut at 200 a.u.:

bad_solver = wp.solver.ChebychevSolver(equation, math.pi/10, (0, 200))

for t, psi in bad_solver.propagate(psi0, t0=0.0, num_steps=10):
    wp.log(t, psi)

-----------------------------------------------
t = 0.0,     trace = 1.0

<x_0> = -5.0  +/- 0.707107

-----------------------------------------------
t = 0.314159,     trace = 86877.9

<x_0> = -0.702138  +/- 9.1446

-----------------------------------------------
t = 0.628319,     trace = 8.93325e+23

<x_0> = -0.702112  +/- 9.14464

-----------------------------------------------
t = 0.942478,     trace = 9.18575e+42

<x_0> = -0.702112  +/- 9.14464

-----------------------------------------------
t = 1.25664,     trace = 9.44538e+61

<x_0> = -0.702112  +/- 9.14464

-----------------------------------------------
t = 1.5708,     trace = 9.71235e+80

<x_0> = -0.702112  +/- 9.14464

-----------------------------------------------
t = 1.88496,     trace = 9.98687e+99

<x_0> = -0.702112  +/- 9.14464

-----------------------------------------------
t = 2.19911,     trace = 1.02691e+119

<x_0> = -0.702112  +/- 9.14464

-----------------------------------------------
t = 2.51327,     trace = 1.05594e+138

<x_0> = -0.702112  +/- 9.14464

-----------------------------------------------
t = 2.82743,     trace = 1.08579e+157

<x_0> = -0.702112  +/- 9.14464

-----------------------------------------------
t = 3.14159,     trace = 1.11648e+176

<x_0> = -0.702112  +/- 9.14464

The footprint of our wrong spectrum is most apparent in the trace. So that is something you should always do in practice. The expectation values become fixed, because the wave function is dominated by the constantly diverging components.

Are there better ways to increase the efficiency of the solver?

So far, we have taken the Hamiltonian for granted and tried to estimate bounds of the spectrum. However, we may just as well define the spectral bounds by truncating the operator. To motivate this approach, let us print the average and standard deviation of the initial state’s energy:

psi0 = wp.builder.product_wave_function(grid, wp.special.Gaussian(-5, 0, rms=1))
energy = wp.expectation_value(hamiltonian, psi0).real
energy2 = wp.expectation_value(hamiltonian*hamiltonian, psi0).real

print(f"E = {energy},   dE^2 = {energy2 - energy**2}")

E = 13.000000000006736,   dE^2 = 12.500000001972978

From this data, it seems safe to assume that the state is well represented using only energy eigenstates with energies <= 30 a.u. This is much less than the estimate of 280 a.u. that we have guessed in the preceding section. In other words, 90% of our computing time is wasted on the correct propagation of high-energy eigenstates. These are irrelevant for our dynamics, but would blow up the calculation if we chose the bounds too small. Such an imbalance between the spectrum that the Hamiltonian can support and the spectrum that we actually need is unfortunately common.

To get rid of this waste, we truncate the Hamiltonian. Because this is difficult, we actually truncate the kinetic and potential energy operators separately. The question is: At what values? Here we need some intuition again (or guesswork or plain trial-and-error).

Our Gaussian (x0 = -5, rms = 1) extends maybe three rms up to x=8, where the potential value is $8^2/2 = 32 \ \mathrm{a.u.}$ . Hence, we might truncate at 35 a.u., which is just a little less than the maximum potential of 50 a.u. For the kinetic energy, we do the same, which has a much larger effect. The system setup changes to

truncated_kinetic = wp.operator.CartesianKineticEnergy(grid, 0, mass=1, cutoff=35)
truncated_potential = wp.operator.Potential1D(grid, 0, lambda x: 0.5 * x ** 2, cutoff=35)
truncated_hamiltonian = truncated_kinetic + truncated_potential
truncated_equation = wp.expression.SchroedingerEquation(truncated_hamiltonian)

We do not know nor care about the exact spectrum of this Hamiltonian, but we know that it cannot be larger than 35 + 35 = 70, which is one fourth of the original Hamiltonian. This allows us to increase the time step by a factor of 4 with the same computational effort.

solver = wp.solver.ChebychevSolver(truncated_equation, 4 * math.pi/10, (0, 70))

for t, psi in solver.propagate(psi0, t0=0.0, num_steps=10):
    wp.log(t, psi)

-----------------------------------------------
t = 0.0,     trace = 1.0

<x_0> = -5.0  +/- 0.707107

-----------------------------------------------
t = 1.25664,     trace = 1.0

<x_0> = -1.54509  +/- 0.707115

-----------------------------------------------
t = 2.51327,     trace = 1.0

<x_0> = 4.04508  +/- 0.707092

-----------------------------------------------
t = 3.76991,     trace = 1.0

<x_0> = 4.04509  +/- 0.707136

-----------------------------------------------
t = 5.02655,     trace = 1.0

<x_0> = -1.54507  +/- 0.7071

-----------------------------------------------
t = 6.28319,     trace = 1.0

<x_0> = -5.0  +/- 0.707104

-----------------------------------------------
t = 7.53982,     trace = 1.0

<x_0> = -1.54511  +/- 0.707158

-----------------------------------------------
t = 8.79646,     trace = 1.0

<x_0> = 4.04506  +/- 0.707069

-----------------------------------------------
t = 10.0531,     trace = 1.0

<x_0> = 4.0451  +/- 0.707174

-----------------------------------------------
t = 11.3097,     trace = 1.0

<x_0> = -1.54505  +/- 0.707109

-----------------------------------------------
t = 12.5664,     trace = 1.0

<x_0> = -4.99999  +/- 0.707107

Using Chebychev solvers

Setting up a Chebychev solver

Determining the spectrum

How is the spectrum of the Liouvillian related to that of the Hamiltonian?

How can I estimate the bounds of a Hamiltonian’s spectrum?

Are there better ways to increase the efficiency of the solver?