2f:T5d3,# Numerical method of relaxation Relaxation is a numerical technique to solve Poisson's equation ${tex.block`\Delta u = f(x,y)`} over a region ${tex`\Omega \subset \mathbb R^n`} with prescribed values on the boundary of ${tex`\Omega`}. It's particularly convenient when ${tex`\Omega`} is a hypercube, i.e. the generalization of a rectangle in ${tex`\mathbb R^2`}. We present here an implementation to solve Laplace's over the unit square in the plane to simply illustrate the technique. Since we only deal with Laplace's equation, we have ${tex`f(x,y) = 0`}. The basic idea is to Since ${tex.block`u_{xx} \approx \frac{u(x+h,y) - 2u(x,y) + u(x-h,y)}{h}`} and ${tex.block`u_{yy} \approx \frac{u(x,y+h) - 2u(x,y) + u(x,y-h)}{h},`} we have ${tex.block`\Delta u \approx \frac{u(x+h,y) + u(x,y+h) - 4u(x,y) + u(x-h,y) + u(x,y-h)}{h}.`} If we set this equal to zero (to approximate Laplace's equation) and solve for ${tex`u(x,y)`} we obtain ${tex.block`u(x,y) \approx \frac{u(x+h,y) + u(x,y+h) + u(x-h,y) + u(x,y-h)}{4}.`} Thus, we expect the value of a harmonic function at a point ${tex`(x,y)`} to be about the average of the four points to the left, right, above, and below ${tex`(x,y)`}. That suggests the following strategy to solve Laplace's equation over the unit square: We first choose ${tex`N\in\mathbb N`} and discretize the square by choosing points ${tex.block`(x_j,y_i) = (j/N,i/N) \text{ where } 0 \leq i,j \leq N.`} This type of discretization is shown in the image below:30:T80a,{ let container = d3.create("div").style("position", "relative"); let N = 10; let pts = d3 .range(N + 1) .map((i) => d3.range(N + 1).map((j) => [j / N, i / N])) .flat(); let w = d3.min([800, width]); let h = d3.min([800, width]); let svg = container .append("svg") .attr("viewBox", [0, 0, w, h]) .attr("width", "100%") .style("max-width", `${w}px`) .style("border", "solid black 2px"); let pad = 20; let [xmin, xmax, ymin, ymax] = [0, 1, 0, 1]; let x_scale = d3 .scaleLinear() .domain([xmin, xmax]) .range([pad, w - pad]); let y_scale = d3 .scaleLinear() .domain([ymin, ymax]) .range([h - pad, pad]); let point_group = svg.append("g"); point_group .selectAll("circle") .data(pts) .join("circle") .attr("cx", (d) => x_scale(d[0])) .attr("cy", (d) => y_scale(d[1])) .attr("r", (d) => d[0] == 0 || d[0] == 1 || d[1] == 0 || d[1] == 1 ? 4 : 2 ); let i = 4; let j = 3; let atom = [ { x: j, y: i, display: tex`u_{i,j}`.outerHTML }, { x: j - 1, y: i, display: tex`u_{i,j-1}`.outerHTML }, { x: j + 1, y: i, display: tex`u_{i,j+1}`.outerHTML }, { x: j, y: i - 1, display: tex`u_{i-1,j}`.outerHTML }, { x: j, y: i + 1, display: tex`u_{i+1,j}`.outerHTML } ]; let atom_group = svg.append("g"); atom_group .selectAll("circle") .data(atom) .join("circle") .attr("cx", (d) => x_scale(d.x / N)) .attr("cy", (d) => y_scale(d.y / N)) .attr("r", (d) => (d.x == 3 && d.y == 4 ? 20 : 25)) .attr("fill", "white") .attr("stroke", "black"); let overlay = container .append("div") .style("position", "absolute") .style("left", "0px") .style("top", "0px"); overlay .selectAll("div") .data(atom) .join("div") .style("position", "absolute") .style( "left", (d) => `${x_scale(d.x / N - (d.x == 3 && d.y == 4 ? 0.015 : 0.025))}px` ) .style("top", (d) => `${y_scale(d.y / N + 0.015)}px`) .html((d) => d.display); return container.node(); }31:T457,We'll associate a value ${tex`u_{i,j}`} with each point. The idea is that ${tex`u_{i,j}`} should approximate ${tex`u(x_j,y_i)`}. Those values are prescribed by the boundary conditions at the larger points on the boundary; they are unknown at the smaller points in the interior. We do know, though, that if ${tex`h=1/N`}, then our averaging formula for the value of ${tex`u(x,y)`} becomes ${tex.block`u_{i,j} \approx \frac{u_{i+1,j} + u_{i,j+1} + u_{i,j+1} + u_{i,j-1}}{4}.`} Thus, we'll assign arbitrary initial values (${tex`u_{i,j} = 0`} throughout should do) and iteratively improve those values by recomputing the value at ${tex`u_{i,j}`} to be the average of its neighbors at each step. In the example below, the boundary conditions are - ${tex`u(0,y) = 4y^3(1-y)`}, - ${tex`u(1,y) = 4y(1-y)^3`}, and - ${tex`u(x,0) = u(x,y) = 0`}. The square is discretized into ${tex`50\times50=2500`} pieces. The initial values of the ${tex`u_{i,j}`}s are set to agree with the boundary conditions but to be zero otherwise. The iterative averaging process is applied 5000 times with the results saved every 50 steps.32:T737,pic = { let color = "#40bfbf"; let k = 0; let quads = []; let lines = []; let cut = Math.floor(N / 10); for (let i = 0; i <= N; i++) { for (let j = 0; j <= N; j++) { let y = -1 + (2 * i) / N; let x = -1 + (2 * j) / N; if (i < N && j < N) { quads.push([k, k + 1, k + N + 2, k + N + 1]); } if (j % cut == 0 && i < N) { lines.push([k, k + N + 1]); } if (i % cut == 0 && j < N) { lines.push([k, k + 1]); } k++; } } let grid = grids[0]; let face_index_string = quads.map((a) => a.toString() + ",-1").toString(); let line_index_string = lines.map((a) => a.toString() + ",-1").toString(); let surface_transform = d3.create("transform").attr("id", "surface"); let surface_shape = surface_transform.append("shape"); let appearance = surface_shape.append("appearance"); let material = appearance.append("material").attr("diffuseColor", color); // .attr("transparency", transparency); let indexedFaceSet = surface_shape .append("indexedFaceset") .attr("coordIndex", face_index_string) .attr("solid", "false") .attr("creaseAngle", "3.14156"); indexedFaceSet.append("coordinate").attr("point", grids[0]); let mesh = surface_transform.append("shape"); mesh.append("appearance").append("material"); mesh .append("IndexedLineSet") .attr("coordIndex", line_index_string) .append("coordinate") .attr("point", grids[0]); return show_x3d([surface_transform.node()], { extent: { xmin: -0.5, xmax: 1.2, ymin: -0.5, ymax: 1.2, zmin: -0.5, zmax: 0.5 }, viewpoint: { position: "2.5732409603459807 1.1740144760644737 0.8426150401498438", orientation: "0.36744975073861025 0.5553225082283657 0.7460546846827583,2.228538657977855" } }); }33:T4e1,### Why "relaxation"? The formulation ${tex.block`\frac{u_{i+1,j} + u_{i,j+1} - 4u_{i,j} + u_{i,j+1} + u_{i,j-1}}{h} = 0`} for ${tex`1\leq i,j < N`}, makes it clear that the system can be expressed as a system of linear equations. We could thus write ${tex.block`M\vec{u} = \vec{0}.`} The iterative steps can be obtained by "relaxing" the requirement that the right hand side be zero. At least so says [this MathOverflow post](https://mathoverflow.net/questions/385228/origin-of-the-term-relaxation-method-in-numerical-analysis-for-iteratively-solvi). ### Other approaches The finite element method, of course, is The Bomb of all numerical methods in PDEs. Its main advantage is the flexibility that it offers in dealing with regions of all shapes. You can play with the finite element method as applied to steady state heat flow using [this Steady State Heat Explorer](https://marksmath.org/classes/Spring2024PDE/demos/SteadyStateHeatExplorer2D.html). It's worth mentioning that Laplace's equation on a square can be solved analytically via separation of variables and Fourier's technique. [This Mathematica notebook](https://www.wolframcloud.com/obj/mcmcclur/Published/LaplaceSeparation.nb) illustrates the technique for the equation here. 34:T5f8,U = { // Initial grid with values determined solely by u0 // let U0 = []; // for (let i = 0; i <= N; i++) { // let y = i / N; // U0.push( // d3.range(N + 1).map(function (j) { // let x = j / N; // return u0(x, y); // }) // ); // } // At the first step, we simply set the value on the boundary let U1 = []; for (let i = 0; i <= N; i++) { let y = i / N; U1.push( d3.range(N + 1).map(function (j) { let x = j / N; if (i == 0 || i == N || j == 0 || j == N) { return uBoundary(x, y); } else { return u0(x, y); } }) ); } // We now take T more steps defining the value at each interior point // to be the average of the values of its neighbors. // Every 50 steps, we push the result to the list of grids. let U = [U1]; for (let k = 1; k < T + 1; k++) { let UU = []; for (let i = 0; i <= N; i++) { let y = i / N; UU.push( d3.range(0, N + 1).map(function (j) { if (i == 0 || i == N || j == 0 || j == N) { let x = j / N; return uBoundary(x, y); } else { return ( (U[k - 1][i - 1][j] + U[k - 1][i + 1][j] + U[k - 1][i][j - 1] + U[k - 1][i][j + 1]) / 4 ); } }) ); } U.push(UU); } let UThin = [U1]; for (let i = 50; i < U.length; i = i + 50) { UThin.push(U[i]); } return UThin; }35:[] 2d:{"title":"Numerical method of relaxation","cells":[{"id":0,"value":"$2f","pinned":false,"hidden":false,"mode":"md"},{"id":289,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":399,"value":"$31","pinned":false,"hidden":false,"mode":"md"},{"id":158,"value":"viewof t = Inputs.range([0, grids.length - 1], {\n step: 1,\n value: 0,\n label: \"Steps (by 50):\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":83,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":417,"value":"$33","pinned":false,"hidden":false,"mode":"md"},{"id":248,"value":"## Code","pinned":false,"hidden":false,"mode":"md"},{"id":251,"value":"### Boundary and initial conditions","pinned":false,"hidden":false,"mode":"md"},{"id":9,"value":"// Boundary condition\nfunction uBoundary(x, y) {\n if (x == 0) {\n return 4 * y ** 3 * (1 - y);\n } else if (y == 0) {\n return 0;\n } else if (x == 1) {\n return 4 * y * (1 - y) ** 3;\n } else if (y == 1) {\n return 0;\n }\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":6,"value":"// Initial condition\nu0 = (x, y) => 0","pinned":true,"hidden":false,"mode":"ojs"},{"id":255,"value":"### The action\n\nHere's where the main algorithm is implemented.","pinned":false,"hidden":false,"mode":"md"},{"id":11,"value":"$34","pinned":true,"hidden":false,"mode":"ojs"},{"id":268,"value":"### Javascripty stuff","pinned":false,"hidden":false,"mode":"md"},{"id":273,"value":"#### Translate the U grids into something that X3Dom understands","pinned":false,"hidden":false,"mode":"md"},{"id":143,"value":"grids = {\n return d3.range(0, U.length).map(function (t) {\n let grid = [];\n let k = 0;\n for (let i = 0; i <= N; i++) {\n for (let j = 0; j <= N; j++) {\n let y = i/ N;\n let x = j / N;\n grid.push([x, y, U[t][i][j]]);\n k++;\n }\n }\n return pt_list_to_coordString(grid);\n });\n}","pinned":true,"hidden":false,"mode":"ojs"},{"id":275,"value":"#### React to input value of `t`","pinned":false,"hidden":false,"mode":"md"},{"id":168,"value":"d3.select(pic).selectAll(\"coordinate\").attr(\"point\", grids[t])","pinned":true,"hidden":false,"mode":"ojs"},{"id":278,"value":"#### Parameters","pinned":false,"hidden":false,"mode":"md"},{"id":145,"value":"T = 5000","pinned":true,"hidden":false,"mode":"ojs"},{"id":4,"value":"N = 50","pinned":true,"hidden":false,"mode":"ojs"},{"id":280,"value":"#### Imports","pinned":false,"hidden":false,"mode":"md"},{"id":74,"value":"import { pt_list_to_coordString, show_x3d } from \"@mcmcclur/parametric-surfaces\"","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W35","stdlib":"1","titling":"auto"}