Posted Saturday November 4 2023.
The idea of Green’s function is that if you know the unit impulse response of a linear system, then you can find solutions with arbitrary inputs by integrating the impulse response multiplied by the input. (I will also assume the system is translation-invariant, incl. time-invariant.)
The most well known Green’s function is probably the heat kernel, which is the Green’s function for the heat equation. It is the pdf of a normal distribution with mean zero and variance 2t. As t → 0+, the heat kernel “converges” pointwise to the Dirac delta, a singular distribution on the real line with unit mass concentrated at x = 0 and zero elsewhere - an “infinitely sharp spike” at zero. That is the unit impulse for real smooth functions.
Despite being a little nebulous and involving distributions for continuous real functions, the Green’s function for discrete spaces is relatively simple. For example, here is a numerical simulation on the circle where the unit impulse is just a one-hot vector :3 we can do diffusion via convolution and using a finite difference method.
The unit impulse response: the behavior of a small segment of the rod heated to unit temperature.
Now that we have the system response to the unit impulse, we can take some initial condition. How about a step function initial condition.
To get the response of the system to this initial data, we will look at the unit impulse response centered at each point of the initial data, multiplied by its height.
I debated multiple ways to display this, but how about this one: there is a white dot moving along at the top, which is the center of the impulse. At each point, it adds the unit impulse response centered at that point, multiplied by the height of the initial data.
Summing up (integrating) all the impulse responses, we obtain a surface telling us the heat of the rod at each point in space, at each point in time:
Indeed, this is a solution to the heat equation, and it satisfies our initial conditoin at t=0.
import numpy as np
import scipy.ndimage
# I will call this impulse_response
Nx = 50
Nt = 1000
x_max = 1
t_max = .01
dx = x_max/Nx
dt = t_max/Nt
impulse_response = np.zeros((Nt, Nx))
impulse_response[0, 0] = 1.0
eps = dt
K = dt/(dx**2)*np.array([1,-2,1])
for i in range(Nt-1):
impulse_response[i+1,:] = impulse_response[i,:] + scipy.ndimage.convolve(impulse_response[i,:], K, mode="wrap")
import matplotlib.pyplot as plt
plt.imshow(np.roll(impulse_response, Nx // 2), aspect="auto", extent=[0,1,t_max,0], cmap="hot")
plt.colorbar(label="temperature (farenheit)")
plt.xlabel("position along ring")
plt.ylabel("time")
plt.title("unit impulse response")
plt.show()
X = np.arange(0, 1, dx)
initial_data = 50 * np.heaviside(X - 0.2, 1.0) + 50 * np.heaviside(X - 0.7, 1.0) - 100 * np.heaviside(X - 0.8, 1.0)
plt.plot(X, initial_data, color="red")
plt.ylabel("temperature (farenheit)")
plt.xlabel("position along ring")
plt.title("initial data")
plt.show()
U = np.zeros((Nt, Nx))
images = []
from PIL import Image
for i in range(0,Nx):
IR_rolled = initial_data[i] * np.roll(impulse_response, shift=(i,0), axis=(1,))
U += IR_rolled
if i%(Nx//10)==0:
plt.imshow(U, vmin=0, vmax=100, aspect="auto", cmap="hot")
plt.show()
image = Image.fromarray(255/100*U).resize((256,256), Image.Resampling.NEAREST)
image.putpixel((int(i/Nx*256),0), 255)
images.append(image)
images[0].save("out.gif", save_all=True, append_images=images[1:], duration=100, loop=0)
from mpl_toolkits.mplot3d import Axes3D
x = np.linspace(0, U.shape[1]-1, U.shape[1])
y = np.linspace(0, U.shape[0]-1, U.shape[0])
X, Y = np.meshgrid(x, y)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surface = ax.plot_surface(X, Y, U)
ax.view_init(elev=30, azim=70)
ax.set_ylabel('Y-axis')
ax.set_zlabel('Z-axis')
# Remove all ticks
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])
# Remove tick labels
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_zticklabels([])
ax.set_xlabel('position')
ax.set_ylabel('time')
ax.set_zlabel('')
# Remove all ticks
ax.set_xticks([])
ax.set_yticks([])
ax.set_zticks([])
# Remove tick labels
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_zticklabels([])
plt.show()