p17: Helmholtz Equation in 2-D

u_{xx} + u_{yy} + k^{2}u = f, \qquad \mbox{on} \qquad [-1,1] \times [-1,1]

(1)

A minor modification of p16 to solve such problem for the particular choices as follows:

k = 9, \qquad f(x,y) = \exp (-10[(y-1)^2 +(x-1/2)^2])

(2)

#%config InlineBackend.figure_format='svg'
from chebPy import cheb
from numpy import meshgrid,sin,dot,eye,kron,zeros,reshape,exp,linspace
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.pyplot import figure,subplot,plot,title,axis,xlabel,ylabel,contour
from matplotlib import cm
from scipy.linalg import solve
from scipy.interpolate import RectBivariateSpline

N = 24; D,x = cheb(N); y = x;
xx,yy = meshgrid(x[1:N],y[1:N],indexing='ij')
xx = reshape(xx,(N-1)**2,order='F')
yy = reshape(yy,(N-1)**2,order='F')
f = exp(-10*((yy-1)**2 + (xx - 0.5)**2 ))
D2 = dot(D,D); D2 = D2[1:N,1:N]; I = eye(N-1)
k = 9
L = kron(I,D2) + kron(D2,I) + k**2*eye((N-1)**2)
# Solve Lu=f
u = solve(L,f)
# Convert 1-d vectors to 2-d
uu = zeros((N+1,N+1)); uu[1:N,1:N] = reshape(u,(N-1,N-1),order='F')
[xx,yy] = meshgrid(x,y,indexing='ij')
value = uu[N//2,N//2]

f = RectBivariateSpline(x,y,uu)
xxx = linspace(-1.0,1.0,60)
uuu = f(xxx,xxx)

fig = figure(figsize=(8,8))
ax = fig.add_subplot(111, projection='3d')
[X ,Y] = meshgrid(xxx,xxx,indexing='ij')
ax.plot_surface(X,Y,uuu,rstride=1,cstride=1,cmap=cm.jet,edgecolor='black')
title("$u(0,0)$="+str(value))
xlabel("x"); ylabel("y");

figure(figsize = (8,8))
contour(X,Y,uuu,levels=10);