<jupyter_start><jupyter_text>Finite differences in infinite domains        This notebook originally appeared as a post on my blog.    Because of my friend, [<NAME>](http://cosmoscalibur.com/), I ended up thinking about using a change of variableswhen solving an eigenvalue problem with finite difference. The problemLet's say that we want to solve a differential equation over an infinite domain. A common caseis to solve the [Time independent Schrödinger equation](https://en.wikipedia.org/wiki/Schr%C3%B6dinger_equationTime-independent_equation) subject to a potential $V(x)$. For example$$-\frac{1}{2}\frac{\mathrm{d}^2}{\mathrm{d}x^2}\psi(x) + V(x) \psi(x) = E\psi(x),\quad \forall x\in (-\infty, \infty),$$where we want to find the pairs of eigenvalues/eigenfunctions $(E_n, \psi_n(x))$.What I normally do when using [finite differences](https://en.wikipedia.org/wiki/Finite_difference_method) is to regularly divide the domain. Where I take a _large enough_ domain, so the solution have decayedclose to zero. What I do in this post is to make a change of variable to render the interval finite firstand then regularly divide the transformed domain in finite intervals.  My usual approachMy usual approach is to approach the second derivative with a[centered difference](https://en.wikipedia.org/wiki/Finite_difference_coefficient)for the point $x_i$ like this$$\frac{\mathrm{d}^2 f(x)}{\mathrm{d}x^2} \approx \frac{f(x + \Delta x) - 2 f(x_i) + f(x - \Delta x)}{\Delta x^2}\, ,$$with $\Delta x$ the separation between points.We can solve this in Python with the following snippet:<jupyter_code>import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import eigs


def regular_FD(pot, npts=101, x_max=10, nvals=6):
    """
    Find eigenvalues/eigenfunctions for Schrodinger
    equation for the given potential `pot` using
    finite differences
    """
    x = np.linspace(-x_max, x_max, npts)
    dx = x[1] - x[0]
    D2 = diags([1, -2, 1], [-1, 0, 1], shape=(npts, npts))/dx**2
    V = diags(pot(x))
    H = -0.5*D2 + V
    vals, vecs = eigs(H, k=nvals, which="SM")
    return x, np.real(vals), vecs
# Jupyter notebook plotting setup & imports
%matplotlib notebook
import matplotlib.pyplot as plt

gray = '#757575'
plt.rcParams["figure.figsize"] = 6, 4
plt.rcParams["mathtext.fontset"] = "cm"
plt.rcParams["text.color"] = gray
fontsize = plt.rcParams["font.size"] = 12
plt.rcParams["xtick.color"] = gray
plt.rcParams["ytick.color"] = gray
plt.rcParams["axes.labelcolor"] = gray
plt.rcParams["axes.edgecolor"] = gray
plt.rcParams["axes.spines.right"] = False
plt.rcParams["axes.spines.top"] = False<jupyter_output><empty_output><jupyter_text>Let's consider the [Quantum harmonic oscillator](https://en.wikipedia.org/wiki/Quantum_harmonic_oscillator),that has as eigenvalues$$E_n = n + \frac{1}{2},\quad \forall n = 0, 1, 2 \cdots$$Using the finite difference method we have values that are really close to the analytic ones.<jupyter_code>x, vals, vecs = regular_FD(lambda x: 0.5*x**2, npts=201)
vals
anal_vals = np.array(range(6)) + 0.5
anal_vals
plt.figure()
plt.plot(anal_vals)
plt.plot(vals, "o")
plt.xlabel(r"$n$", fontsize=16)
plt.ylabel(r"$E_n$", fontsize=16)
plt.legend(["Analytic", "Finite differences"])
plt.tight_layout();<jupyter_output><empty_output><jupyter_text>Let's see the eigenfunctions.<jupyter_code>plt.figure()
plt.plot(x, np.abs(vecs[:, :3])**2)
plt.xlim(-6, 6)
plt.xlabel(r"$x$", fontsize=16)
plt.ylabel(r"$|\psi_n(x)|^2$", fontsize=16)
plt.yticks([])
plt.tight_layout();<jupyter_output><empty_output><jupyter_text>One inconvenient with this method is the redundant sampling towards the extreme of the intervals, while under sample the middle part.  Changing the domain Let's now consider the case where we transform the infinite domain to a finite one with a change of variable$$\xi = \xi(x)$$with $\xi \in (-1, 1)$. Two options for this transformation are:- $\xi = \tanh x$; and- $\xi = \frac{2}{\pi} \arctan x$.Making this change of variable the equation we need to solve the following equation$$-\frac{1}{2}\left(\frac{\mathrm{d}\xi}{\mathrm{d}x}\right)^2\frac{\mathrm{d}^2}{\mathrm{d}\xi^2}\psi(\xi) - \frac{\mathrm{d}^2\xi}{\mathrm{d}x^2}\frac{\mathrm{d}}{\mathrm{d}\xi}\psi(\xi) + V(x) \psi(x) = E\psi(x)$$The following snippet solve the eigenproblem for the mapped domain:<jupyter_code>def mapped_FD(pot, fun, dxdxi, dxdxi2, npts=101, nvals=6, xi_tol=1e-6):
    """
    Find eigenvalues/eigenfunctions for Schrodinger
    equation for the given potential `pot` using
    finite differences over a mapped domain on (-1, 1)
    """
    xi = np.linspace(-1 + xi_tol, 1 - xi_tol, npts)
    x = fun(xi)
    dxi = xi[1] - xi[0]
    D2 = diags([1, -2, 1], [-1, 0, 1], shape=(npts, npts))/dxi**2
    D1 = 0.5*diags([-1, 1], [-1, 1], shape=(npts, npts))/dxi
    V = diags(pot(x))
    fac1 = diags(dxdxi(xi)**2)
    fac2 = diags(dxdxi2(xi))
    H = -0.5*fac1.dot(D2) - 0.5*fac2.dot(D1) + V
    vals, vecs = eigs(H, k=nvals, which="SM")
    return x, np.real(vals), vecs<jupyter_output><empty_output><jupyter_text>First transformation: $\xi = \tanh(x)$Let's consider first the transformation$$\xi = \tanh(x)\, .$$In this case$$\frac{\mathrm{d}\xi}{\mathrm{d}x} = 1 - \tanh^2(x) = 1 - \xi^2\, ,$$and$$\frac{\mathrm{d}^2\xi}{\mathrm{d}x^2} = -2\tanh(x)[1 - \tanh^2(x)] = -2\xi[1 - \xi^2]\, .$$<jupyter_code>pot = lambda x: 0.5*x**2
fun = lambda xi: np.arctanh(xi)
dxdxi = lambda xi: 1 - xi**2
dxdxi2 = lambda xi: -2*xi*(1 - xi**2)

x, vals, vecs = mapped_FD(pot, fun, dxdxi, dxdxi2, npts=201)
vals
plt.figure()
plt.plot(anal_vals)
plt.plot(vals, "o")
plt.legend(["Analytic", "Finite differences"])
plt.xlabel(r"$n$", fontsize=16)
plt.ylabel(r"$E_n$", fontsize=16)
plt.tight_layout();
plt.figure()
plt.plot(x, np.abs(vecs[:, :3])**2)
plt.xlim(-6, 6)
plt.xlabel(r"$x$", fontsize=16)
plt.ylabel(r"$|\psi_n(x)|^2$", fontsize=16)
plt.yticks([])
plt.tight_layout();<jupyter_output><empty_output><jupyter_text>Second transformation: $\xi = \frac{2}{\pi}\mathrm{atan}(x)$Let's consider first the transformation$$\xi = \frac{2}{\pi}\mathrm{atan}(x)\, .$$In this case$$\frac{\mathrm{d}\xi}{\mathrm{d}x} = \frac{2}{\pi(1 + x^2)} = \frac{2 \cos^2\xi}{\pi} \, ,$$and$$\frac{\mathrm{d}^2\xi}{\mathrm{d}x^2} = -\frac{4x}{\pi(1 + x^2)^2} = -\frac{4 \cos^4\xi \tan \xi}{\pi}\, .$$<jupyter_code>pot = lambda x: 0.5*x**2
fun = lambda xi: np.tan(0.5*np.pi*xi)
dxdxi = lambda xi: 2/np.pi * np.cos(0.5*np.pi*xi)**2
dxdxi2 = lambda xi: -4/np.pi * np.cos(0.5*np.pi*xi)**4 * np.tan(0.5*np.pi*xi)
x, vals, vecs = mapped_FD(pot, fun, dxdxi, dxdxi2, npts=201)
vals
plt.figure()
plt.plot(anal_vals)
plt.plot(vals, "o")
plt.legend(["Analytic", "Finite differences"])
plt.xlabel(r"$n$", fontsize=16)
plt.ylabel(r"$E_n$", fontsize=16)
plt.tight_layout();
plt.figure()
plt.plot(x, np.abs(vecs[:, :3])**2)
plt.xlabel(r"$x$", fontsize=16)
plt.ylabel(r"$|\psi|^2$", fontsize=16)
plt.xlim(-6, 6)
plt.xlabel(r"$x$", fontsize=16)
plt.ylabel(r"$|\psi_n(x)|^2$", fontsize=16)
plt.yticks([])
plt.tight_layout();<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18883998296180604, 'fraction_numerical': 0.029390884566235977, 'mean_word_length': 4.879799666110183, 'pattern_counts': {'":': 0, '<': 26, '<?xml version=': 0, '>': 26, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '328614', 'license': 'MIT', 'n_tokens_mistral': 2894, 'n_tokens_neox': 2705, 'n_words': 659, 'path': 'downloads/notebooks/.ipynb_checkpoints/finite_diff_map-checkpoint.ipynb', 'repo_name': 'nicoguaro/nicoguaro.github.io', 'size': '7052'}