Numerical Simulation of the Gravitational Potential Generated by a Stellar Mass with the 2D Poisson Equation

Papadopol Lucian Ioan
l.i.papadopol@gmail.com

1. Objective of the simulation

In this notebook we numerically simulate the gravitational potential generated by a mass concentrated in space, idealized as a stellar mass.

The goal is not to simulate the motion of a star or a planet over time, but to compute the potential field produced by a given mass distribution.

The guiding question is: given a massive body in space, what gravitational potential does it generate around itself?

2. Connection with elliptic PDEs

Elliptic partial differential equations often describe equilibrium problems or stationary fields.

the gravitational potential;
the electric potential;
the stationary temperature in a plate;
stationary pressure distributions;
elastic equilibrium problems.

In these problems we are not necessarily studying how the system evolves in time. Instead, we look for a spatial function that satisfies an equation throughout the domain.

In this case the unknown function is:

\phi(x,y)

where $\phi$ represents the gravitational potential.

3. Difference from evolutionary PDEs

In the previous notebooks we studied PDEs involving functions such as:

u(x,t)

that is, functions depending on space and time. For example, in 1D convection and diffusion we had:

u = u(x,t)

and therefore we built a space-time grid.

In a stationary elliptic PDE, instead, the unknown function depends only on space. In the 2D case:

\phi = \phi(x,y)

So we do not build a space-time grid, but a two-dimensional spatial grid. The problem is not to advance in time, but to find a configuration satisfying the equation at every point of the grid.

4. Simulated physical system

We consider a square region of space. At the center of the domain we place a concentrated mass distribution, which represents a simplified stellar mass.

The mass generates a gravitational potential in the surrounding space. In the Newtonian model, the gravitational field can be obtained from the potential as:

\vec{g} = -\nabla \phi

where:

$\phi(x,y)$ is the gravitational potential;
$\vec{g}$ is the gravitational field;
$\nabla \phi$ is the gradient of the potential.

The gravitational field points toward the mass, that is, toward the region where the potential has a minimum.

5. Poisson equation

The mathematical model used in this notebook is the Poisson equation:

\nabla^2 \phi = \rho

where:

$\phi(x,y)$ is the gravitational potential;
$\rho(x,y)$ represents the mass distribution;
$\nabla^2$ is the Laplacian operator.

For the complete Newtonian gravitational potential one would write:

\nabla^2 \Phi = 4\pi G \rho

where $G$ is the universal gravitational constant.

In this notebook we use the normalized form:

\nabla^2 \phi = \rho

This avoids very large or very small physical constants and lets us focus on the numerical method.

6. Meaning of the Laplacian

In two dimensions the Laplacian is defined as:

\nabla^2 \phi = \frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2}

Therefore the Poisson equation becomes:

\frac{\partial^2 \phi}{\partial x^2} + \frac{\partial^2 \phi}{\partial y^2} = \rho

Intuitively, the Laplacian measures how much the value of a function at a point differs from the values around it. In an elliptic PDE, the potential at each point is linked to the values at neighboring points and to the source term located there.

In this notebook, the source is the mass distribution.

7. Spatial domain

We consider a square domain:

0 \le x \le L

0 \le y \le L

The spatial grid is made of points:

x_i = i\Delta x

y_j = j\Delta y

where:

i = 0,1,2,\dots,N_x-1

j = 0,1,2,\dots,N_y-1

The numerical value:

\phi_{i,j}

approximates:

\phi(x_i,y_j)

that is, the gravitational potential at the grid point with coordinates $(x_i,y_j)$ .

8. Boundary conditions

To solve an elliptic PDE, it is necessary to specify what happens on the boundary of the domain.

Here we use Dirichlet boundary conditions. We impose:

\phi = 0

on the whole boundary of the domain.

This means that the potential is fixed to zero along the outer contour of the simulated region. In the code this corresponds to:

\phi_{0,j}=0,\qquad \phi_{N_x-1,j}=0,\qquad \phi_{i,0}=0,\qquad \phi_{i,N_y-1}=0

for all indices $i$ and $j$ .

This choice is numerically simple and represents the idea that the boundary of the domain is far enough from the mass that we can set a reference value there.

9. Mass distribution

An ideal point mass would produce a mathematical singularity. To avoid numerical problems, we represent the stellar mass with a smooth distribution concentrated at the center of the domain.

A simple choice is a Gaussian function:

\rho(x,y) = A \exp\left(-\frac{(x-x_c)^2+(y-y_c)^2}{2\sigma^2}\right)

where:

$A$ controls the intensity of the mass;
$(x_c,y_c)$ is the position of the center of the mass;
$\sigma$ controls how concentrated the mass is.

If $\sigma$ is small, the mass is highly concentrated. If it is larger, the mass is spread over a wider region.

10. Discretization of the Laplacian

To solve the Poisson equation numerically, we use finite differences. Assuming $\Delta x = \Delta y$ for simplicity, the second derivative with respect to $x$ is approximated by:

\frac{\partial^2 \phi}{\partial x^2}\approx\frac{\phi_{i+1,j}-2\phi_{i,j}+\phi_{i-1,j}}{\Delta x^2}

The second derivative with respect to $y$ is approximated by:

\frac{\partial^2 \phi}{\partial y^2}\approx\frac{\phi_{i,j+1}-2\phi_{i,j}+\phi_{i,j-1}}{\Delta x^2}

Adding the two contributions gives:

\nabla^2 \phi \approx \frac{\phi_{i+1,j}+\phi_{i-1,j}+\phi_{i,j+1}+\phi_{i,j-1}-4\phi_{i,j}}{\Delta x^2}

11. Numerical scheme for Poisson equation

Starting from:

\nabla^2 \phi = \rho

and substituting the discrete Laplacian, we get:

\frac{\phi_{i+1,j}+\phi_{i-1,j}+\phi_{i,j+1}+\phi_{i,j-1}-4\phi_{i,j}}{\Delta x^2}=\rho_{i,j}

Solving for the central value gives:

\phi_{i,j}=\frac{1}{4}\left(\phi_{i+1,j}+\phi_{i-1,j}+\phi_{i,j+1}+\phi_{i,j-1}-\rho_{i,j}\Delta x^2\right)

This formula says that the potential at one point depends on the four neighboring points and on the source term at that point.

12. Jacobi method

To solve the numerical system we use the Jacobi iterative method. The idea is:

start from an initial guess for the potential;
compute a new grid using the values from the previous grid;
repeat many times;
stop when the solution changes very little.

The iterative formula is:

\phi_{i,j}^{(k+1)}=\frac{1}{4}\left(\phi_{i+1,j}^{(k)}+\phi_{i-1,j}^{(k)}+\phi_{i,j+1}^{(k)}+\phi_{i,j-1}^{(k)}-\rho_{i,j}\Delta x^2\right)

The index $k$ does not represent physical time. It only represents the number of iterations of the numerical method.

13. Physical time versus numerical iteration

In the convection and diffusion problems we had values such as:

u_i^n

where $n$ represented a true physical time level.

In this problem, instead, we have:

\phi_{i,j}^{(k)}

where $k$ is only the Jacobi iteration index.

In evolutionary PDEs, the time step has physical meaning.
In stationary elliptic PDEs, iteration is only a numerical procedure used to reach the solution.

The Jacobi method does not simulate the real evolution of the system. It is only a way to find the stationary field satisfying the Poisson equation.

14. Stopping criterion

During the iterations we can monitor how much the solution changes. For example, the iterative error can be defined as:

\text{error}=\max\left|\phi^{(k+1)}-\phi^{(k)}\right|

If this value becomes very small, the solution is no longer changing significantly. We can stop the iterative cycle when:

\text{error}<\text{tolerance}

15. Gravitational field

Once the potential $\phi(x,y)$ has been computed, we can obtain the gravitational field:

\vec{g}=-\nabla \phi

In two dimensions:

g_x=-\frac{\partial \phi}{\partial x}

g_y=-\frac{\partial \phi}{\partial y}

Numerically, the gradient can be approximated with finite differences. The field will be represented with arrows pointing toward the central mass.

16. Implementation

Jacobi iterator

import numpy as np

def poisson_jacobi(phi0, rhs, dx, dy, max_iter=10000, tol=1e-6):
    """
    Solve the 2D Poisson equation with the Jacobi iterative method.

    The equation solved is:

        ∇² phi = rhs

    Parameters
    ----------
    phi0 : array
        Initial numerical guess for the potential.
    rhs : array
        Right-hand side of the Poisson equation.
    dx, dy : float
        Grid spacing in the x and y directions.
    max_iter : int
        Maximum number of Jacobi iterations.
    tol : float
        Stopping tolerance based on the maximum change between iterations.
    """

    phi = np.asarray(phi0, dtype=float).copy()

    dx2 = dx * dx
    dy2 = dy * dy

    denom = 2.0 * (dx2 + dy2)

    for it in range(max_iter):
        old = phi.copy()

        phi[1:-1, 1:-1] = (
            (old[1:-1, 2:] + old[1:-1, :-2]) * dy2
            +
            (old[2:, 1:-1] + old[:-2, 1:-1]) * dx2
            -
            rhs[1:-1, 1:-1] * dx2 * dy2
        ) / denom

        err = np.max(np.abs(phi - old))

        if err < tol:
            break

    return phi, it + 1, err

Domain, source and boundary conditions

import numpy as np
import matplotlib.pyplot as plt

# Size of the square domain
L = 1.0

# Number of grid points
nx = 101
ny = 101

# Spatial grid
x = np.linspace(0.0, L, nx)
y = np.linspace(0.0, L, ny)

dx = x[1] - x[0]
dy = y[1] - y[0]

# 2D grid
X, Y = np.meshgrid(x, y)

# Center of the stellar mass
x_c = 0.5
y_c = 0.5

# Amplitude and width of the mass distribution
A = 100.0
sigma = 0.05

# Source term of the Poisson equation
rho = A * np.exp(
    -((X - x_c)**2 + (Y - y_c)**2) / (2 * sigma**2)
)

rhs = rho

# Initial numerical guess for the potential
phi0 = np.zeros((ny, nx))

# Dirichlet boundary conditions:
# phi = 0 on the boundary of the domain
phi0[0, :] = 0.0
phi0[-1, :] = 0.0
phi0[:, 0] = 0.0
phi0[:, -1] = 0.0

print("dx =", dx)
print("dy =", dy)
print("rho max =", np.max(rho))

Mass density distribution used as source term of the Poisson equation.

Numerical solution

# Parameters of the iterative method

max_iter = 20000
tol = 1e-7

# Solution of the Poisson equation with the Jacobi method

phi, n_iter, err = poisson_jacobi(
    phi0,
    rhs,
    dx,
    dy,
    max_iter=max_iter,
    tol=tol
)

print("Number of iterations =", n_iter)
print("Final error =", err)

if err < tol:
    print("The Jacobi method converged with the chosen tolerance.")
else:
    print("The method reached the maximum number of iterations.")

Numerical output.
dx = 0.01
dy = 0.01
rho max = 100.0
Number of iterations = 14862
Final error = 9.998478e-08

Gravitational potential — Computed gravitational potential. With the normalized sign convention used here, the positive mass source produces a negative potential well near the center while the boundary is fixed to zero.

Gravitational field

# Computation of the gravitational field from the potential

# np.gradient returns the derivative with respect to y first, then x
dphi_dy, dphi_dx = np.gradient(phi, dy, dx)

# Gravitational field: g = -grad(phi)
g_x = -dphi_dx
g_y = -dphi_dy