Partial Differential Equation (PDE)

I’ll help you solve a Partial Differential Equation (PDE) example using $Python$.

Let’s solve a 2D heat equation, which is a classic PDE problem.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def solve_2d_heat_equation(dx=0.1, dy=0.1, dt=0.0001, t_final=0.01):
"""
Solve 2D heat equation using finite difference method
∂u/∂t = α(∂²u/∂x² + ∂²u/∂y²)
with initial conditions: u(x,y,0) = sin(πx)sin(πy)
and boundary conditions: u = 0 at boundaries
"""
# Thermal diffusivity
alpha = 1.0

# Create grid
x = np.arange(0, 1+dx, dx)
y = np.arange(0, 1+dy, dy)
t = np.arange(0, t_final+dt, dt)

# Initialize temperature grid
nx, ny = len(x), len(y)
u = np.zeros((nx, ny))

# Set initial conditions
X, Y = np.meshgrid(x, y)
u = np.sin(np.pi*X) * np.sin(np.pi*Y)

# Create arrays for new and old solutions
u_new = np.zeros((nx, ny))

# Stability condition check
stability = alpha * dt * (1/(dx**2) + 1/(dy**2))
print(f"Stability parameter: {stability:.3f} (should be < 0.5)")

# Time stepping
for n in range(len(t)-1):
# Interior points
for i in range(1, nx-1):
for j in range(1, ny-1):
u_new[i,j] = u[i,j] + alpha*dt*(
(u[i+1,j] - 2*u[i,j] + u[i-1,j])/dx**2 +
(u[i,j+1] - 2*u[i,j] + u[i,j-1])/dy**2
)

# Update solution
u = u_new.copy()

# Enforce boundary conditions
u[0,:] = u[-1,:] = u[:,0] = u[:,-1] = 0

return X, Y, u

# Solve the equation
X, Y, u_final = solve_2d_heat_equation()

# Create visualization
fig = plt.figure(figsize=(12, 5))

# Surface plot
ax1 = fig.add_subplot(121, projection='3d')
surf = ax1.plot_surface(X, Y, u_final, cmap=cm.viridis)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('Temperature')
ax1.set_title('2D Heat Distribution (3D View)')
fig.colorbar(surf)

# Contour plot
ax2 = fig.add_subplot(122)
contour = ax2.contourf(X, Y, u_final, levels=20, cmap=cm.viridis)
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_title('2D Heat Distribution (Top View)')
fig.colorbar(contour)

plt.tight_layout()
plt.show()

Let me explain this code and the problem we’re solving:

  1. The Problem: We’re solving the 2D heat equation:

    • Equation: $∂u/∂t = α(∂²u/∂x² + ∂²u/∂y²)$
    • Initial condition: $u(x,y,0) = sin(πx)sin(πy)$
    • Boundary conditions: $u = 0$ at all boundaries
    • $α$ is the thermal diffusivity (set to $1.0$ in our case)
  2. Solution Method:

    • We use the Finite Difference Method (FDM)
    • Space is discretized into a grid with steps $dx$ and $dy$
    • Time is discretized with step $dt$
    • The solution uses an explicit scheme for time-stepping
  3. Code Structure:

    • Creates a mesh grid for $x$ and $y$ coordinates
    • Sets initial conditions using sine functions
    • Implements time-stepping loop with boundary conditions
    • Uses $numpy$ for efficient array operations
  4. Visualization:

    • Left plot: 3D surface plot showing temperature distribution
    • Right plot: 2D contour plot showing temperature from top view
    • Colors indicate temperature (darker blue = cooler, yellow = warmer)
Stability parameter: 0.020 (should be < 0.5)

  1. Results Interpretation:
    • The initial sine wave pattern gradually diffuses
    • Heat spreads from warmer regions to cooler regions
    • The solution maintains zero temperature at boundaries
    • The pattern becomes smoother over time due to diffusion

The visualization shows both a 3D surface plot and a 2D contour plot to help understand how heat distributes across the 2D surface.

The color gradient represents temperature variations, with warmer colors indicating higher temperatures.