Numerical Solutions to Non-linear Partial Differential Equations with Code [Python]

In this article, we will delve into the numerical solutions of non-linear partial differential equations, with easy-to-understand concrete examples. We will learn an approach using Python.

Through this article, let’s master one of the numerical solutions for partial differential equations!

Note: The discussion is limited to the finite difference method. Details on errors and boundary conditions are omitted to avoid redundancy.

What is Numerical Solution of Partial Differential Equations?

The numerical solution of partial differential equations (PDEs) refers to methods used to approximate the solutions of PDEs. Since analytical solutions are often not possible, numerical methods become necessary. Below, we introduce several main numerical methods.

1. Finite Difference Method: This method discretizes space and time into a grid and approximates derivatives with differences. While intuitive and relatively easy to implement, the choice of grid significantly affects the solution’s accuracy.

2. Finite Element Method: The problem domain is divided into small “elements,” and equations are approximated within each element. This method is suitable for problems with complex shapes or boundary conditions.

3. Finite Volume Method: Particularly suited for conservation law equations, the domain is divided into small “volumes,” and average values within each volume are calculated.

4. Spectral Method: The equation is represented as a sum of waves with different frequencies, and solving involves calculating the coefficients of these waveforms. High accuracy can be achieved, but smooth solutions are required.

5. Finite-Difference Time-Domain (FDTD) Method: A type of finite difference method used extensively in electromagnetics, applied in the time domain.

These methods are used in physics, engineering, meteorology, and other fields to solve PDEs. Each method has its characteristics and applicable problem types, so choosing the optimal method according to the problem’s nature is necessary. In advanced computer simulations, these methods are often combined.

This article introduces the Finite Difference Method for numerical solutions.

Difference Approximation

In the numerical solution of partial differential equations, the finite difference method is very commonly used.

The finite difference method utilizes difference approximations. The details of the concept of difference approximation for partial derivatives are covered in another article, and here we will briefly introduce only the difference approximation used for parabolic equations, which is the main topic of this article.

Let $u = u(x,y)$ be given in the domain $0 \leqq x \leqq a, 0 \leqq y \leqq b$. Divide this domain into a grid (difference grid) as shown below. The grid spacing can be equal or, as shown in the figure, unequal.

Label the grid points in the $x$ direction as $i=0 \sim m$ and in the $y$ direction as $j=0 \sim n$. Label the point $(x,y)$ on the grid as $(i,j)$ and the value of $u$ at this point as $u_{i,j}$. Label the neighboring grid points as shown in the figure with $\Delta x_i,\Delta y_i$, etc. The difference approximations for partial derivatives $∂u/∂x, ∂u/∂y$ can be represented as follows, fixing $j$ for $∂u/∂x$ and $i$ for $∂u/∂y$.

When the grid spacing is uniform,

$$
\begin{aligned}
\Delta x_{i-1} &= \Delta x_{i} =h \\
\Delta y_{i-1} &= \Delta y_{i} =k \end{aligned}
$$

thus, the first-order partial derivative coefficients are obtained using the forward difference approximation as follows:

$$
\begin{aligned}
\frac{∂u}{∂x} &= \frac{1}{h}(u_{i+1,j} – u_{i,j}) \\
\frac{∂u}{∂y} &= \frac{1}{k}(u_{i,j+1} – u_{i,j})
\end{aligned}
$$

Furthermore, the second-order partial derivative coefficients are obtained using the central difference approximation as follows:

$$
\begin{align*}
\frac{∂^2u}{∂x∂y} &= \frac{1}{4hk}(u_{i+1,j+1} – u_{i+1,j-1} – u_{i-1,j+1} + u_{i-1,j-1}) \\
\frac{∂^2u}{∂x^2} &= \frac{1}{h^2}(u_{i+1,j} – 2u_{i,j} + u_{i-1,j}) \\
\frac{∂^2u}{∂y^2} &= \frac{1}{k^2}(u_{i,j+1} – 2u_{i,j} + u_{i,j-1})
\end{align*}
$$

These difference approximations are used to convert partial differential equations into difference equations.

Overview of Numerical Solution using the Finite Difference Method

The steps for numerically solving a partial differential equation using the finite difference method are:

1. Convert the partial differential equation into a difference equation.

2. Calculate the values of each element using the difference equation.

Let’s explain each step below.

Converting Partial Differential Equations to Difference Equations

Please bear with me as the introduction is a bit lengthy.

As an example of non-linear partial differential equations, let’s consider the case of a non-steady-state heat conduction equation where the conductivity $\lambda$ is a function of temperature $u(x,t)$, i.e., $\lambda=\lambda(u)$.
The non-steady-state heat conduction equation is expressed as

$$
c \rho \frac{\partial u}{\partial t} = \frac{\partial}{\partial x} \Big( \lambda \frac{\partial u}{\partial x} \Big) + \frac{\partial}{\partial y} \Big( \lambda \frac{\partial u}{\partial y} \Big) + \frac{\partial}{\partial z} \Big( \lambda \frac{\partial u}{\partial z} \Big)
$$

and in the one-dimensional case, it simplifies to

$$
\begin{align*} c \rho \frac{\partial u}{\partial t} &= \frac{\partial}{\partial x} \Big( \lambda \frac{\partial u}{\partial x} \Big) \\ \therefore \frac{\partial u}{\partial t} &= \frac{1}{c \rho} \Big( \frac{\partial \lambda}{\partial x} \frac{\partial u}{\partial x} + \lambda \frac{\partial^2 u}{\partial x^2} \Big) \end{align*}
$$

where $c, \rho$ are the specific heat and density, respectively. If we let $\lambda = au + b$ (where $a, b$ are constants),

$$
\frac{\partial u}{\partial t} = \frac{1}{c \rho} \Big\{ a \Big( \frac{\partial u}{\partial x}\Big)^2 + (au+b) \frac{\partial^2 u}{\partial x^2} \Big\} \ \ \ \ \ (1)
$$

Initial condition: $u(x, 0) = f(x)$

Boundary conditions: $u(0,t) = g_1(t), u(l,t) = g_2(t)$

with $0 \leqq x \leqq l, 0 \leqq t < \infty$.

Let’s convert this non-linear partial differential equation $(1)$ into a difference equation.
The conversion utilizes the difference approximations for partial derivatives as discussed in the previous section, namely

$$
\begin{align*} \frac{u_{i,j+1} – u_{i,j}}{k} &= \frac{1}{c \rho} \Big\{ \Big( \frac{u_{i+1,j} – u_{i-1,j} }{2h}\Big)^2 + (a u_{i,j}+b) \frac{(u_{i+1,j} – 2u_{i,j} + u_{i-1,j}) }{h^2} \Big\} \\ u_{i,j+1} &= u_{i,j} + \frac{k}{c \rho} \Big\{ \frac{a}{4h^2} ( u_{i+1,j} – u_{i-1,j})^2 + \frac{1}{h^2} (a u_{i,j}+b) (u_{i+1,j} – 2u_{i,j} + u_{i-1,j}) \Big\} \ \ \ \ (2) \end{align*}
$$

Calculating the Values of Each Element Using Difference Equations

Let’s take a closer look at difference equation $(2)$.

By setting $i=1,2, …, l-1$, the value of $u_{i,j}$ at $j=0$ is known from the initial and boundary conditions, allowing us to successively calculate the values of $u_{i,j}$ for times beyond $j=0$.

This method of directly obtaining the values of $u$ on the next time row from the known values of $u$ on the current time row is called the explicit method.

Let’s also reorganize difference equation $(2)$ as follows:

$$
\begin{align*} u_{i,j+1} &= u_{i,j} + \frac{ak}{4h^2} \frac{1}{c \rho} ( u_{i+1,j} – u_{i-1,j})^2 + \frac{k}{c \rho h^2} \lambda(u) (u_{i+1,j} – 2u_{i,j} + u_{i-1,j}) \\ &= \Big( 1 – \frac{2k}{h^2}\frac{1}{c \rho} \lambda(u) \Big)u_{i,j} + \frac{ak}{4h^2}\frac{1}{c \rho}(u_{i+1,j} – u_{i-1,j})^2 + \frac{k}{h^2}\frac{1}{c \rho}\lambda(u)(u_{i+1,j} + u_{i-1,j}) \ \ \ \ (3) \end{align*}
$$

Assuming all temperatures are positive, the second and third terms on the right-hand side of equation $(3)$ are positive, and the first term has the following condition:

$$
\begin{align*} 1 – \frac{2k}{h^2}\frac{1}{c \rho} \lambda(u) &\geqq 0 \\ \therefore \frac{k}{h^2}\frac{\lambda(u)}{c \rho} &\leqq \frac{1}{2} \ \ \ \ \ \ (4) \end{align*}
$$

While the explicit method is straightforward and easy to understand, the condition in equation $(4)$ necessitates taking $k$ to be very small relative to $h$, which can significantly increase computation time.

Even if $h$ and $k$ are chosen to satisfy equation $(4)$, since $\lambda(u)$ changes over time, it’s possible to eventually violate this condition. Therefore, it might be wise to program a check at each time step to ensure the condition is met and to adjust the value of $k$ or stop the calculation if it’s not. However, if the maximum value of $\lambda$, $\lambda_{max}$, can be estimated, then by setting $h$, $k$ can be determined from equation $(4)$ as

$$
k \leqq \frac{h^2 c \rho}{2 \lambda _{max}}
$$

to be effective for all $t$.

Example Problem

Background

When the thermal conductivity $\lambda$ is a function of temperature, the one-dimensional non-steady-state thermoelectric equation is represented as follows:

$$
\frac{\partial u}{\partial t} = \frac{1}{c \rho} \Big( \frac{\partial \lambda}{\partial x} \frac{\partial u}{\partial x} + \lambda \frac{\partial^2 u}{\partial x^2} \Big)
$$

Now, by non-dimensionalizing this equation using a reference temperature $u_0$, reference length $l_0$, and reference thermal conductivity $\lambda_0$, we get:

$$
\frac{\partial U}{\partial \tau} = \frac{\partial R}{\partial X} \frac{\partial U}{\partial X} + R \frac{\partial^2 U}{\partial X^2} \ \ \ \ (5)
$$

Where,

$$
U=u/u_0, \ \ \ X=x/l_0, \ \ \ R=\lambda/\lambda_0, \ \ \ \tau=\lambda_0 t/ c \rho l_0^2 \ \ \ \ (6)
$$

These are all non-dimensional quantities representing non-dimensional temperature, non-dimensional coordinates, non-dimensional thermal conductivity, and non-dimensional time, respectively.

Problem

Consider the following scenario as the problem for equation $(5)$: The temperature $U(x, \tau)$ of a thin rod was initially $0$. When the temperature at both ends of the rod is maintained at $U=1+\sin2\pi \tau$, investigate the temporal change in the temperature at the center of the rod, assuming the length of the rod is $1$ and there is no heat transfer from the surroundings of the rod. The non-dimensional thermal conductivity $R$ is given as a linear function of $U$

$$
R=U+1 \ \ \ \ (7)
$$

Hint

The fundamental equation and conditions for the temperature distribution $U(X,\tau)$ are as follows:

$$
\frac{\partial U}{\partial \tau} = \Big( \frac{\partial U}{\partial X}\Big)^2 + (U + 1)\frac{\partial^2 U}{\partial X^2}; \ \ \ \ \ \ \ 0 \leqq X \leqq 1, 0 \leqq \tau < \infty
$$

Initial condition: $U(X, 0) = 0$

Boundary conditions: $U(0,\tau) = U(1,\tau) = 1 + \sin 2 \pi \tau$

Let the grid numbers be $i=0 \sim m, j=0 \sim n$, and the value of $U$ at any grid point be $U_{i,j}$.

As introduced earlier, let’s first convert this partial differential equation into a difference equation, and then implement a program to solve the difference equation.

Implementation Method

The code can be implemented as follows:

The creation of graphs and output to text files is also included, making the code somewhat lengthy, but the crucial part is the calculate_equation function.

Save the following implementation as main.py in a folder (here, ~/Desktop/labcode/python/numerical_calculation/non_linear_equation).

import matplotlib.pyplot as plt
from matplotlib.animation import ArtistAnimation
import numpy as np

def _set_initial_condition(
    *,
    array_2d: list,
    grid_counts_x: int,
) -> None:
    """
    Set initial conditions

    Set all U values except for the boundaries to 0.0
    """
    for i in range(1, grid_counts_x):
        array_2d[i][0] = 0.0

def _set_boundary_condition(
    *,
    array_2d: list,
    grid_counts_x: int,
    grid_counts_t: int,
    time_delta: float,
) -> None:
    """Set boundary conditions"""
    for j in range(grid_counts_t + 1):
        array_2d[0][j] = 1 + np.sin(2 * np.pi * j * time_delta)
        array_2d[grid_counts_x][j] = 1 + np.sin(2 * np.pi * j * time_delta)

def calculate_equation(*, M: int, N: int, H: float, K: float) -> list:
    """
    Calculate the difference equation

    The core function of this program
    """
    # List for U values
    U: list = [[0.0 for _ in range(N + 1)] for _ in range(M + 1)]

    # Set initial values
    _set_initial_condition(
        array_2d=U,
        grid_counts_x=M,
    )
    _set_boundary_condition(
        array_2d=U,
        grid_counts_x=M,
        grid_counts_t=N,
        time_delta=K,
    )

    # Perform calculation
    KH = K / (H**2)
    K4 = KH / 4
    for j in range(N):
        for i in range(1, M):
            A = U[i][j] + K4 * (U[i + 1][j] - U[i - 1][j]) ** 2
            B = KH * (U[i][j] + 1) * (U[i + 1][j] - 2 * U[i][j] + U[i - 1][j])
            U[i][j + 1] = A + B
    return U

def animate_time_evolution(
    *,
    array_2d: list,
    grid_counts_x: int,
    grid_counts_t: int,
    grid_space: float,
    time_delta: float,
) -> None:
    """
    Create an animation of the dot-line plot
    Plot the data in array_1d and show the time evolution in an animation
    """
    # Transpose 2D array
    array_2d = np.array(array_2d).T
    # Generate x-axis values
    x_values = np.arange(grid_counts_x + 1) * grid_space

    # Create the animation
    path_gif = "./animation.gif"
    frames = []

    fig, ax = plt.subplots(figsize=(7, 5))
    plt.xlabel("X (a.u.)")
    plt.ylabel("Temperature (a.u.)")
    plt.grid(True)

    for j in range(grid_counts_t + 1):
        frame = plt.plot(
            x_values,
            array_2d[j],
            linestyle="--",
            marker="o",
            color="b",
        )

        time = j * time_delta
        text = ax.text(
            0.05,
            1.05,
            f"tau = {time:.3f} (a.u.)",
            transform=ax.transAxes,
            fontsize=14,
            verticalalignment="top",
        )
        frames.append(frame + [text])

    anim = ArtistAnimation(fig, frames, interval=100)
    anim.save(path_gif, writer="pillow")

def output_result_file(
    array_2d: list,
    grid_counts_x: int,
    grid_counts_t: int,
    grid_space: float,
    time_delta: float,
) -> None:
    """Write a 2D array to a text file"""
    with open("./calculated_result.txt", "w") as file:
        file.write(f"# This file shows calculated result as below.\n\n")
        file.write(f"# Calculation Parameters.\n")
        file.write(f"grid_counts_x: {grid_counts_x}\n")
        file.write(f"grid_counts_t: {grid_counts_t}\n")
        file.write(f"grid_space: {grid_space}\n")
        file.write(f"time_delta: {time_delta}\n\n")

        file.write(f"# Calculated Matrix U.\n")
        for row in array_2d:
            line = " ".join(map(str, row))
            file.write(line + "\n")

def validate_parameters(
    grid_counts_t: int, grid_space: float, time_delta: float
) -> None:
    """Check the validity of parameters"""

    def get_max_u_value() -> float:
        max_u = 0.0
        for j in range(grid_counts_t + 1):
            u = 1 + np.sin(2 * np.pi * j * time_delta)
            if max_u < u:
                max_u = u
        return max_u

    max_r_value = get_max_u_value() + 1
    r = (grid_space**2) / (2 * max_r_value)

    if r > 0.5:
        raise ValueError("R must be less than 0.5.")

if __name__ == "__main__":
    """
    Numerical solution of non-linear equations (explicit method)
    """
    grid_counts_x: int = 10  # Value of grid point number m
    grid_counts_t: int = 1000  # Value of grid point number n
    grid_space: float = 0.1  # Step size (a.u.)
    time_delta: float = 0.001  # Step size (a.u.)

    validate_parameters(
        grid_counts_t=grid_counts_t,
        grid_space=grid_space,
        time_delta=time_delta,
    )

    array_2d = calculate_equation(
        M=grid_counts_x,
        N=grid_counts_t,
        H=grid_space,
        K=time_delta,
    )
    animate_time_evolution(
        array_2d=array_2d,
        grid_counts_x=grid_counts_x,
        grid_counts_t=grid_counts_t,
        grid_space=grid_space,
        time_delta=time_delta,
    )
    output_result_file(
        array_2d=array_2d,
        grid_counts_x=grid_counts_x,
        grid_counts_t=grid_counts_t,
        grid_space=grid_space,
        time_delta=time_delta,
    )

Executing the Program

Open the terminal and enter:

$ cd Desktop/labcode/python/numerical_calculation/non_linear_equation

Then, execute the following command (ignore the $ symbol):

$ python main.py

Execution Results

If a file named animation.gif has been created in ~/Desktop/LabCode/python/numerical_calculation/non_linear_equation and it turns into a gif image like the following, then it’s a success!

The visualization shows how the temperature in the inner regions changes in response to the gradual changes in temperature at both ends.

Additionally, a text file named calculated_result.txt should have been created as follows. Although only basic information is outputted this time, by outputting the calculation results along with the parameters used for the calculation, it becomes possible for others to evaluate the results. For your own benefit as well as for the benefit of others, make it a habit to output the results.

# This file shows calculated result as below.

# Calculation Parameters.
grid_counts_x: 10
grid_counts_t: 1000
grid_space: 0.1
time_delta: 0.001

# Calculated Matrix U.
1.0 1.0062831439655588 1.0125660398833527 1.0188484397154083 1.0251300954433376 1.0314107590781283 1.0376901826699345 1.0439681183178648 1.0502443181797696 1.0565185344820245 1.0627905195293135 1.0690600257144058 1.0753268055279328 1.0815906115681575 1.
...

In Conclusion

We have introduced the numerical solution of non-linear partial differential equations using the finite difference method.

The Python computation program is also generously shared, so please deepen your understanding by running it.

In the future, we will write new articles on numerically solving other types of partial differential equations, so we hope you will look forward to them.