Skip to content

Instantly share code, notes, and snippets.

@Rodbourn
Created August 25, 2024 23:53
Show Gist options
  • Select an option

  • Save Rodbourn/cd66ce14afd4e7768fb3ff4606524c6f to your computer and use it in GitHub Desktop.

Select an option

Save Rodbourn/cd66ce14afd4e7768fb3ff4606524c6f to your computer and use it in GitHub Desktop.
import numpy as np
import matplotlib.pyplot as plt
def heat_equation_1d(L, nx, nt, alpha, T_initial, T_left, T_right):
"""
Solve the 1D heat equation using explicit finite difference method.
Parameters:
L : float
Length of the domain
nx : int
Number of spatial points
nt : int
Number of time steps
alpha : float
Thermal diffusivity
T_initial : float
Initial temperature
T_left : float
Left boundary temperature
T_right : float
Right boundary temperature
Returns:
T : ndarray
Temperature distribution over time and space
"""
# Calculate step sizes
dx = L / (nx - 1)
dt = 0.5 * dx**2 / alpha # Ensure stability
# Initialize temperature array
T = np.ones((nt, nx)) * T_initial
# Set boundary conditions
T[:, 0] = T_left
T[:, -1] = T_right
# Compute finite difference
for n in range(1, nt):
for i in range(1, nx-1):
T[n, i] = T[n-1, i] + alpha * dt / dx**2 * \
(T[n-1, i+1] - 2*T[n-1, i] + T[n-1, i-1])
return T
# Set parameters
L = 1.0 # Length of the domain
nx = 50 # Number of spatial points
nt = 1000 # Number of time steps
alpha = 0.01 # Thermal diffusivity
T_initial = 0 # Initial temperature
T_left = 100 # Left boundary temperature
T_right = 0 # Right boundary temperature
# Solve the heat equation
T = heat_equation_1d(L, nx, nt, alpha, T_initial, T_left, T_right)
# Plot results
x = np.linspace(0, L, nx)
plt.figure(figsize=(10, 6))
plt.plot(x, T[0, :], label='t = 0')
plt.plot(x, T[int(nt/4), :], label=f't = {int(nt/4)}')
plt.plot(x, T[int(nt/2), :], label=f't = {int(nt/2)}')
plt.plot(x, T[-1, :], label=f't = {nt-1}')
plt.xlabel('Position')
plt.ylabel('Temperature')
plt.title('1D Heat Equation - Finite Difference Solution')
plt.legend()
plt.grid(True)
plt.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment