{
"cells": [
{
"cell_type": "markdown",
"id": "d22ffe10",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"# Demo - Mixed bases for the Helmholtz problem\n",
"**Mikael Mortensen** (email: `mikaem@math.uio.no`), Department of Mathematics, University of Oslo.\n",
"\n",
"Date: **March 22, 2021**\n",
"\n",
"**Summary.** This demo shows how to solve the Helmholtz equation using different\n",
"bases for test and trial spaces. The use of different bases leads for\n",
"some optimal combinations to highly sparse and well-conditioned\n",
"coefficient matrices."
]
},
{
"cell_type": "markdown",
"id": "e66d8420",
"metadata": {
"editable": true
},
"source": [
"## The Helmholtz problem\n",
"\n",
"We will consider Helmholtz equation with homogeneous Dirichlet boundary conditions"
]
},
{
"cell_type": "markdown",
"id": "fc92d314",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" \\alpha u - u^{''} = f \\quad \\text{in} \\, {I}=(-1, 1), \\quad u(\\pm 1) = 0,\n",
"\\label{_auto1} \\tag{1}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "df92a93f",
"metadata": {
"editable": true
},
"source": [
"where $\\alpha \\in \\mathbb{R^+}$. The relevant function space for the Dirichlet problem is"
]
},
{
"cell_type": "markdown",
"id": "c0ee2909",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
" S_N=\\text{span}\\{T_k\\}_{k=0}^{N-1}, \\quad V_{N} = \\{v \\in {S}_N\\,|\\, v(\\pm 1) = 0\\},\n",
"\\label{_auto2} \\tag{2}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "cfb1e1a0",
"metadata": {
"editable": true
},
"source": [
"and the Chebyshev-Galerkin (CG) method is to find $u_N \\in V_N$ such that"
]
},
{
"cell_type": "markdown",
"id": "c0137bd8",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\label{eq:dirGalerkin} \\tag{3}\n",
" \\alpha (u_N, v)_{\\omega^{\\sigma}} -(u^{''}_N, v)_{\\omega^{\\sigma}} = (f, v)_{\\omega^{\\sigma}}, \\forall \\, v \\in V_N,\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "118ec8b5",
"metadata": {
"editable": true
},
"source": [
"where $(u,v)_{\\omega^{\\sigma}}=\\int_{{I}}uv\\omega^{\\sigma} dx$ is the scalar product in the weighted space $L^2_{\\omega^{\\sigma}}({I})$.\n",
"\n",
"Shenfun has implemented three different Chebyshev Dirichlet\n",
"basis functions"
]
},
{
"cell_type": "markdown",
"id": "8a6ae7d4",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\label{eq:shen} \\tag{4}\n",
"\\phi_k = T_k-T_{k+2}, \\quad k=0,1, \\ldots, N-3,\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "56b8b2b8",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\label{eq:heinrichs} \\tag{5}\n",
"\\varphi_k = (1-x^2)T_k, \\quad k=0,1, \\ldots, N-3,\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "4bb15362",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\label{eq:dirichletU} \\tag{6}\n",
"\\psi_k = U_k-\\frac{k+1}{k+3}U_{k+2}, \\quad k=0,1, \\ldots, N-3.\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "f0bf61eb",
"metadata": {
"editable": true
},
"source": [
"These three bases are all linearly dependent and they are all bases\n",
"for $V_N$."
]
},
{
"cell_type": "markdown",
"id": "0df870c4",
"metadata": {
"editable": true
},
"source": [
"## Implementation\n",
"\n",
"We can get all three function spaces as"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "f1f9ca5d",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:39.361760Z",
"iopub.status.busy": "2024-02-19T13:48:39.361483Z",
"iopub.status.idle": "2024-02-19T13:48:40.276844Z",
"shell.execute_reply": "2024-02-19T13:48:40.276127Z"
}
},
"outputs": [],
"source": [
"from shenfun import *\n",
"N = 40\n",
"V0 = FunctionSpace(N, 'C', basis='ShenDirichlet')\n",
"V1 = FunctionSpace(N, 'C', basis='Heinrichs')\n",
"V2 = FunctionSpace(N, 'U', basis='CompactDirichlet')"
]
},
{
"cell_type": "markdown",
"id": "e31306bc",
"metadata": {
"editable": true
},
"source": [
"where $V0 = \\text{span}\\{\\phi_k\\}_{k=0}^{N-3}$,\n",
"$V1 = \\text{span}\\{\\varphi_k\\}_{k=0}^{N-3}$ and\n",
"$V2 = \\text{span}\\{\\psi_k\\}_{k=0}^{N-3}$. Now, to solve the Helmholtz problem we simply need to choose\n",
"test and trial bases. Shen's original method is using\n",
"`V0` for both. To assemble the stiffness and mass matrices\n",
"for this choice do"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "c0a4b817",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:40.280230Z",
"iopub.status.busy": "2024-02-19T13:48:40.279834Z",
"iopub.status.idle": "2024-02-19T13:48:40.285055Z",
"shell.execute_reply": "2024-02-19T13:48:40.284589Z"
}
},
"outputs": [],
"source": [
"u = TrialFunction(V0)\n",
"v = TestFunction(V0)\n",
"A = inner(v, div(grad(u)))\n",
"B = inner(v, u)"
]
},
{
"cell_type": "markdown",
"id": "ebc3bfe4",
"metadata": {
"editable": true
},
"source": [
"A manufactured solution can be chosen using [Sympy](https://www.sympy.org)\n",
"We choose"
]
},
{
"cell_type": "markdown",
"id": "116f1647",
"metadata": {
"editable": true
},
"source": [
"\n",
"\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(x) = \\sin \\left( 2 \\pi \\cos \\left( 2 \\pi x \\right) \\right)\n",
"\\label{_auto3} \\tag{7}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"id": "3ac72ba9",
"metadata": {
"editable": true
},
"source": [
"implemented as"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "b621fdd7",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:40.288102Z",
"iopub.status.busy": "2024-02-19T13:48:40.287531Z",
"iopub.status.idle": "2024-02-19T13:48:40.294925Z",
"shell.execute_reply": "2024-02-19T13:48:40.294417Z"
}
},
"outputs": [],
"source": [
"import sympy as sp\n",
"x = sp.Symbol('x', real=True)\n",
"ue = sp.sin(2*sp.pi*sp.cos(2*sp.pi*x))"
]
},
{
"cell_type": "markdown",
"id": "8ed01c45",
"metadata": {
"editable": true
},
"source": [
"The right hand side $f$ of Helmholtz equation is"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "b7648fc4",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:40.297685Z",
"iopub.status.busy": "2024-02-19T13:48:40.297462Z",
"iopub.status.idle": "2024-02-19T13:48:40.923369Z",
"shell.execute_reply": "2024-02-19T13:48:40.922562Z"
}
},
"outputs": [],
"source": [
"alpha = 1\n",
"f = sp.simplify(alpha*ue-ue.diff(x, 2))"
]
},
{
"cell_type": "markdown",
"id": "5c02f31e",
"metadata": {
"editable": true
},
"source": [
"To solve the problem we can do"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "341119de",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:40.926180Z",
"iopub.status.busy": "2024-02-19T13:48:40.925963Z",
"iopub.status.idle": "2024-02-19T13:48:40.941262Z",
"shell.execute_reply": "2024-02-19T13:48:40.940562Z"
}
},
"outputs": [],
"source": [
"fj = Array(V0, buffer=f) # Get f on quadrature mesh\n",
"f_hat = inner(v, fj) # Compute right hand side\n",
"M = alpha*B - A # Get coefficient matrix\n",
"u_hat = Function(V0) # Container for the solution\n",
"sol = la.Solver(M) # Solver\n",
"u_hat = sol(f_hat, u_hat) # Solve"
]
},
{
"cell_type": "markdown",
"id": "ff687ebc",
"metadata": {
"editable": true
},
"source": [
"Compare with exact solution."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "c29b3b52",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:40.944831Z",
"iopub.status.busy": "2024-02-19T13:48:40.944580Z",
"iopub.status.idle": "2024-02-19T13:48:40.983675Z",
"shell.execute_reply": "2024-02-19T13:48:40.982789Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Error = 0.41890381687693135\n"
]
}
],
"source": [
"uj = Array(V0, buffer=ue)\n",
"error = inner(1, (u_hat.backward()-uj)**2)\n",
"print('Error =', error)"
]
},
{
"cell_type": "markdown",
"id": "d3104bb4",
"metadata": {
"editable": true
},
"source": [
"Now that was the solution for test and trial bases from the same\n",
"basis ([4](#eq:shen)). Let us create a function that takes any\n",
"test and any trial basis, any manufactured solution and any $\\alpha$\n",
"in the Helmholtz equation. We let the function return either\n",
"the L2-error norm, the condition number of the Helmholtz\n",
"coefficient matrix, or the matrix itself."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "e7fc2aff",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:40.986318Z",
"iopub.status.busy": "2024-02-19T13:48:40.986123Z",
"iopub.status.idle": "2024-02-19T13:48:41.038794Z",
"shell.execute_reply": "2024-02-19T13:48:41.037921Z"
}
},
"outputs": [],
"source": [
"def main(N, test, trial, alpha=1, method=0, ue=sp.sin(2*sp.pi*sp.cos(2*sp.pi*x))):\n",
" \"\"\"Solve Helmholtz problem and return L2-error, condition number or matrix\n",
"\n",
" Parameters\n",
" ----------\n",
" N : int\n",
" Number of quadrature points\n",
" test, trial : int\n",
" Test and trial functions.\n",
" 0 = :math:`T_k-T_{k+2}`\n",
" 1 = :math:`(1-x^2)T_k`\n",
" 2 = :math:`U_k-\\frac{k+1}{k+3}U_{k+2}`\n",
" alpha : Helmholtz parameter\n",
" method : int\n",
" 0 = Return L2-error norm\n",
" 1 = Return condition number of matrix\n",
" 2 = Return Helmholtz matrix\n",
" ue : Sympy function, optional\n",
" The manufactured solution with homogeneous boundary conditions.\n",
"\n",
" Note\n",
" ----\n",
" Inhomogeneous boundary conditions require a small rewrite, but is\n",
" not difficult.\n",
"\n",
" \"\"\"\n",
" bases = {0: ('C', 'ShenDirichlet'), 1: ('C', 'Heinrichs'), 2: ('U', 'CompactDirichlet')}\n",
" test = FunctionSpace(N, bases[test][0], basis=bases[test][1])\n",
" trial= FunctionSpace(N, bases[trial][0], basis=bases[trial][1])\n",
" # Check that boundary conditions are homogeneous\n",
" assert abs(ue.subs(x, -1)) < 1e-8 and abs(ue.subs(x, 1)) < 1e-8\n",
" u = TrialFunction(trial)\n",
" v = TestFunction(test)\n",
" f = sp.simplify(alpha*ue-ue.diff(x, 2))\n",
" fj = Array(test, buffer=f) # Get f on quadrature mesh\n",
" f_hat = inner(v, fj) # Compute right hand side\n",
" B = inner(v, u)\n",
" A = inner(v, div(grad(u)))\n",
" M = alpha*B-A\n",
" if method == 1:\n",
" return np.linalg.cond(M.diags('csr').toarray())\n",
" if method == 2:\n",
" return M\n",
"\n",
" u_hat = Function(trial)\n",
" sol = la.Solver(M)\n",
" u_hat = sol(f_hat, u_hat)\n",
" uj = Array(trial, buffer=ue)\n",
" error = np.sqrt(inner(1, (u_hat.backward()-uj)**2))\n",
" return error"
]
},
{
"cell_type": "markdown",
"id": "311641c8",
"metadata": {
"editable": true
},
"source": [
"Let us first try basis ([4](#eq:shen)) as test function and\n",
"([5](#eq:heinrichs)) as trial function. Use otherwise\n",
"default parameters."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "6a61d76e",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:41.042436Z",
"iopub.status.busy": "2024-02-19T13:48:41.042140Z",
"iopub.status.idle": "2024-02-19T13:48:41.474999Z",
"shell.execute_reply": "2024-02-19T13:48:41.473797Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"3.8866412036007366e-07\n"
]
}
],
"source": [
"error = main(100, 0, 1)\n",
"print(error)"
]
},
{
"cell_type": "markdown",
"id": "19538c44",
"metadata": {
"editable": true
},
"source": [
"So the error is small in deed. Perhaps more interesting, let's\n",
"look at the sparsity pattern of the coefficient matrix"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "21ff1456",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:41.479228Z",
"iopub.status.busy": "2024-02-19T13:48:41.478935Z",
"iopub.status.idle": "2024-02-19T13:48:42.514067Z",
"shell.execute_reply": "2024-02-19T13:48:42.513423Z"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"application/vnd.plotly.v1+json": {
"config": {
"plotlyServerURL": "https://plot.ly"
},
"data": [
{
"hovertemplate": "x: %{x}
y: %{y}
color: %{z[0]}",
"name": "0",
"source": "data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAGIAAABiCAAAAACOzlnSAAABY0lEQVR4XrXRORIEMQwCQMn//7M3mGNtjw8dQEjSRaFStWoVYkrVqlXHGpmidONewTTK9QXTKEI3itCNIsI2ioiQjYugGjfBNB6CaLwEz/gTNKMhWEZLkIyO4Bg9QTEGgmGMBMH4EHjjS8CNCYE2ZgTYmBJYY05AjQWBNFYE0FgSOGNNwIwNgTJ2BMjYEhhjT0CMA4EwTgTAOBJ540ykDQORNSxE0jAROcNGpAwjkTGsRMIwE3HDToQNBxE1PETQcBExw0eEDCcRMbxEwHATfsNPuI0A4TUihNMIET4jRriMIOExooTDCBN2I06YjQRhNTKE0UgRNiNHmIwkYTGyhMFIE2cjTxwNAHEyEMTBgBB7A0NsDRCxM1DExoARawNHLA0gsTKQxMKAEnMDS0wNMDEz0MTEgBNfA098DAIxGgxiMChEb3CIziARrcEiGoNG/A0e8RpE4jGYxG1QicvgElK16g+xo4HRIs8XKAAAAABJRU5ErkJggg==",
"type": "image",
"xaxis": "x",
"yaxis": "y"
}
],
"layout": {
"margin": {
"t": 60
},
"template": {
"data": {
"bar": [
{
"error_x": {
"color": "#2a3f5f"
},
"error_y": {
"color": "#2a3f5f"
},
"marker": {
"line": {
"color": "#E5ECF6",
"width": 0.5
},
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "bar"
}
],
"barpolar": [
{
"marker": {
"line": {
"color": "#E5ECF6",
"width": 0.5
},
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "barpolar"
}
],
"carpet": [
{
"aaxis": {
"endlinecolor": "#2a3f5f",
"gridcolor": "white",
"linecolor": "white",
"minorgridcolor": "white",
"startlinecolor": "#2a3f5f"
},
"baxis": {
"endlinecolor": "#2a3f5f",
"gridcolor": "white",
"linecolor": "white",
"minorgridcolor": "white",
"startlinecolor": "#2a3f5f"
},
"type": "carpet"
}
],
"choropleth": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "choropleth"
}
],
"contour": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "contour"
}
],
"contourcarpet": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "contourcarpet"
}
],
"heatmap": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "heatmap"
}
],
"heatmapgl": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "heatmapgl"
}
],
"histogram": [
{
"marker": {
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "histogram"
}
],
"histogram2d": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "histogram2d"
}
],
"histogram2dcontour": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "histogram2dcontour"
}
],
"mesh3d": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "mesh3d"
}
],
"parcoords": [
{
"line": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "parcoords"
}
],
"pie": [
{
"automargin": true,
"type": "pie"
}
],
"scatter": [
{
"fillpattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
},
"type": "scatter"
}
],
"scatter3d": [
{
"line": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatter3d"
}
],
"scattercarpet": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattercarpet"
}
],
"scattergeo": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattergeo"
}
],
"scattergl": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattergl"
}
],
"scattermapbox": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattermapbox"
}
],
"scatterpolar": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterpolar"
}
],
"scatterpolargl": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterpolargl"
}
],
"scatterternary": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterternary"
}
],
"surface": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "surface"
}
],
"table": [
{
"cells": {
"fill": {
"color": "#EBF0F8"
},
"line": {
"color": "white"
}
},
"header": {
"fill": {
"color": "#C8D4E3"
},
"line": {
"color": "white"
}
},
"type": "table"
}
]
},
"layout": {
"annotationdefaults": {
"arrowcolor": "#2a3f5f",
"arrowhead": 0,
"arrowwidth": 1
},
"autotypenumbers": "strict",
"coloraxis": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"colorscale": {
"diverging": [
[
0,
"#8e0152"
],
[
0.1,
"#c51b7d"
],
[
0.2,
"#de77ae"
],
[
0.3,
"#f1b6da"
],
[
0.4,
"#fde0ef"
],
[
0.5,
"#f7f7f7"
],
[
0.6,
"#e6f5d0"
],
[
0.7,
"#b8e186"
],
[
0.8,
"#7fbc41"
],
[
0.9,
"#4d9221"
],
[
1,
"#276419"
]
],
"sequential": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"sequentialminus": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
]
},
"colorway": [
"#636efa",
"#EF553B",
"#00cc96",
"#ab63fa",
"#FFA15A",
"#19d3f3",
"#FF6692",
"#B6E880",
"#FF97FF",
"#FECB52"
],
"font": {
"color": "#2a3f5f"
},
"geo": {
"bgcolor": "white",
"lakecolor": "white",
"landcolor": "#E5ECF6",
"showlakes": true,
"showland": true,
"subunitcolor": "white"
},
"hoverlabel": {
"align": "left"
},
"hovermode": "closest",
"mapbox": {
"style": "light"
},
"paper_bgcolor": "white",
"plot_bgcolor": "#E5ECF6",
"polar": {
"angularaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"bgcolor": "#E5ECF6",
"radialaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
}
},
"scene": {
"xaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
},
"yaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
},
"zaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
}
},
"shapedefaults": {
"line": {
"color": "#2a3f5f"
}
},
"ternary": {
"aaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"baxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"bgcolor": "#E5ECF6",
"caxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
}
},
"title": {
"x": 0.05
},
"xaxis": {
"automargin": true,
"gridcolor": "white",
"linecolor": "white",
"ticks": "",
"title": {
"standoff": 15
},
"zerolinecolor": "white",
"zerolinewidth": 2
},
"yaxis": {
"automargin": true,
"gridcolor": "white",
"linecolor": "white",
"ticks": "",
"title": {
"standoff": 15
},
"zerolinecolor": "white",
"zerolinewidth": 2
}
}
},
"xaxis": {
"anchor": "y",
"domain": [
0.0,
1.0
]
},
"yaxis": {
"anchor": "x",
"domain": [
0.0,
1.0
]
}
}
},
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"M = main(100, 0, 1, method=2)\n",
"import plotly.express as px\n",
"z = np.where(abs(M.diags().toarray()) > 1e-6, 0, 1).astype(bool)\n",
"fig = px.imshow(z, binary_string=True)\n",
"fig.show()\n",
"#plt.spy(M.diags(), markersize=0.2) # or use matplotlib"
]
},
{
"cell_type": "markdown",
"id": "07f47d49",
"metadata": {
"editable": true
},
"source": [
"The coefficient matrix has 4 non-zero diagonals. You can now experiment\n",
"with different test and trial functions, but you will not get a better\n",
"result than that. Try basis ([5](#eq:heinrichs)) for both test and trial\n",
"function, and you'll get 5 nonzero diagonals.\n",
"\n",
"To see the convergence rate call `main` for a range of\n",
"different mesh sizes"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "2e593550",
"metadata": {
"collapsed": false,
"editable": true,
"execution": {
"iopub.execute_input": "2024-02-19T13:48:42.540426Z",
"iopub.status.busy": "2024-02-19T13:48:42.540202Z",
"iopub.status.idle": "2024-02-19T13:48:44.199472Z",
"shell.execute_reply": "2024-02-19T13:48:44.198821Z"
}
},
"outputs": [
{
"data": {
"application/vnd.plotly.v1+json": {
"config": {
"plotlyServerURL": "https://plot.ly"
},
"data": [
{
"hovertemplate": "x=%{x}
y=%{y}",
"legendgroup": "",
"line": {
"color": "#636efa",
"dash": "solid"
},
"marker": {
"symbol": "circle"
},
"mode": "lines",
"name": "",
"orientation": "v",
"showlegend": false,
"type": "scatter",
"x": [
16,
64,
256,
1024
],
"xaxis": "x",
"y": [
142.6329066363412,
0.0012020230030466178,
3.6479943647595484e-14,
1.8722985578999504e-14
],
"yaxis": "y"
}
],
"layout": {
"legend": {
"tracegroupgap": 0
},
"margin": {
"t": 60
},
"template": {
"data": {
"bar": [
{
"error_x": {
"color": "#2a3f5f"
},
"error_y": {
"color": "#2a3f5f"
},
"marker": {
"line": {
"color": "#E5ECF6",
"width": 0.5
},
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "bar"
}
],
"barpolar": [
{
"marker": {
"line": {
"color": "#E5ECF6",
"width": 0.5
},
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "barpolar"
}
],
"carpet": [
{
"aaxis": {
"endlinecolor": "#2a3f5f",
"gridcolor": "white",
"linecolor": "white",
"minorgridcolor": "white",
"startlinecolor": "#2a3f5f"
},
"baxis": {
"endlinecolor": "#2a3f5f",
"gridcolor": "white",
"linecolor": "white",
"minorgridcolor": "white",
"startlinecolor": "#2a3f5f"
},
"type": "carpet"
}
],
"choropleth": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "choropleth"
}
],
"contour": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "contour"
}
],
"contourcarpet": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "contourcarpet"
}
],
"heatmap": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "heatmap"
}
],
"heatmapgl": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "heatmapgl"
}
],
"histogram": [
{
"marker": {
"pattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
}
},
"type": "histogram"
}
],
"histogram2d": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "histogram2d"
}
],
"histogram2dcontour": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "histogram2dcontour"
}
],
"mesh3d": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"type": "mesh3d"
}
],
"parcoords": [
{
"line": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "parcoords"
}
],
"pie": [
{
"automargin": true,
"type": "pie"
}
],
"scatter": [
{
"fillpattern": {
"fillmode": "overlay",
"size": 10,
"solidity": 0.2
},
"type": "scatter"
}
],
"scatter3d": [
{
"line": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatter3d"
}
],
"scattercarpet": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattercarpet"
}
],
"scattergeo": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattergeo"
}
],
"scattergl": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattergl"
}
],
"scattermapbox": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scattermapbox"
}
],
"scatterpolar": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterpolar"
}
],
"scatterpolargl": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterpolargl"
}
],
"scatterternary": [
{
"marker": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"type": "scatterternary"
}
],
"surface": [
{
"colorbar": {
"outlinewidth": 0,
"ticks": ""
},
"colorscale": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"type": "surface"
}
],
"table": [
{
"cells": {
"fill": {
"color": "#EBF0F8"
},
"line": {
"color": "white"
}
},
"header": {
"fill": {
"color": "#C8D4E3"
},
"line": {
"color": "white"
}
},
"type": "table"
}
]
},
"layout": {
"annotationdefaults": {
"arrowcolor": "#2a3f5f",
"arrowhead": 0,
"arrowwidth": 1
},
"autotypenumbers": "strict",
"coloraxis": {
"colorbar": {
"outlinewidth": 0,
"ticks": ""
}
},
"colorscale": {
"diverging": [
[
0,
"#8e0152"
],
[
0.1,
"#c51b7d"
],
[
0.2,
"#de77ae"
],
[
0.3,
"#f1b6da"
],
[
0.4,
"#fde0ef"
],
[
0.5,
"#f7f7f7"
],
[
0.6,
"#e6f5d0"
],
[
0.7,
"#b8e186"
],
[
0.8,
"#7fbc41"
],
[
0.9,
"#4d9221"
],
[
1,
"#276419"
]
],
"sequential": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
],
"sequentialminus": [
[
0.0,
"#0d0887"
],
[
0.1111111111111111,
"#46039f"
],
[
0.2222222222222222,
"#7201a8"
],
[
0.3333333333333333,
"#9c179e"
],
[
0.4444444444444444,
"#bd3786"
],
[
0.5555555555555556,
"#d8576b"
],
[
0.6666666666666666,
"#ed7953"
],
[
0.7777777777777778,
"#fb9f3a"
],
[
0.8888888888888888,
"#fdca26"
],
[
1.0,
"#f0f921"
]
]
},
"colorway": [
"#636efa",
"#EF553B",
"#00cc96",
"#ab63fa",
"#FFA15A",
"#19d3f3",
"#FF6692",
"#B6E880",
"#FF97FF",
"#FECB52"
],
"font": {
"color": "#2a3f5f"
},
"geo": {
"bgcolor": "white",
"lakecolor": "white",
"landcolor": "#E5ECF6",
"showlakes": true,
"showland": true,
"subunitcolor": "white"
},
"hoverlabel": {
"align": "left"
},
"hovermode": "closest",
"mapbox": {
"style": "light"
},
"paper_bgcolor": "white",
"plot_bgcolor": "#E5ECF6",
"polar": {
"angularaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"bgcolor": "#E5ECF6",
"radialaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
}
},
"scene": {
"xaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
},
"yaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
},
"zaxis": {
"backgroundcolor": "#E5ECF6",
"gridcolor": "white",
"gridwidth": 2,
"linecolor": "white",
"showbackground": true,
"ticks": "",
"zerolinecolor": "white"
}
},
"shapedefaults": {
"line": {
"color": "#2a3f5f"
}
},
"ternary": {
"aaxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"baxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
},
"bgcolor": "#E5ECF6",
"caxis": {
"gridcolor": "white",
"linecolor": "white",
"ticks": ""
}
},
"title": {
"x": 0.05
},
"xaxis": {
"automargin": true,
"gridcolor": "white",
"linecolor": "white",
"ticks": "",
"title": {
"standoff": 15
},
"zerolinecolor": "white",
"zerolinewidth": 2
},
"yaxis": {
"automargin": true,
"gridcolor": "white",
"linecolor": "white",
"ticks": "",
"title": {
"standoff": 15
},
"zerolinecolor": "white",
"zerolinewidth": 2
}
}
},
"xaxis": {
"anchor": "y",
"domain": [
0.0,
1.0
],
"title": {
"text": "x"
}
},
"yaxis": {
"anchor": "x",
"domain": [
0.0,
1.0
],
"exponentformat": "e",
"showexponent": "all",
"title": {
"text": "y"
},
"type": "log"
}
}
},
"text/html": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"error = []\n",
"N = (2**4, 2**6, 2**8, 2**10)\n",
"for n in N:\n",
" error.append(main(n, 0, 1))\n",
"fig = px.line(x=N, y=error, log_y=True)\n",
"fig.update_layout(yaxis=dict(showexponent='all', exponentformat='e'))\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"id": "3e83de81",
"metadata": {
"editable": true
},
"source": [
""
]
}
],
"metadata": {
"kernelspec": {
"display_name": "shenfun24",
"language": "python3",
"name": "shenfun24"
},
"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.10.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}