# @CX-REMOVE
!pip install Pillow
Requirement already satisfied: Pillow in /opt/conda/lib/python3.11/site-packages (10.1.0)
Convolution Algebras and Discrete Geometry - a Study¶
We start our exploration of numerical geometry with a function on a plane.
We represent this function by an array: F[x,y] where x,y are points on an equidistant grid with step size h.
# Example Function
#@CX-TOGGLE
import numpy as np
import matplotlib.pyplot as plt
def f(x, y):
return (np.sin(x) + np.cos(y) + 0.3 + 3*x - 0.2*y) / (1 + x**2 + y**2 + 20*(0.2*x*y + 0.03*x + 0.02*y)**2)
# return np.arctan(1*x - 2*y + 0.1 * (x**4 + y**4) )
# Generate a grid of (x, y) coordinates
N = 100
W = 6
x = np.linspace(-W, W, 2*N)
y = np.linspace(-W, W, 2*N)
X, Y = np.meshgrid(x, y) # This creates a 2D grid of x, y values
h = x[1]-x[0]
# Compute f(x,y) for each coordinate in the grid
F = f(X, Y)
# Plot the heatmap
def PLOT(F):
plt.figure(figsize=(5,5))
plt.imshow(F, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r')
plt.colorbar(shrink=0.8)
plt.xlabel('x')
plt.ylabel('y')
PLOT(F)
plt.title('Heatmap of F');
Directional Derivatives¶
We compute the "balanced" directional derivatives by convolution on the array F:
DxF[x,y] = (F[x+1] - F[x-1])/2/h = [1, 0, -1]/2/h @ F
Here @ denotes the discrete convolution operation.
The formula is obtained as the arithmetic mean of the forward differential (F[x+1]-F[x])/h and the backward differential (F[x]-F[x-1])/h.
# Kernels and Derivatives
#@CX-TOGGLE
import numpy as np
from scipy.signal import convolve2d
class Kernel(np.ndarray):
def __new__(cls, input_array):
obj = np.asarray(input_array).view(cls)
return obj
def __matmul__(self, other):
"""Overloads @ operator for convolution."""
if isinstance(other, Kernel):
return Kernel(convolve2d(self, other, mode='full', boundary='fill'))
elif isinstance(other, np.ndarray):
return convolve2d(other, self, mode='same', boundary='symm')
else:
raise ValueError("Unsupported type for convolution.")
def __add__(self, other):
# Get the shape of the larger kernel
shape = (max(self.shape[0], other.shape[0]), max(self.shape[1], other.shape[1]))
result = np.zeros(shape)
# x: [ ----------------------------- ]
# [ pad_0 ][ self ][ pad_1 ]
padx = (shape[0] - self.shape[0]) // 2
pady = (shape[1] - self.shape[1]) // 2
result[ padx : padx + self.shape[0], pady : pady + self.shape[1] ] += self
padx = (shape[0] - other.shape[0]) // 2
pady = (shape[1] - other.shape[1]) // 2
result[ padx : padx + other.shape[0], pady : pady + other.shape[1] ] += other
return Kernel(result)
from scipy.signal import convolve2d
import numpy as np
import matplotlib.pyplot as plt
# Compute the partial derivatives using convolution
Dx = Kernel(np.array([[1, 0, -1]]))/2/h
Dy = Dx.T
DxF = Dx @ F
DyF = Dy @ F
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
# Plot F
im1 = axes[0].imshow(F, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r')
axes[0].set_title('F')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
fig.colorbar(im1, ax=axes[0])
# Plot DxF
im2 = axes[1].imshow(DxF, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r')
axes[1].set_title('Dx F')
axes[1].set_xlabel('x')
axes[1].set_ylabel('y')
fig.colorbar(im2, ax=axes[1])
# Plot DyF
im3 = axes[2].imshow(DyF, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r')
axes[2].set_title('Dy F')
axes[2].set_xlabel('x')
axes[2].set_ylabel('y')
fig.colorbar(im3, ax=axes[2])
plt.tight_layout()
plt.show()
Discrete Laplace Operator¶
The Laplace Operator is defined as , can be expressed via convolution:
LF = Dx @ Dx @ F + Dy @ Dy @ F = (Dx @ Dx + Dy @ Dy) @ F
The convolution kernel L = (Dx @ Dx + Dy @ Dy) looks as follows:
Dxx = Dx @ Dx
Dyy = Dy @ Dy
L0 = (Dxx + Dyy)
print(L0*h*h*4)
[[ 0. 0. 1. 0. 0.] [ 0. 0. 0. 0. 0.] [ 1. 0. -4. 0. 1.] [ 0. 0. 0. 0. 0.] [ 0. 0. 1. 0. 0.]]
We can "compress" this into the standard Laplace Kernel:
L = Kernel(np.array([[0,1,0],[1,-4,1],[0,1,0]]))/h/h/4
print(L*h*h*4)
[[ 0. 1. 0.] [ 1. -4. 1.] [ 0. 1. 0.]]
Note that, convolution with this kernel, compares each F[x,y] to the average over the adjacent values.
# Plot Laplace of F
# @CX_TOGGLE
def PLOT2(F, G, imshow_args={}):
# Create subplots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# Plot the original matrix F
ax1.imshow(F, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r', **imshow_args)
ax1.set_xlabel('x')
ax1.set_ylabel('y')
# Plot the result of applying Laplacian to F
ax2.imshow(G, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r', **imshow_args)
ax2.set_xlabel('x')
ax2.set_ylabel('y')
# Adjust spacing between subplots
plt.tight_layout()
return (ax1, ax2)
ax1, ax2 = PLOT2(F, L@F, {"vmin":-2, "vmax":3})
ax1.set_title('Original F')
ax2.set_title('Laplacian of F');
Heat Flow¶
The heat equation is given by $\partial F / \partial t = \Delta F$, heuristically, we changethe value of F along in the direction of the Laplacian at that point.
It's straight forward to simulate the heat flow numerically by subsituting:
F = F + eps * L @ F
in an interative way.
When trying this out, naively, we run into numeric instabilities with checker patterns "blowing up". In order to get a numerically stable version, we need to replace the Laplacian with another version, that takes into account diagonal entries. The heuristic of averaging the neighbours remains the same:
# Stable 2D Laplace Kernel
LS = Kernel(np.array([[0.25,.5,.25],[.5,-3,.5],[.25,.5,.25]]))/h/h/4
print(LS*h*h*4)
[[ 0.25 0.5 0.25] [ 0.5 -3. 0.5 ] [ 0.25 0.5 0.25]]
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.signal import convolve2d
G = F.copy()
# Create a figure and axis for the animation
fig, ax = plt.subplots(figsize=(5, 5))
cax = ax.imshow(G, extent=[-W, W, -W, W], origin='lower', cmap='RdYlBu_r', vmin=-2, vmax=2)
plt.colorbar(cax, shrink=0.8)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Heat Equation Simulation')
# Define a function to update the animation
def update(frame):
global G
for _ in range(30):
G = G + 0.001 * convolve2d(G, LS, mode="same", boundary="symm")
cax.set_array(G)
return cax,
# Create the animation
animation = FuncAnimation(fig, update, frames=100, repeat=True, blit=True)
# Display the animation in the Jupyter notebook
from IPython.display import HTML
html = animation.to_jshtml(fps=20)
HTML(html)