xbe.at

Data visualisation

Visualizing Partial Differential Equations with Python

9 min read, 16 sections

Contents
  1. Introduction to Partial Differential Equations (PDEs)
  2. Types of PDEs
  3. Elliptic PDEs
  4. Parabolic PDEs
  5. Hyperbolic PDEs
  6. Finite Difference Method
  7. Finite Element Method
  8. Method of Characteristics
  9. Spectral Methods
  10. Numerical Stability and Convergence
  11. Boundary Value Problems
  12. Initial Value Problems
  13. Real-life Example: Heat Conduction in a Rod
  14. Real-life Example: Wave Propagation in a String
  15. Numerical Methods for Nonlinear PDEs
  16. Additional Resources

Introduction to Partial Differential Equations (PDEs)

Partial Differential Equations (PDEs) are equations involving functions of multiple variables and their partial derivatives. They are essential in modeling various physical phenomena.

import sympy as sp

# Define variables and function
x, y, t = sp.symbols('x y t')
u = sp.Function('u')(x, y, t)

# Example of a PDE: Heat equation
heat_equation = sp.Eq(sp.diff(u, t), sp.diff(u, x, x) + sp.diff(u, y, y))
print(f"Heat Equation: {heat_equation}")

Types of PDEs

PDEs are classified based on their order and linearity. Common types include:

  1. First-order PDEs
  2. Second-order PDEs (e.g., elliptic, parabolic, hyperbolic)
# First-order PDE example: Transport equation
transport_equation = sp.Eq(sp.diff(u, t) + sp.diff(u, x), 0)
print(f"Transport Equation: {transport_equation}")

# Second-order PDE example: Wave equation
wave_equation = sp.Eq(sp.diff(u, t, t), sp.diff(u, x, x))
print(f"Wave Equation: {wave_equation}")

Elliptic PDEs

Elliptic PDEs describe steady-state phenomena. A classic example is Laplace's equation.

import numpy as np
import matplotlib.pyplot as plt

def laplace_equation(nx, ny, iterations):
    u = np.zeros((ny, nx))
    u[0, :] = 100  # Top boundary condition
    
    for _ in range(iterations):
        u[1:-1, 1:-1] = 0.25 * (u[1:-1, 2:] + u[1:-1, :-2] + u[2:, 1:-1] + u[:-2, 1:-1])
    
    return u

nx, ny = 50, 50
u = laplace_equation(nx, ny, 1000)
plt.imshow(u, cmap='hot')
plt.colorbar()
plt.title("Solution to Laplace's Equation")
plt.show()

Parabolic PDEs

Parabolic PDEs model diffusion processes. The heat equation is a prime example.

def heat_equation_1d(nx, nt, alpha):
    dx = 1 / (nx - 1)
    dt = 0.5 * dx**2 / alpha
    u = np.zeros(nx)
    u[0] = 100  # Left boundary condition
    
    for _ in range(nt):
        u[1:-1] += alpha * dt / dx**2 * (u[2:] - 2*u[1:-1] + u[:-2])
    
    return u

nx, nt = 50, 1000
alpha = 0.01
u = heat_equation_1d(nx, nt, alpha)
plt.plot(np.linspace(0, 1, nx), u)
plt.title("1D Heat Equation Solution")
plt.xlabel("Position")
plt.ylabel("Temperature")
plt.show()

Hyperbolic PDEs

Hyperbolic PDEs describe wave-like phenomena. The wave equation is a classic example.

def wave_equation_1d(nx, nt, c):
    dx = 1 / (nx - 1)
    dt = 0.5 * dx / c
    u = np.zeros((3, nx))
    u[0, nx//2] = 1  # Initial pulse
    
    for n in range(1, nt):
        u[2, 1:-1] = 2*u[1, 1:-1] - u[0, 1:-1] + c**2 * dt**2 / dx**2 * (u[1, 2:] - 2*u[1, 1:-1] + u[1, :-2])
        u[0], u[1] = u[1], u[2]
    
    return u[1]

nx, nt = 100, 200
c = 1
u = wave_equation_1d(nx, nt, c)
plt.plot(np.linspace(0, 1, nx), u)
plt.title("1D Wave Equation Solution")
plt.xlabel("Position")
plt.ylabel("Amplitude")
plt.show()

Finite Difference Method

The finite difference method approximates derivatives using discrete differences.

def central_difference(f, x, h):
    return (f(x + h) - f(x - h)) / (2 * h)

def f(x):
    return np.sin(x)

x = np.pi / 4
h = 0.001
approx_derivative = central_difference(f, x, h)
exact_derivative = np.cos(x)

print(f"Approximate derivative: {approx_derivative}")
print(f"Exact derivative: {exact_derivative}")
print(f"Error: {abs(approx_derivative - exact_derivative)}")

Finite Element Method

The finite element method divides the domain into smaller elements and approximates the solution using basis functions.

from scipy.integrate import quad

def basis_function(x, xi, h):
    if xi - h <= x <= xi:
        return (x - (xi - h)) / h
    elif xi < x <= xi + h:
        return ((xi + h) - x) / h
    else:
        return 0

def assemble_system(n, f):
    h = 1 / n
    A = np.zeros((n-1, n-1))
    b = np.zeros(n-1)
    
    for i in range(n-1):
        for j in range(n-1):
            xi, xj = (i+1)*h, (j+1)*h
            A[i, j] = quad(lambda x: basis_function(x, xi, h) * basis_function(x, xj, h), 0, 1)[0]
        b[i] = quad(lambda x: f(x) * basis_function(x, (i+1)*h, h), 0, 1)[0]
    
    return A, b

# Example usage
n = 10
f = lambda x: np.sin(np.pi * x)
A, b = assemble_system(n, f)
u = np.linalg.solve(A, b)

x = np.linspace(0, 1, 100)
plt.plot(x, f(x), label='Exact')
plt.plot(np.linspace(h, 1-h, n-1), u, 'o-', label='FEM')
plt.legend()
plt.title("Finite Element Method Solution")
plt.show()

Method of Characteristics

The method of characteristics solves first-order PDEs by finding characteristic curves.

def method_of_characteristics(u0, a, T, nx, nt):
    x = np.linspace(0, 1, nx)
    t = np.linspace(0, T, nt)
    u = np.zeros((nt, nx))
    
    u[0, :] = u0(x)
    for i in range(1, nt):
        u[i, :] = u[0, (x - a * t[i]) % 1]
    
    return u, x, t

def initial_condition(x):
    return np.sin(2 * np.pi * x)

a = 1
T = 1
nx, nt = 100, 50
u, x, t = method_of_characteristics(initial_condition, a, T, nx, nt)

plt.imshow(u, extent=[0, 1, 0, T], aspect='auto', origin='lower')
plt.colorbar()
plt.title("Solution using Method of Characteristics")
plt.xlabel("x")
plt.ylabel("t")
plt.show()

Spectral Methods

Spectral methods use Fourier series or other orthogonal functions to approximate solutions.

def spectral_method(u0, N, T, nt):
    x = np.linspace(0, 2*np.pi, N, endpoint=False)
    t = np.linspace(0, T, nt)
    u = np.zeros((nt, N))
    
    u[0] = u0(x)
    u_hat = np.fft.fft(u[0])
    k = np.fft.fftfreq(N, 2*np.pi/N)
    
    for i in range(1, nt):
        u_hat *= np.exp(-1j * k**2 * t[i])
        u[i] = np.real(np.fft.ifft(u_hat))
    
    return u, x, t

def initial_condition(x):
    return np.exp(-(x - np.pi)**2)

N = 128
T = 0.1
nt = 50
u, x, t = spectral_method(initial_condition, N, T, nt)

plt.imshow(u, extent=[0, 2*np.pi, 0, T], aspect='auto', origin='lower')
plt.colorbar()
plt.title("Solution using Spectral Method")
plt.xlabel("x")
plt.ylabel("t")
plt.show()

Numerical Stability and Convergence

Numerical stability and convergence are crucial for accurate PDE solutions.

def von_neumann_stability(scheme, dx, dt, k):
    return abs(scheme(k * dx, dt))

def ftcs_scheme(kdx, dt):
    return 1 - 2 * np.sin(kdx/2)**2

dx = 0.1
dt = 0.001
k_values = np.linspace(0, np.pi/dx, 100)
stability = [von_neumann_stability(ftcs_scheme, dx, dt, k) for k in k_values]

plt.plot(k_values, stability)
plt.axhline(y=1, color='r', linestyle='--')
plt.title("Von Neumann Stability Analysis")
plt.xlabel("k*dx")
plt.ylabel("|G(k)|")
plt.show()

Boundary Value Problems

Boundary value problems involve solving PDEs with specified boundary conditions.

def solve_bvp(a, b, c, f, ya, yb, N):
    h = (b - a) / (N - 1)
    x = np.linspace(a, b, N)
    
    A = np.zeros((N, N))
    B = np.zeros(N)
    
    A[0, 0] = 1
    B[0] = ya
    A[-1, -1] = 1
    B[-1] = yb
    
    for i in range(1, N-1):
        A[i, i-1] = 1 / h**2
        A[i, i] = -2 / h**2 + c(x[i])
        A[i, i+1] = 1 / h**2
        B[i] = f(x[i])
    
    y = np.linalg.solve(A, B)
    return x, y

a, b = 0, 1
c = lambda x: 0
f = lambda x: np.exp(-x)
ya, yb = 0, 1
N = 100

x, y = solve_bvp(a, b, c, f, ya, yb, N)
plt.plot(x, y)
plt.title("Solution to Boundary Value Problem")
plt.xlabel("x")
plt.ylabel("y")
plt.show()

Initial Value Problems

Initial value problems involve solving PDEs with specified initial conditions.

def solve_ivp(f, y0, t_span, dt):
    t = np.arange(t_span[0], t_span[1], dt)
    y = np.zeros((len(t), len(y0)))
    y[0] = y0
    
    for i in range(1, len(t)):
        k1 = dt * f(t[i-1], y[i-1])
        k2 = dt * f(t[i-1] + dt/2, y[i-1] + k1/2)
        k3 = dt * f(t[i-1] + dt/2, y[i-1] + k2/2)
        k4 = dt * f(t[i-1] + dt, y[i-1] + k3)
        y[i] = y[i-1] + (k1 + 2*k2 + 2*k3 + k4) / 6
    
    return t, y

def f(t, y):
    return np.array([y[1], -np.sin(y[0])])

y0 = [np.pi/2, 0]
t_span = [0, 10]
dt = 0.01

t, y = solve_ivp(f, y0, t_span, dt)
plt.plot(t, y[:, 0])
plt.title("Solution to Initial Value Problem")
plt.xlabel("t")
plt.ylabel("y")
plt.show()

Real-life Example: Heat Conduction in a Rod

Let's model heat conduction in a metal rod using the heat equation.

def heat_conduction_rod(L, T, nx, nt, alpha, T_left, T_right):
    dx = L / (nx - 1)
    dt = T / (nt - 1)
    x = np.linspace(0, L, nx)
    t = np.linspace(0, T, nt)
    
    u = np.zeros((nt, nx))
    u[0, :] = 0  # Initial temperature
    u[:, 0] = T_left  # Left boundary condition
    u[:, -1] = T_right  # Right boundary condition
    
    for j in range(1, nt):
        for i in range(1, nx-1):
            u[j, i] = u[j-1, i] + alpha * dt / dx**2 * (u[j-1, i+1] - 2*u[j-1, i] + u[j-1, i-1])
    
    return u, x, t

L = 1  # Length of the rod
T = 0.5  # Total time
nx, nt = 50, 1000
alpha = 0.01  # Thermal diffusivity
T_left, T_right = 100, 0  # Boundary temperatures

u, x, t = heat_conduction_rod(L, T, nx, nt, alpha, T_left, T_right)

plt.imshow(u, extent=[0, L, 0, T], aspect='auto', origin='lower', cmap='hot')
plt.colorbar(label='Temperature')
plt.title("Heat Conduction in a Rod")
plt.xlabel("Position (m)")
plt.ylabel("Time (s)")
plt.show()

Real-life Example: Wave Propagation in a String

Let's model wave propagation in a vibrating string using the wave equation.

def wave_propagation_string(L, T, nx, nt, c):
    dx = L / (nx - 1)
    dt = T / (nt - 1)
    x = np.linspace(0, L, nx)
    t = np.linspace(0, T, nt)
    
    u = np.zeros((nt, nx))
    u[0, :] = np.sin(np.pi * x / L)  # Initial displacement
    u[1, :] = u[0, :]  # Initial velocity is zero
    
    r = c * dt / dx
    for j in range(1, nt-1):
        for i in range(1, nx-1):
            u[j+1, i] = 2 * u[j, i] - u[j-1, i] + r**2 * (u[j, i+1] - 2*u[j, i] + u[j, i-1])
    
    return u, x, t

L = 1  # Length of the string
T = 2  # Total time
nx, nt = 100, 200
c = 1  # Wave speed

u, x, t = wave_propagation_string(L, T, nx, nt, c)

plt.imshow(u, extent=[0, L, 0, T], aspect='auto', origin='lower', cmap='viridis')
plt.colorbar(label='Displacement')
plt.title("Wave Propagation in a String")
plt.xlabel("Position (m)")
plt.ylabel("Time (s)")
plt.show()

Numerical Methods for Nonlinear PDEs

Nonlinear PDEs often require specialized numerical methods. Let's look at the Burgers' equation as an example.

def burgers_equation(nx, nt, Lx, Lt, nu):
    dx = Lx / (nx - 1)
    dt = Lt / (nt - 1)
    x = np.linspace(0, Lx, nx)
    t = np.linspace(0, Lt, nt)
    
    u = np.zeros((nt, nx))
    u[0, :] = np.sin(2 * np.pi * x / Lx)  # Initial condition
    
    for n in range(nt - 1):
        u[n+1, 0] = u[n+1, -1] = 0  # Boundary conditions
        for i in range(1, nx - 1):
            u[n+1, i] = u[n, i] - u[n, i] * dt / dx * (u[n, i] - u[n, i-1]) + \
                        nu * dt / dx**2 * (u[n, i+1] - 2 * u[n, i] + u[n, i-1])
    
    return u, x, t

nx, nt = 100, 500
Lx, Lt = 2, 1
nu = 0.01  # Viscosity

u, x, t = burgers_equation(nx, nt, Lx, Lt, nu)

plt.imshow(u, extent=[0, Lx, 0, Lt], aspect='auto', origin='lower', cmap='coolwarm')
plt.colorbar(label='Velocity')
plt.title("Solution of Burgers' Equation")
plt.xlabel("Position")
plt.ylabel("Time")
plt.show()

Additional Resources

For further study on Partial Differential Equations, consider these resources:

  1. ArXiv.org papers:
    • "Numerical Methods for Partial Differential Equations" by J. S. Hesthaven (arXiv:1803.07454)
    • "Machine Learning for Partial Differential Equations" by M. Raissi et al. (arXiv:1706.02242)
  2. Online courses:
    • MIT OpenCourseWare: "Partial Differential Equations" (18.152)
    • Coursera: "Numerical Methods for Partial Differential Equations" by The Hong Kong University of Science and Technology
  3. Textbooks:
    • "Partial Differential Equations for Scientists and Engineers" by Stanley J. Farlow
    • "Numerical Solution of Partial Differential Equations" by K. W. Morton and D. F. Mayers

Remember to verify the accuracy and relevance of these resources, as they may have been updated or changed since my last knowledge update.

More in Data visualisation

  1. Visualizing Pride and Prejudice with Python12 min

    "Pride and Prejudice" is Jane Austen's beloved novel about the Bennet family in Regency-era England. This presentation will explore key events and themes...

  2. Visualizing Machine Learning Insights8 min

    Visualization is a crucial tool in machine learning that helps us understand complex data, interpret models, and gain insights into the learning process. It...

  3. Visualizing R-Squared and Regression Variance with Python7 min

    R² and regression variance are fundamental concepts in statistical analysis. This presentation explores these concepts using Euler/Venn diagrams, providing a...

  4. Visualizing K-Means Clustering with Animation8 min

    K-Means is a popular unsupervised machine learning algorithm used for clustering data points into groups based on similarity. This animation explores the...

  5. Visualizing Riemannian Manifolds with Python11 min

    Riemannian manifolds are fundamental structures in differential geometry, combining the concepts of smooth manifolds with inner product spaces. They provide...

  6. Visualizing Jordan Canonical Form in Python8 min

    Jordan Canonical Form is a powerful concept in linear algebra that allows us to represent a matrix in a simplified, block-diagonal structure. This form is...