{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# p17: Helmholtz Equation in 2-D\n",
"\n",
"$$ \n",
"u_{xx} + u_{yy} + k^{2}u = f, \\qquad \\mbox{on} \\qquad [-1,1] \\times [-1,1]\n",
"$$\n",
"\n",
"A minor modification of p16 to solve such problem for the particular choices as follows:\n",
"\n",
"$$\n",
"k = 9, \\qquad f(x,y) = \\exp (-10[(y-1)^2 +(x-1/2)^2])\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"%config InlineBackend.figure_format='svg'\n",
"from chebPy import cheb\n",
"from numpy import meshgrid,sin,dot,eye,kron,zeros,reshape,exp,linspace\n",
"from mpl_toolkits.mplot3d import Axes3D\n",
"from matplotlib.pyplot import figure,subplot,plot,title,axis,xlabel,ylabel,contour\n",
"from matplotlib import cm\n",
"from scipy.linalg import solve\n",
"from scipy.interpolate import RectBivariateSpline"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"N = 24; D,x = cheb(N); y = x;\n",
"xx,yy = meshgrid(x[1:N],y[1:N],indexing='ij')\n",
"xx = reshape(xx,(N-1)**2,order='F')\n",
"yy = reshape(yy,(N-1)**2,order='F')\n",
"f = exp(-10*((yy-1)**2 + (xx - 0.5)**2 ))\n",
"D2 = dot(D,D); D2 = D2[1:N,1:N]; I = eye(N-1)\n",
"k = 9\n",
"L = kron(I,D2) + kron(D2,I) + k**2*eye((N-1)**2)\n",
"# Solve Lu=f\n",
"u = solve(L,f)\n",
"# Convert 1-d vectors to 2-d\n",
"uu = zeros((N+1,N+1)); uu[1:N,1:N] = reshape(u,(N-1,N-1),order='F')\n",
"[xx,yy] = meshgrid(x,y,indexing='ij')\n",
"value = uu[N//2,N//2]\n",
"\n",
"f = RectBivariateSpline(x,y,uu)\n",
"xxx = linspace(-1.0,1.0,60)\n",
"uuu = f(xxx,xxx)\n",
"\n",
"fig = figure(figsize=(8,8))\n",
"ax = fig.add_subplot(111, projection='3d')\n",
"[X ,Y] = meshgrid(xxx,xxx,indexing='ij')\n",
"ax.plot_surface(X,Y,uuu,rstride=1,cstride=1,cmap=cm.jet,edgecolor='black')\n",
"title(\"$u(0,0)$=\"+str(value))\n",
"xlabel(\"x\"); ylabel(\"y\");\n",
"\n",
"figure(figsize = (8,8))\n",
"contour(X,Y,uuu,levels=10);"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.8"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 4
}