This notebook was prepared as part of an answer to this physics stackexchange question.
It it demonstrating an example application of the method described in the paper: Reviving the Method of Particular Solutions by Timo Betcke and Lloyd N. Trefethen from 2005. [doi] [pdf]
The idea behind the method is to choose some parametric basis of solutions to some part of the broader Helmholtz equation, and then to find the eigenvalues for which there exist a linear combination of this basis that satisfies the boundary conditions.
Schematically, if we have some parametric basis: $ \phi_k(x;\lambda)$ we see a solution of the form $ \sum_k c_k \phi_k(x;\lambda) $ and we do this by choose a set of points on the boundary, for which we will ensure that our solution vanishes, and a set of points on the interior for which the solution doesn't vanish.
The details are worked out in the paper, but this means you choose a set of points $x_i$ on both the boundary and interior, then we define $\sigma(\lambda)$ in the following way:
# just some standard imports
import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
from scipy.special import jv as jn
import matplotlib as mpl
mpl.cm.register_cmap(name='cubehelix2', data=mpl._cm.cubehelix(gamma=1.0, s=0.3, r=-0.5, h=1.5))
%matplotlib inline
In particular, in this example we will be interested in finding the eigenvalues and eigenfunctions for an L shaped domain, let's choose our points, both on the interior and boundary randomly.
N = 36
NP = 2*N
k = (np.arange(N)+1.0)
alpha = 2/3.0
t1 = 1.5 * np.pi * (np.arange(NP)+0.5)/NP
r1 = 1./np.maximum( np.abs(np.sin(t1)), np.abs(np.cos(t1)))
t2 = 1.5 * np.pi*sp.rand(NP)
r2 = sp.rand(NP)/np.maximum( np.abs(np.sin(t2)), np.abs(np.cos(t2)) )
t = np.r_[t1,t2]
r = np.r_[r1,r2]
Just so we can see what's going on, these are the points we've sampled:
xs = r*np.cos(t)
ys = r*np.sin(t)
plt.scatter(xs,ys);
plt.axis('equal');
plt.axis('off');
Now we will define $\sigma(\lambda)$ as we described above.
For our basis of representations, we will use the functions
$$ \phi_k(x; \lambda) = j_{\alpha k} ( \sqrt{\lambda} r ) \sin(\alpha \theta k ) $$since these are the solution to the helmholtz equation on a semiinfinite wedge. If we use these functions at the interior angle of our $L$ shape, this will guarentee we satisfy the boundary conditions along the interior walls of the $L$, which will help us out quite a bit with the accuracy. Here, since our interior angle is $3\pi/2$, we take $\alpha = 2/3$
def mat(lam):
return np.sin(alpha*t[:,None]*k[None,:])*jn(alpha*k[None,:], np.sqrt(lam)*r[:,None])
def sigma(lam):
A = mat(lam)
Q,R = sp.linalg.qr(A, mode='economic')
s = sp.linalg.svdvals(Q[:NP]).min()
return s
To see how this function behaves, let's look at it as a function of $\lambda$.
lamvec = np.r_[0.2:40:0.2]
plt.plot(lamvec, map(sigma, lamvec));
plt.xlabel(r"$\lambda$")
plt.ylabel(r"$\sigma(\lambda)$");
We can clearly see those points where it is vanishing, those correspond to our eigenvalues. Let's use a scalar minimizer to find the first 6 values.
from scipy.optimize import minimize_scalar
lams = [ minimize_scalar(sigma,(z-0.1,z+0.1), tol=1e-16).x for z in [9.6, 15.2, 19.7, 29.5, 31.9, 41.5] ]
print lams
[9.639723844634462, 15.197251925760987, 19.739208802179888, 29.521481117131952, 31.912635948136323, 41.474509890573152]
Great, these agree very well with the results in the paper. Now that we know what the eigenvalues are to good numerical precision, let's determine the eigenfunctions. Following the method, these are determined by
$$ R(\lambda) c = \hat v $$where $R$ is the R from the QR decomposition of $A(\lambda)$ as above, $c$ is the vector of coefficients defining our solution, and $\hat v$ is the right singular vector associated with the smallest eigenvalue of the boundary part of $Q$.
We can in fact plot the lowest 6 eigenfunctions and give their associated eigenvalues.
x, y = np.ogrid[-1:1:200j, -1:1:200j]
rr = np.sqrt( x**2 + y**2 )
tt = (sp.arctan2(y,x)+2*np.pi)%(2*np.pi)
fig,axs = plt.subplots(2,3, figsize=(10,5))
for pk,ax in enumerate(axs.flat):
Q,R = sp.linalg.qr(mat(lams[pk]), mode='economic')
u,s,vt = sp.linalg.svd( Q[:NP,:], full_matrices=False )
c = sp.linalg.solve(R, vt[-1])
ans = (np.sin(alpha*tt[:,:,None]*k[None,None,:])*jn(alpha*k[None,None,:], np.sqrt(lams[pk])*rr[:,:,None])).dot(
c)
ans[100:,:100] = np.nan
#ax.imshow( ans[:,::-1],cmap='cubehelix2', extent=[-1,1,-1,1]);
ax.contourf( ans[:,::-1], 10, cmap='cubehelix2')
ax.set_title(r"$\lambda_{} = {}$".format(pk, lams[pk]))
ax.axis('off');
ax.axis('equal');
plt.tight_layout()