Skip to content

Instantly share code, notes, and snippets.

Created September 9, 2016 01:50
Show Gist options
  • Save anonymous/e30eb765c2e2b5bf6675f3f2b0714ddb to your computer and use it in GitHub Desktop.
Save anonymous/e30eb765c2e2b5bf6675f3f2b0714ddb to your computer and use it in GitHub Desktop.
Compute and plot Pearcey integral
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate
# Dimension of image in pixels
N = 256
# Number of samples to use for integration
M = 32
# Magnitude of contour displacement.
# Determined by trial and error with a bunch of contour plots
# of the integrand.
rate = 0.01
# Tabulate Pe(α,β) = ∫(x = -∞,∞) exp(i(x⁴+αx²+β))dx
table = np.zeros((N, N), dtype=np.float64)
xmin, xmax = (-10.0, 10.0)
ymin, ymax = (4.0, -16.0)
alphas = np.linspace(ymin, ymax, N)
betas = np.linspace(xmin, xmax, N)
# Integration range.
# The integrand drops to zero so fast this range
# is fine for now.
x = np.linspace(-4.0, 4.0, M)
for i, alpha in enumerate(alphas):
for j, beta in enumerate(betas):
# f(z) = z⁴+αz²+β
# g(z) = exp(if(z))
# Instead of integrating along x-axis we're
# going to integrate along a contour displaced
# vertically from the x-axis.
# A good choice of displacement is the gradient
# d/(Im f(x+iy))/dy.
# That way, we're displacing in a direction that makes
# |exp(if(x+iy))| smaller.
y = rate*(4*x**3+2*alpha*x+beta)
z = x+1j*y
f = z**4+alpha*z**2+beta*z
g = np.exp(1j*f)
# ∫f(z)dz = ∫f(z)dz/dx dz
dz = 1.0+1j*rate*(12*x**2+2*alpha)
I = integrate.simps(g*dz, x)
table[i, j] = np.abs(I)
plt.imshow(table, cmap='gray', extent=[xmin, xmax, ymax, ymin])
plt.title('|Pe(x,y)|')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar()
plt.show()
# Refs.
# [1] Berry and Klein, Colored diffraction catastrophes
# http://www.pnas.org/content/93/6/2614.full.pdf
# [2] https://en.wikipedia.org/wiki/Pearcey_integral
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment