Heat equation in 1D

To verify our Laplacian calculation, we solve the heat equation in one dimension with zero flux at the boundaries. The initial condition is a cosine, which yields the trivial analytical solution:

\(f(x, t) = \cos(2\pi x)e^{-4\pi^2 t}\).

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

x, dx = np.linspace(0, 1, 60, retstep=True)
dt = 0.5 * dx**2
print(f"{dt = :.3e}")

# Initial condition
n_periods = 1
f_0 = np.cos(2 * n_periods * np.pi * x)


def f_ana(t: float) -> np.ndarray:
    """
    Retunrs the analytical solution:
    f(x, t) = cos(2πx) * exp(-4π²t)
    """
    return np.exp(-(2 * n_periods * np.pi)**2 * t) * f_0

dt = 1.436e-04

A temperature update function is created.

def update(f):
    """Met à jour la solution numérique en calculant le laplacien 1D"""
    f_moins = np.concatenate((f[1:2], f[:-1]))
    f_plus = np.concatenate((f[1:], f[-2:-1]))
    laplacien = (f_moins - 2 * f + f_plus) / dx**2
    f += dt * laplacien

We define a time generator.

def time_generator(N):
    """
    Générateur de temps qui s'arrête à N ou quand t vaut
    1 / (2 * n_periodes * np.pi)**2
    """
    for i in range(N):
        t = i * dt
        if t > 1 / (2 * n_periods * np.pi)**2:
            return
        else:
            yield t

Initialize the figure and run the animation.

%matplotlib notebook

fig = plt.figure()
N = 5000

f = f_0.copy()
line, = plt.plot(x, f, label="numerical")
ana, = plt.plot(x, f_ana(0.), '+', label=r'$\cos(2\pi x)e^{-4\pi^2 t}$')
ax = fig.get_axes()[0]
ax.legend()
ax.set_xlabel("x")


def animate(t):

    # Flux is zero at the edges, so the integral must remain constant
    # We calculate the integral using the midpoint method
    integral = (f[1:-1].sum() + 0.5 * (f[0] + f[-1])) * dx
    error = np.linalg.norm(f - f_ana(t), ord=2)
    title = f"{t = :.3e} {integral = :.3e} {error = :.3e}"
    fig.suptitle(title)

    update(f)

    line.set_ydata(f)
    ana.set_ydata(f_ana(t))

    return line, ana


ani = animation.FuncAnimation(fig, animate, frames=time_generator(N),
                              interval=1, blit=True, repeat=False,
                              cache_frame_data=False)

from IPython.display import HTML
HTML(ani.to_jshtml())

Once Loop Reflect