Technical Origin of OpenML - A Quantum Physics Problem

2161 words

11 minutes

Technical Origin of OpenML - A Quantum Physics Problem

2026-01-26

Technology

This post is a comprehensive complement of the Quantum Physics paperMatrix Numerov method for solving Schrödinger’s equation, which present both an intellectual interest and an important historical root of underneath OpenML

Prerequisites
General mastery of Fundamental Calculus and Differential Equations are assumed. We recommend the following texts:

Calculus With Analytic Geometry, George Simmons

Differential Equations and Boundary Value Problems, Edwards, Penney, & Calvis

The Numerov method is a specialized integration formula for numerically integrating differential equations of the form

\psi''(x) = f(x)\psi(x)

It’s an ordinary differential equations of second order in which the first-order term does not appear

The Numerov method was developed by the Russian astronomer and geophysicist Boris Vasilievich Numerov (1891–1941) who originally developed the method in the 1920s to solve problems in celestial mechanics, specifically for calculating the orbits of celestial bodies (integration of perturbations). He published the method in paper “A method of extrapolation of perturbations”, in which Numerov originally applied the algorithm to celestial mechanics for calculating the orbits of astronomical bodies, which involves solving second-order differential equations similar in form to the Schrödinger equation (where the first derivative term is absent).

Boris Numerov was a prominent scientist who served as the director of the Central Astronomical Observatory of Geophysics. Tragically, he was arrested in 1936 during the Stalinist purges and executed in 1941 🕯️. The method remains his most enduring legacy in computational physics 🌺.

TIP
While the paper applies the method to the Schrödinger equation (a problem in quantum mechanics), Numerov designed it to solve the same type of differential equation found in astronomy: second-order equations where the first-derivative term is missing ( $y'' = f(x)y$ )
In his 1924 original paper, B. V. Numerov assumes the reader is already familiar with the standard equations of celestial mechanics (which look like $\frac{d^2x}{dt^2} = -\frac{k^2x}{r^3} + \text{perturbations}$ ) and refers to it in the very first sentence: “numerical integration of the equations of perturbed motion”.

The “equations of perturbed motion” that Numerov refers to in the 1924 paper are mathematically identical in structure to the Schrödinger equation. This shared structure is exactly why the Numerov method works for both astronomy and quantum mechanics.

Both problems involve a second-order linear differential equation where the first-derivative term is missing.

Component	Celestial Mechanics (Numerov, 1924)	Quantum Mechanics (Schrödinger)
Independent Variable	Time ( $t$ )	Position ( $x$ )
Dependent Variable	Planet Coordinate ( $x$ )	Wavefunction ( $\psi$ )
Function $f(\text{var})$	Gravity Factor $\left( -\frac{k^2}{r^3} \right)$	Scaled Potential $\frac{2m}{\hbar^2}[V(x) - E]$
Equation	$\frac{d^2x}{dt^2} = -\frac{k^2}{r^3}x$	$\frac{d^2\psi}{dx^2} = \frac{2m}{\hbar^2}(V-E)\psi$

In the table above, $r$ is the radius vector which changes with time, making the gravity factor a function of $t$ , just as the potential $V$ is a function of $x$ .

WARNING
For the context of Schrödinger’s equation only, the Matrix Numerov method works only for time-independent 1D Schrödinger’s equation of the form
$-\frac{\hbar ^{2}}{2m}\frac{d^{2}\psi (x)}{dx^{2}}+V(x)\psi (x)=E\psi (x)$

The standard result known as the Numerov’s method (or Numerov scheme) is exactly what we see presented as equation (2) from the paper:

It is a fourth-order linear multistep method used to solve second-order differential equations of the form $\psi''(x) = f(x)\psi(x)$ where the first-derivative term is absent. The derivation relies on the Taylor expansion of the function $\psi(x)$ and the specific differential equation given in equation

\psi''(x) = f(x)\psi(x)

Deriving Equation (2)#

The paper proceeds with the result of the equation (2) as shown below:

We patch the missing derivation of this equation starting with Taylor Expansion

Step 1: Taylor Expansion#

Taylor Expansion
Let $f: \mathbb{R} \to \mathbb{R}$ be a function that is infinitely differentiable at a real number $a$ . The Taylor series of $f$ centered at $a$ is the power series:
$f(x) = \sum_{n=0}^{\infty} \frac{f^{(n)}(a)}{n!} (x-a)^n$
where:

$n!$ denotes the factorial of $n$ .

$f^{(n)}(a)$ denotes the $n$ -th derivative of $f$ evaluated at the point $a$ .

The zeroth derivative $f^{(0)}$ is defined as the function $f$ itself, and $0! = 1$ .

In the context of numerical analysis, we often use Taylor’s Theorem, which states that for a function $f$ with $k+1$ derivatives on an interval containing $a$ and $x$ , the function can be approximated by a polynomial of degree $k$ (the $k$ -th order Taylor polynomial) plus a remainder term $R_k(x)$ :
$f(x) = f(a) + f'(a)(x-a) + \frac{f''(a)}{2!}(x-a)^2 + \dots + \frac{f^{(k)}(a)}{k!}(x-a)^k + R_k(x)$
The remainder term, often expressed in Peano form as $h_k(x)(x-a)^k$ where $\lim_{x\to a} h_k(x) = 0$ , represents the error of the approximation and is commonly denoted using Big O notation as $O(|x-a|^{k+1})$ .

To derive the equation used in the Numerov method paper, we apply the Taylor series definition to the wave function $\psi(x)$ . Let $f(x) = \psi(x)$

CAUTION
We assume $\psi(x)$ is smooth and differentiable.

Then we expand around the current grid point, so let $a = x_i$ . We want to find the value at the next grid point, so let $x = x_i + d$ . The difference term $(x - a)$ , therefore, becomes:

(x - a) = (x_i + d) - x_i = d

Substituting these values into the formal definition up to the 4th order ( $k=4$ ) gives us:

\psi(x_i + d) = \frac{\psi^{(0)}(x_i)}{0!} d^0 + \frac{\psi^{(1)}(x_i)}{1!} d^1 + \frac{\psi^{(2)}(x_i)}{2!} d^2 + \frac{\psi^{(3)}(x_i)}{3!} d^3 + \frac{\psi^{(4)}(x_i)}{4!} d^4 + R_4(x)

Now we expand $\psi(x)$ around a point $x_i$ to find the values at the neighboring points $x_{i+1} = x_i + d$ and $x_{i-1} = x_i - d$ .

\psi(x_i + d) = \psi(x_i) + d\psi'(x_i) + \frac{d^2}{2!}\psi''(x_i) + \frac{d^3}{3!}\psi'''(x_i) + \frac{d^4}{4!}\psi^{(4)}(x_i) + O(d^5)

\psi(x_i - d) = \psi(x_i) - d\psi'(x_i) + \frac{d^2}{2!}\psi''(x_i) - \frac{d^3}{3!}\psi'''(x_i) + \frac{d^4}{4!}\psi^{(4)}(x_i) + O(d^5)

When we sum these two expansions, the odd-power terms (first and third derivatives) cancel out:

\psi_{i+1} + \psi_{i-1} = 2\psi_i + d^2 \psi''_i + \frac{d^4}{12} \psi^{(4)}_i + O(d^6)

Step 2: Substitution of the Differential Equation#

We use the differential equation provided in the paper, $\psi''(x) = f(x)\psi(x)$ . This allows us to replace the second derivative $\psi''_i$ with $f_i \psi_i$ :

\psi_{i+1} + \psi_{i-1} = 2\psi_i + d^2 f_i \psi_i + \frac{d^4}{12} \psi^{(4)}_i

We still have the fourth derivative $\psi^{(4)}_i$ to deal with. We can find an expression for it by differentiating the original equation $\psi'' = f\psi$ twice:

\psi^{(4)}(x) = \frac{d^2}{dx^2} \left( f(x)\psi(x) \right)

Step 3: Discretizing the Fourth Derivative#

Instead of calculating $(f\psi)''$ analytically (which would be complex), we approximate it using the same central difference formula we essentially derived in Step 1. The second derivative of any function $y(x)$ can be approximated as $\frac{y_{i+1} - 2y_i + y_{i-1}}{d^2}$ with an error of $O(d^2)$ .

Applying this to $y = f\psi$ :

\psi^{(4)}_i = (f\psi)''_i \approx \frac{f_{i+1}\psi_{i+1} - 2f_i\psi_i + f_{i-1}\psi_{i-1}}{d^2}

Substituting this approximation back into our sum equation:

\psi_{i+1} + \psi_{i-1} = 2\psi_i + d^2 f_i \psi_i + \frac{d^4}{12} \left[ \frac{f_{i+1}\psi_{i+1} - 2f_i\psi_i + f_{i-1}\psi_{i-1}}{d^2} \right]

Simplifying the $d^4$ term:

\psi_{i+1} + \psi_{i-1} = 2\psi_i + d^2 f_i \psi_i + \frac{d^2}{12} (f_{i+1}\psi_{i+1} - 2f_i\psi_i + f_{i-1}\psi_{i-1})

Step 4: Rearrangement to Equation (2)#

Now we group the terms by their indices ( $i+1$ , $i$ , $i-1$ ).

Group $\psi_{i+1}$ terms:

\psi_{i+1} - \frac{d^2}{12} f_{i+1}\psi_{i+1} = \psi_{i+1} \left( 1 - \frac{d^2}{12} f_{i+1} \right) = \psi_{i+1} \left( \frac{12 - d^2 f_{i+1}}{12} \right)

Group $\psi_{i-1}$ terms:

\psi_{i-1} - \frac{d^2}{12} f_{i-1}\psi_{i-1} = \psi_{i-1} \left( \frac{12 - d^2 f_{i-1}}{12} \right)

Group $\psi_{i}$ terms:

2\psi_i + d^2 f_i \psi_i - \frac{2d^2}{12} f_i \psi_i = \psi_i \left( 2 + \frac{10d^2}{12} f_i \right) = \psi_i \left( \frac{24 + 10d^2 f_i}{12} \right)

The equation becomes:

\psi_{i+1} \left( \frac{12 - d^2 f_{i+1}}{12} \right) + \psi_{i-1} \left( \frac{12 - d^2 f_{i-1}}{12} \right) = \psi_i \left( \frac{24 + 10d^2 f_i}{12} \right)

Multiplying the entire equation by 12 to remove denominators:

\psi_{i+1} (12 - d^2 f_{i+1}) + \psi_{i-1} (12 - d^2 f_{i-1}) = 2\psi_i (12 + 5d^2 f_i)

Solving for $\psi_{i+1}$ :

\psi_{i+1} (12 - d^2 f_{i+1}) = 2\psi_i (12 + 5d^2 f_i) - \psi_{i-1} (12 - d^2 f_{i-1})

Finally, divide by $(12 - d^2 f_{i+1})$ . To exactly match the paper’s notation in equation (2), we multiply the numerator and denominator of the fraction by $-1$ (flipping the signs):

\psi_{i+1} = \frac{\psi_{i-1}(12 - d^2 f_{i-1}) - 2\psi_i(12 + 5d^2 f_i)}{d^2 f_{i+1} - 12}

This matches equation (2) in the paper.

Matrix Numerov Simulator#

Here is a Python implementation of the Matrix Numerov method to solve the dimensionless Schrödinger equation for the Simple Harmonic Oscillator (SHO) with the potential of $V(x) = \frac{1}{2}x^2$ (Harmonic Oscillator in natural units $\hbar=m=\omega=1$ ). The expected energy $E_n = n + \frac{1}{2}$ (i.e., 0.5, 1.5, 2.5, …)

Explanation of the Implementation:

Matrix Construction:
- Matrix A represents the second derivative operator $\frac{d^2}{dx^2}$ . In the finite difference method, this leads to the familiar $[1, -2, 1]$ pattern on the diagonals.
- Matrix B is the unique feature of the Numerov method. It modifies the kinetic energy term to account for fourth-order error terms ( $O(d^6)$ locally), making it significantly more accurate than standard finite difference methods which are usually only second-order accurate.
- Matrix V is simply the potential energy evaluated at each grid point on the diagonal.
The Hamiltonian: The code solves equation (5) from the paper: $-\frac{\hbar^2}{2m} B^{-1} A \psi + V \psi = E \psi$ .
- We inverted $B$ explicitly using np.linalg.inv. For 1D problems where $N$ is small (hundreds or thousands), this is computationally cheap and easy to read.
Boundary Conditions: By using a fixed grid size $N$ from $-L$ to $L$ and not connecting the ends, we implicitly enforce “hard wall” boundary conditions ( $\psi=0$ outside the grid). This is why choosing a sufficiently large $L$ (like $L=6.0$ for low energy states) is important—the wavefunction must naturally decay to zero before hitting the edge.

1
import numpy as np
2
import matplotlib.pyplot as plt
3

4
def solve_sho_matrix_numerov(N=200, L=6.0):
5
    """
6
    Solves the 1D Simple Harmonic Oscillator using the Matrix Numerov method.
7

8
    Args:
9
        N (int): Number of grid points.
10
        L (float): Spatial boundary +/- L.
11

12
    Returns:
13
        x (ndarray): The spatial grid.
14
        eigenvalues (ndarray): Calculated energy levels.
15
        eigenvectors (ndarray): Normalized wave functions (columns).
16
    """
17
    # 1. Define the Grid
18
    # We use a linspace from -L to L
19
    x = np.linspace(-L, L, N)
20
    d = x[1] - x[0]  # Grid spacing
21

22
    # 2. Construct the Matrix A (Second Derivative Approximation)
23
    # A = (I_{i-1} - 2I_i + I_{i+1}) / d^2
24
    # Diagonals: 1, -2, 1
25
    main_diag_A = -2.0 * np.ones(N)
26
    off_diag_A = np.ones(N - 1)
27
    A = (np.diag(main_diag_A) + np.diag(off_diag_A, k=1) + np.diag(off_diag_A, k=-1)) / d**2
28

29
    # 3. Construct the Matrix B (Numerov Consistency Matrix)
30
    # B = (I_{i-1} + 10I_i + I_{i+1}) / 12
31
    # Diagonals: 1, 10, 1
32
    main_diag_B = 10.0 * np.ones(N)
33
    off_diag_B = np.ones(N - 1)
34
    B = (np.diag(main_diag_B) + np.diag(off_diag_B, k=1) + np.diag(off_diag_B, k=-1)) / 12.0
35

36
    # 4. Construct the Potential Matrix V
37
    # V is diagonal with V_ii = 0.5 * x_i^2
38
    V_diagonal = 0.5 * x**2
39
    V = np.diag(V_diagonal)
40

41
    # 5. Construct the Hamiltonian Matrix
42
    # Form: H = - (1/2) * B^(-1) * A + V
43
    # Note: We assume h_bar = 1, mass = 1.
44
    # We use np.linalg.inv for B. Since B is strictly diagonally dominant, it is well-conditioned.
45
    B_inv = np.linalg.inv(B)
46
    Kinetic = -0.5 * np.dot(B_inv, A)
47
    H = Kinetic + V
48

49
    # 6. Solve the Eigenvalue Problem
50
    # Since H is not strictly symmetric in this representation (B^-1 * A is not necessarily symmetric),
51
    # we use the general eigensolver `eig` rather than `eigh`.
52
    # However, physically, the eigenvalues must be real.
53
    eigenvalues, eigenvectors = np.linalg.eig(H)
54

55
    # 7. Sort and Clean Results
56
    # Sort by real part of eigenvalues
57
    idx = eigenvalues.argsort()
58
    eigenvalues = eigenvalues[idx]
59
    eigenvectors = eigenvectors[:, idx]
60

61
    # Normalize wavefunctions (Simple sum of squares approx)
62
    # integral |psi|^2 dx = 1  =>  sum |psi|^2 * d = 1
63
    for i in range(N):
64
        norm = np.sqrt(np.sum(eigenvectors[:, i]**2) * d)
65
        eigenvectors[:, i] /= norm
66

67
    # Enforce consistent phase (make max amplitude positive for plotting consistency)
68
    for i in range(N):
69
        if np.max(eigenvectors[:, i]) < abs(np.min(eigenvectors[:, i])):
70
            eigenvectors[:, i] *= -1
71

72
    return x, eigenvalues.real, eigenvectors.real
73

74
def plot_results(x, eigenvalues, eigenvectors, num_states=4):
75
    """Plots the potential and the first few wave functions offset by their energies."""
76
    plt.figure(figsize=(10, 7))
77

78
    # Plot Potential V(x)
79
    plt.plot(x, 0.5 * x**2, 'k--', linewidth=1.5, label=r'$V(x) = \frac{1}{2}x^2$')
80

81
    # Plot Wavefunctions
82
    # We scale the amplitude for visibility: Energy + scale * psi
83
    scale_factor = 0.8
84
    colors = ['#d62728', '#1f77b4', '#ff7f0e', '#2ca02c'] # Red, Blue, Orange, Green
85

86
    for n in range(num_states):
87
        E = eigenvalues[n]
88
        psi = eigenvectors[:, n]
89

90
        # Determine label with calculated energy
91
        label = f'$n={n}, E_{{calc}}={E:.4f}$'
92

93
        # Plot the wavefunction baseline at E
94
        plt.axhline(E, color='gray', linestyle=':', alpha=0.5)
95
        plt.plot(x, E + scale_factor * psi, color=colors[n % len(colors)], label=label, linewidth=2)
96

97
    plt.title('Schrödinger Solutions: Harmonic Oscillator via Matrix Numerov')
98
    plt.xlabel('Position $x$')
99
    plt.ylabel('Energy (and $\psi(x)$ amplitude)')
100
    plt.ylim(0, eigenvalues[num_states-1] + 1.5)
101
    plt.legend(loc='lower right')
102
    plt.grid(True, alpha=0.3)
103
    plt.tight_layout()
104
    plt.show()
105

106
if __name__ == "__main__":
107
    # Parameters
108
    # N=300 provides high accuracy for the first few states
109
    x_grid, energies, wavefuncs = solve_sho_matrix_numerov(N=300, L=6.0)
110

111
    # Print accuracy check
112
    print(f"{'State':<10} | {'Calc. Energy':<15} | {'Exact Energy':<15} | {'Error':<15}")
113
    print("-" * 65)
114
    for n in range(5):
115
        exact = n + 0.5
116
        error = abs(energies[n] - exact)
117
        print(f"n={n:<8} | {energies[n]:<15.6f} | {exact:<15.6f} | {error:.2e}")
118

119
    plot_results(x_grid, energies, wavefuncs)

Simply running the code above with python3 solver.py would give us the following eigenstates printed in terminal:

1
tate      | Calc. Energy    | Exact Energy    | Error
2
-----------------------------------------------------------------
3
n=0        | 0.500000        | 0.500000        | 1.01e-08
4
n=1        | 1.500000        | 1.500000        | 7.10e-08
5
n=2        | 2.500000        | 2.500000        | 2.53e-07
6
n=3        | 3.499999        | 3.500000        | 6.39e-07
7
n=4        | 4.499999        | 4.500000        | 1.31e-06

and the following wave functions:

One can visit matrix-numerov.openml.io for an interactive version:

(To be continued…)

Technical Origin of OpenML - A Quantum Physics Problem

https://blogs.openml.io/posts/matrix-numerov/

Author

OpenML Blogs

Published at

2026-01-26

License

CC BY-NC-SA 4.0

Transformer

1

Deriving Equation (2)

Step 1: Taylor Expansion

Step 2: Substitution of the Differential Equation

Step 3: Discretizing the Fourth Derivative

Step 4: Rearrangement to Equation (2)

2

Matrix Numerov Simulator