I'm currently writing a Python script that'll help me simulate the heat equation in a 2D rectangular room. In order to solve this problem, I had to learn about the Finite Difference Method. By simulating how the heat in the room changes given some boundary conditions on the temperature of the walls, I'm essentially plotting the surface of the temperature inside a loop, clearing the graph and then replotting in the next frame. Here's my code
# Libraries
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
import matplotlib as mpl
def f(X,Y): # f(x,y,0)
return 3 - (X**2+Y**2) * np.sin(Y)
def simulation_heat_eq(rect,hs,BC, c, frames, eps = 1e-10):
""" Simulation of the heat equation on a
2D rectangular grid using the Finite Difference Method """
# Figure settings
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.set_facecolor((0.8,0.8,0.8))
X, Y = np.meshgrid(np.linspace(rect[0],rect[1],hs[0]),
np.linspace(rect[2],rect[3],hs[1])) # 2D meshgrid
# Initial Conditions
Z_init = f(X,Y) # f(x,y,0)
zmax = 3
# Boundary Conditions
Z_init[0] = np.ones(len(Z_init[0])) * BC[0]
Z_init[-1] = np.ones(len(Z_init[-1]))* BC[1]
Z_init[:,0] = np.ones(len(Z_init[:,0]))* BC[2]
Z_init[:, -1] = np.ones(len(Z_init[:, -1])) * BC[3]
# Infitesimals
dx = (rect[1] - rect[0]) / (hs[0] - 1)
dt = dx**2/(10*c**2)
# Looping over Z:s
for iteration in range(frames):
# Clearing graph
ax.cla()
ax.axes.set_xlim3d(rect[0] + eps, rect[1] - eps)
ax.axes.set_ylim3d(rect[2] + eps, rect[3] - eps)
ax.axes.set_zlim3d(0, zmax)
ax.set_axis_off()
# Finite Difference Method
for row in range(1, hs[0]-1):
for col in range(1,hs[1]-1):
Z_init[row,col] = Z_init[row,col] + c**2 * dt / (dx)**2 *(Z_init[row+1,col] - 4*Z_init[row,col]
+ Z_init[row-1,col] + Z_init[row,col+1]
+ Z_init[row,col-1])
# Plotting surface with slight pause
#if iteration > 0:
# ax.view_init(10+5*np.sin(iteration*np.pi/180), iteration/3)
ax.plot([rect[0], rect[1]], [rect[2], rect[2]], [BC[0], BC[0]], color='black', linewidth=2)
ax.plot([rect[1], rect[1]], [rect[2], rect[3]], [BC[3], BC[3]], color='black', linewidth=2)
ax.plot([rect[0], rect[0]], [rect[2], rect[3]], [BC[2], BC[2]], color='black', linewidth=2)
ax.plot([rect[0], rect[1]], [rect[3], rect[3]], [BC[1], BC[1]], color='black', linewidth=2)
ax.plot([rect[0], rect[0]], [rect[2], rect[2]], [BC[0], BC[2]], color='black', linewidth=2)
ax.plot([rect[1], rect[1]], [rect[2], rect[2]], [BC[0], BC[3]], color='black', linewidth=2)
ax.plot([rect[0], rect[0]], [rect[3], rect[3]], [BC[1], BC[2]], color='black', linewidth=2)
ax.plot([rect[1], rect[1]], [rect[3], rect[3]], [BC[1], BC[3]], color='black', linewidth=2)
surf = ax.plot_surface(X, Y, Z_init, alpha = 0.7, cmap =cm.hot)
if iteration == 0:
plt.pause(5)
plt.pause(dt)
simulation_heat_eq(rect = [-5,1,-1,3], hs = [25,25], BC = [2,2,2,2], c = 5, frames = 1000)
What I realized was, however, that when clearing the graph, my heatmap on the surface would "recalibrate" in some sense, and although the temperature (z-value) was a lot lower in the later stages, the heatmap would still indicate low temperature because it's always comparing it to the surrounding temperatures rather than the initial temperature distribution in the room. See the images below
As you can see, it's still very red in the middle in the later stages of the simulation although it's temperature is very close to that of the borders of the room. I've tried looking for some commands that might help me resolve this problem in matplotlib, but haven't really found anything.
