Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
2D Hamilton-Jacobi Equation
import numpy as np
import matplotlib.pyplot as plt
import math
import numba
def buildSphere(cx, cy, r):
return lambda x, y: math.sqrt((x - cx)**2 + (y - cy)**2) - r
@numba.jit
def Convect(X, Y, Z, dt, speed):
hx = X[0, 1] - X[0, 0]
hy = Y[1, 0] - Y[0, 0]
row, column = Z.shape
u = np.empty(Z.shape)
for r in range(row):
for c in range(column):
z = Z[r, c]
dx = (Z[r, (c + 1) % column] - Z[r, (c - 1) % column]) / (2 * hx)
dy = (Z[(r + 1) % row, c] - Z[(r - 1) % row, c]) / (2 * hy)
g = math.sqrt(dx * dx + dy * dy)
u[r, c] = z - dt * speed * g
return u
def step(Z, n, dt, speed):
Z1 = np.array(Z)
for i in range(n):
Z1 = Convect(X, Y, Z1, dt, speed)
return Z1
def plot(x, y, z, z1):
plt.grid(True)
plt.axis('equal')
plt.contour(x, y, z, 0)
plt.contour(x, y, z1, 0)
plt.show()
if __name__ == '__main__':
x = np.arange(-2, 4, 0.01)
y = np.arange(-2, 4, 0.01)
X, Y = np.meshgrid(x, y)
s0 = buildSphere(0.003, 0.003, 1)
s1 = buildSphere(2 - 0.01 + 0.003, 0.003, 1)
f = lambda x, y: min(s0(x, y), s1(x, y))
F = np.vectorize(f)
Z = F(X, Y)
Z1 = step(Z, 700, 0.0001, -2)
plot(X, Y, Z, Z1)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.