Numerov Method (C++)

Differential Equations

This page presents C++ implementations of Numerov's method for solving second-order differential equations. The Numerov method is particularly well-suited for the radial Schrödinger equation.

Standard Numerov

For differential equations of the form:

The standard Numerov method uses the recurrence relation:

where .

Generalized Numerov

For equations with a source term:

The generalized Numerov method handles both terms:

where .

C++ Implementation

The following code implements both standard and generalized Numerov methods with a complete example:

#include <vector>
#include <cmath>
#include <iostream>

// Helper function: linspace
std::vector<double> linspace(double start, double end, size_t num) {
    std::vector<double> result(num);
    double dx = (end - start) / (num - 1);
    for (size_t i = 0; i < num; ++i) {
        result[i] = start + i * dx;
    }
    return result;
}

// Coulomb potential
double coulombPotential(double r, double Z) {
    if (r < 1e-10) return 0.0;  // Avoid singularity
    return -Z / r;
}

/* =======================
   Standard Numerov
   ======================= */
std::vector<double>
Numerov(const std::vector<double>& F,
        double dx,
        double f0,
        double f1)
{
    const std::size_t Nmax = F.size();

    std::vector<double> Solution(Nmax, 0.0);
    Solution[0] = f0;
    Solution[1] = f1;

    const double h2  = dx * dx;
    const double h12 = h2 / 12.0;

    double w0 = (1.0 - h12 * F[0]) * Solution[0];
    double Fx = F[1];
    double w1 = (1.0 - h12 * Fx) * Solution[1];
    double Phi = Solution[1];

    for (std::size_t i = 2; i < Nmax; ++i)
    {
        double w2 = 2.0 * w1 - w0 + h2 * Phi * Fx;

        w0 = w1;
        w1 = w2;

        Fx = F[i];
        Phi = w2 / (1.0 - h12 * Fx);
        Solution[i] = Phi;
    }

    return Solution;
}

/* =======================
   Generalized Numerov
   ======================= */
std::vector<double>
NumerovGen(const std::vector<double>& F,
           const std::vector<double>& U,
           double dx,
           double f0,
           double f1)
{
    const std::size_t Nmax = F.size();

    std::vector<double> Solution(Nmax, 0.0);
    Solution[0] = f0;
    Solution[1] = f1;

    const double h2  = dx * dx;
    const double h12 = h2 / 12.0;

    double w0 = Solution[0] * (1.0 - h12 * F[0]) - h12 * U[0];
    double w1 = Solution[1] * (1.0 - h12 * F[1]) - h12 * U[1];
    double Phi = Solution[1];

    for (std::size_t i = 2; i < Nmax; ++i)
    {
        const double Fx = F[i];
        const double Ux = U[i];

        double w2 = 2.0 * w1 - w0 + h2 * (Phi * Fx + Ux);

        w0 = w1;
        w1 = w2;

        Phi = (w2 + h12 * Ux) / (1.0 - h12 * Fx);
        Solution[i] = Phi;
    }

    return Solution;
}

// Radial Schrödinger equation RHS
std::vector<double> radialRHS(
    double E,
    double l,
    const std::vector<double>& R,
    const std::vector<double>& Veff
) {
    const std::size_t N = R.size();
    std::vector<double> RHS(N, 0.0);

    for (std::size_t i = 0; i < N; ++i) {
        double r = R[i];
        if (r < 1e-10) r = 1e-10;  // Avoid division by zero
        RHS[i] = 2.0 * (-E + 0.5 * l * (l + 1.0) / (r * r) + Veff[i]);
    }

    return RHS;
}

// Example usage for radial Schrödinger equation
int main() {
    std::vector<double> r = linspace(0.01, 20.0, 10000);  // Start slightly away from origin
    std::vector<double> Veff(r.size());

    // Set up potential (e.g., Coulomb potential for hydrogen)
    for (size_t i = 0; i < Veff.size(); i++) {
        Veff[i] = coulombPotential(r[i], 1.0);  // Z = 1
    }

    // Compute RHS for given energy E and angular momentum l
    double E = -0.5;  // Energy in atomic units
    double l = 0.0;   // Angular momentum quantum number
    std::vector<double> rhs = radialRHS(E, l, r, Veff);

    // Solve using Numerov
    double dx = r[1] - r[0];
    std::vector<double> solution = Numerov(rhs, dx, 0.0, 1e-3);

    // Print some results
    std::cout << "Solution at r=0.01: " << solution[0] << std::endl;
    std::cout << "Solution at r=1.0: " << solution[100] << std::endl;
    std::cout << "Solution at r=10.0: " << solution[5000] << std::endl;

    return 0;
}

Key Features