Testing different finite difference approximations for time independent Schrodinger equation

Using harmonic oscillator for testing

from numba import njit # compile python
from scipy.sparse import linalg # Linear algebra for sparse matrix
import matplotlib.pyplot as plt # plotting facility
import numpy as np

Backward finite difference

from backward_diff_qhamiltonian import BackwardDiffQHamiltonian

for omega in [4., 8.]:
    # Find energies of a harmonic oscillator V = 0.5*(omega*x)**2
    
    @njit
    def v(x):
        return 0.5 * (omega * x) **2
    
    harmonic_osc = BackwardDiffQHamiltonian(
                        x_grid_dim=512,
                        x_amplitude=5.,
                        v=v,
                    )
    # get energies as eigenvalues
    energies, _ = linalg.eigs(harmonic_osc.hamiltonian, k=20)

    # sort eigenvalues by real part
    energies = energies[np.argsort(energies.real)]

    print("\n\nFirst energies for harmonic oscillator with omega = {}".format(omega))
    print(energies)

First energies for harmonic oscillator with omega = 4.0
[-4796.38899524-131.99292806j -4796.38899524+131.99292806j
 -4795.63964345+195.46651829j -4795.63964345-195.46651829j
 -4795.46041608 -67.19104663j -4795.46041608 +67.19104663j
 -4794.35894865  +0.j         -4793.06977476+258.3402195j
 -4793.06977476-258.3402195j  -4788.70689683+320.84164922j
 -4788.70689683-320.84164922j -4782.58702139+382.9656478j
 -4782.58702139-382.9656478j  -4774.8193555 -444.54656508j
 -4774.8193555 +444.54656508j -4765.74902378+505.45531119j
 -4765.74902378-505.45531119j -4755.76159677+565.94861715j
 -4755.76159677-565.94861715j -4744.78474843+626.49079258j]


First energies for harmonic oscillator with omega = 8.0
[-4630.97593133 +35.39331154j -4630.97593133 -35.39331154j
 -4624.89292092+104.64946928j -4624.89292092-104.64946928j
 -4613.94956354+168.73235826j -4613.94956354-168.73235826j
 -4604.19782274+339.0845038j  -4604.19782274-339.0845038j
 -4603.9065095 -224.19124381j -4603.9065095 +224.19124381j
 -4602.88217655-278.12680713j -4602.88217655+278.12680713j
 -4602.2665706 +404.20602281j -4602.2665706 -404.20602281j
 -4596.4636516 -471.16645791j -4596.4636516 +471.16645791j
 -4586.36191708-538.79816033j -4586.36191708+538.79816033j
 -4571.49322447+605.50942283j -4571.49322447-605.50942283j]

Forward finite difference

from forward_diff_qhamiltonian import ForwardDiffQHamiltonian

for omega in [4., 8.]:
    # Find energies of a harmonic oscillator V = 0.5*(omega*x)**2
    
    @njit
    def v(x):
        return 0.5 * (omega * x) **2
    
    harmonic_osc = ForwardDiffQHamiltonian(
                        x_grid_dim=512,
                        x_amplitude=5.,
                        v=v,
                    )
    # get energies as eigenvalues 
    energies, _ = linalg.eigs(harmonic_osc.hamiltonian, k=20)

    # sort eigenvalues by real part
    energies = energies[np.argsort(energies.real)]

    print("\n\nFirst energies for harmonic oscillator with omega = {}".format(omega))
    print(energies)

First energies for harmonic oscillator with omega = 4.0
[-4785.94534256 +31.37800436j -4785.94534256 -31.37800436j
 -4785.26390819 +94.10807967j -4785.26390819 -94.10807967j
 -4783.56443732+156.71202734j -4783.56443732-156.71202734j
 -4780.53902614+219.11423295j -4780.53902614-219.11423295j
 -4775.94869943-281.34907035j -4775.94869943+281.34907035j
 -4769.52396221+343.41451144j -4769.52396221-343.41451144j
 -4761.0000548 +405.15670292j -4761.0000548 -405.15670292j
 -4750.1431072 +466.21042042j -4750.1431072 -466.21042042j
 -4736.93984194+525.85754297j -4736.93984194-525.85754297j
 -4722.11942452+583.24004566j -4722.11942452-583.24004566j]


First energies for harmonic oscillator with omega = 8.0
[-4616.92456002 +34.71415564j -4616.92456002 -34.71415564j
 -4616.44332502+104.37638145j -4616.44332502-104.37638145j
 -4614.50858608+174.51366206j -4614.50858608-174.51366206j
 -4610.44929953+245.20442304j -4610.44929953-245.20442304j
 -4604.10640925-316.202231j   -4604.10640925+316.202231j
 -4595.68153876+386.6120605j  -4595.68153876-386.6120605j
 -4584.98769126+455.55000735j -4584.98769126-455.55000735j
 -4570.92577202+521.76016143j -4570.92577202-521.76016143j
 -4541.12451   +700.40220835j -4541.12451   -700.40220835j
 -4529.53240905+767.27619173j -4529.53240905-767.27619173j]

Central finite difference

from central_diff_qhamiltonian import CentralDiffQHamiltonian

for omega in [4., 8.]:
    # Find energies of a harmonic oscillator V = 0.5*(omega*x)**2
    
    @njit
    def v(x):
        return 0.5 * (omega * x) **2
    
    harmonic_osc = CentralDiffQHamiltonian(
                        x_grid_dim=512,
                        x_amplitude=5.,
                        v=v,
                    )
    plt.figure(dpi=350)
    # plot eigenfunctions
    for n in range(4):
        plt.plot(harmonic_osc.x, harmonic_osc.get_eigenstate(n), label=str(n))

    print("\n\nFirst energies for harmonic oscillator with omega = {}".format(omega))

    # set precision for printing arrays
    np.set_printoptions(precision=4)
    print(harmonic_osc.energies)


    plt.title("Eigenfunctions for harmonic oscillator with omega = {} (a.u.)".format(omega))
    plt.xlabel('$x$ (a.u.)')
    plt.ylabel('wave functions ($\\psi_n(x)$)')
    plt.legend()
    plt.show()

First energies for harmonic oscillator with omega = 4.0
[ 1.9998  5.999   9.9975 13.9952 17.9922 21.9884 25.9838 29.9784 33.9723
 37.9654 41.9578 45.9494 49.9402 53.9303 57.9196 61.9081 65.8959 69.8829
 73.8691 77.8546]

First energies for harmonic oscillator with omega = 8.0
[  3.9992  11.9962  19.9901  27.9809  35.9687  43.9534  51.9351  59.9137
  67.8892  75.8617  83.831   91.7974  99.7606 107.7208 115.6779 123.6319
 131.5829 139.5307 147.4755 155.4172]