{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# p16: Poisson equation in 2-D"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We solve the following Poisson problem\n",
"\n",
"$$\n",
"u_{xx} + u_{yy} = 10\\sin(8x(y-1)), \\qquad -1 < x,y < 1, \\qquad u=0 \\quad \\mbox{on boundary}\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%config InlineBackend.figure_format='svg'\n",
"from chebPy import cheb\n",
"from numpy import meshgrid,sin,dot,eye,kron,zeros,reshape,linspace\n",
"from matplotlib.pyplot import figure,subplot,plot,title,axis,xlabel,ylabel,spy\n",
"from matplotlib import cm\n",
"from scipy.linalg import solve\n",
"from scipy.interpolate import RegularGridInterpolator"
]
},
{
"cell_type": "code",
"execution_count": 2,
"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 = 10*sin(8*xx*(yy-1))\n",
"D2 = dot(D,D); D2 = D2[1:N,1:N]; I = eye(N-1)\n",
"L = kron(I,D2) + kron(D2,I)\n",
"# Plot sparsity pattern\n",
"figure(figsize=(6,6)), spy(L)\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[3*N//4,3*N//4]\n",
"\n",
"# Interpolate to finer mesh just for visualization\n",
"f = RegularGridInterpolator((x,y),uu,method='cubic')\n",
"xxx = linspace(-1.0,1.0,50); X,Y = meshgrid(xxx,xxx,indexing='ij')\n",
"uuu = f((X,Y))\n",
"fig = figure(figsize=(6,6))\n",
"ax = fig.add_subplot(111, projection='3d')\n",
"ax.plot_surface(X,Y,uuu,rstride=1,cstride=1,cmap=cm.jet,edgecolor='black')\n",
"title(\"$u(2^{-1/2},2^{-1/2})$=\"+str(value))\n",
"xlabel(\"x\"), ylabel(\"y\");"
]
}
],
"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"
}
},
"nbformat": 4,
"nbformat_minor": 4
}