2a:null 2b:["$","title",null,{"children":"Laplace Equation | Renato Paes Leme | Observable"}] 30:T744,md` ### Rubber sheet Let's look at a very concrete example: suppose you have a piece of rubber sheet in the shape of the set ${tex`U`}. Initially the rubber piece is on the plane, but then you get its boundary (let's call ${tex`\partial U`}) and fix its height (${tex`z`}-coordinate). For each point ${tex`(x,y) \in U`} define ${tex`u(x,y)`} to be its ${tex`z`}-coordinate. To see how ${tex`u`} will look like, we look at a little square ${tex`S_0 = [x,x+ \delta] \times [y,y+\delta]`} and we can see the forces acting on it. It is useful to view the rubber sheet as a bunch of little ${tex`\delta \times \delta`} squares connected by [springs](https://en.wikipedia.org/wiki/Spring_(device%29). According to [Hooke's Law](https://en.wikipedia.org/wiki/Hooke%27s_law) the force is proportional to the amount it is being stretched. For example, the square ${tex`S_x^+ = [x+\delta,x+ 2 \delta] \times [y,y+\delta]`} exerts a force on ${tex`S_0`} proportional to ${tex`u(x+\delta,y) - u(x,y) \approx \delta \frac{\partial u}{\partial x} (x+\delta, y)`}. The force exerted by the square ${tex`S_x^- = [x-\delta,x] \times [y,y+\delta]`} is given by ${tex`u(x-\delta,y) - u(x-,y) \approx - \delta \frac{\partial u}{\partial x} (x, y)`}. We can do the same with the squares ${tex`S_y^+`} and ${tex`S_y^-`}. In the end, the total force on the square ${tex`S_0`} will be: ${tex.block`\sum F = \delta \frac{\partial u}{\partial x} (x+\delta, y) - \delta \frac{\partial u}{\partial x} (x, y) + \delta \frac{\partial u}{\partial y} (x, y+\delta) - \delta \frac{\partial u}{\partial y} (x, y)`} which can be approximated by: ${tex.block`\delta^2 \frac{\partial^2 u}{\partial x^2} + \delta \frac{\partial^2 u}{\partial y^2} = \delta^2 \Delta u(x,y)`} Since the system is equilibrium the sum of forces in each square should be zero, so ${tex`0 = \sum F = \Delta u(x,y)`}. `31:T66c,md` ### Holomorphic Functions Another important application of Laplace's equation is in the study of [holomorphic functions](https://en.wikipedia.org/wiki/Holomorphic_function) in Complex Analysis. An holomorphic function is a map ${tex`f: U \subseteq \mathbb{C} \rightarrow \mathbb{C}`} that is differentiable, i.e, the following limit exists for every point ${tex`z \in U`}: ${tex.block`f'(z) = \lim_{\vert w \vert \rightarrow 0} \frac{f(z+w) - f(z)}{w}`} or in other words, for ${tex`w \in \mathbb{C}`} of tiny norm ${tex`\vert w \vert`}, the expression ${tex`(f(z+w)-f(z))/w`} (where the / is a complex division) is close to ${tex`f'(z)`}. One can show that a function ${tex`f`} is holomorphic iff its components satisfy the [Cauchy-Riemann equations](https://en.wikipedia.org/wiki/Cauchy%E2%80%93Riemann_equations). Let the components be functions ${tex`u,v : U \rightarrow \mathbb{R}`} such that: ${tex.block`f(z) = u(z) + i v(z)`} where ${tex`z=x+yi \in \mathbb{C}`}. Often it will be useful to identify ${tex`z`} with ${tex`(x,y) \in \mathbb{R}^2`}. With that identification we can write the Cauchy-Riemann equations as follows: ${tex.block`\frac{\partial v}{\partial x} = - \frac{\partial u}{\partial y} \qquad \frac{\partial v}{\partial y} = \frac{\partial u}{\partial x}`} From those equations it is easy to check that each of the components is a solution to Laplace's equation: ${tex.block`\Delta v(x)= \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} = \frac{\partial}{\partial x} \left[ -\frac{\partial u}{\partial y} \right] + \frac{\partial}{\partial y} \left[ \frac{\partial u}{\partial x} \right] = 0`} `32:T42b,md`Given a certain set in the plane, we will construct an holomorphic bijection between this set and the ball. Let's start by drawing a simple curve around this point (make sure the curve contains the point and that it has no self-loops. Let's call this region ${tex`U`}. We will try to find a map: ${tex.block`u : U \rightarrow \mathbb{R}`} to be the real part of an holomorphic function we want to construct. So a good thing is to define it as a solution of a Laplace equation: ${tex.block`\Delta u(z) = 0 \quad \text{for } z \in U`} and we will add the following boundary condition: ${tex.block`u(z) = -\log \vert z - z_0 \vert \quad \text{for } z \in \partial U`} where ${tex`\vert z - z_0 \vert`} is the norm in ${tex`\mathbb{C}`}, i.e. the norm when we view ${tex`z - z_0`} as a complex number. This problem has an unique solution by the [Dirichlet Principle](https://en.wikipedia.org/wiki/Dirichlet%27s_principle). After you finish drawing the boundary of ${tex`U`} (containing the point) you will see the solution of the Laplace equation for that boundary.`33:Td8f,{ const height = 320; const width = 640; const context = DOM.context2d(width, height, 1); var frame; context.strokeStyle = "red"; var active = false; var canvas = d3.select(context.canvas); var points = []; function dist(q, r) { return Math.sqrt( (q[0] - r[0]) * (q[0] - r[0]) + (q[1] - r[1]) * (q[1] - r[1]) ); } function reset(dot = true) { context.clearRect(0, 0, width, height); if (dot) { context.fillStyle = "blue"; context.fillRect(200,150,5,5); } context.beginPath(); } reset(); function updatepoint(p,q,r,k) { return [p[0] + q[0]/k - r[0]/k, p[1] + q[1]/k - r[1]/k]; } function mouseup() { reset(false); var n = points.length; var k = Math.min(12, n); var q = [0,0] for (var i = 0; i < k; ++i) { q = updatepoint(q, points[i], [0,0], k); } var upoints = []; context.moveTo(q[0], q[1]); for (var i = 1; i <= n; ++i) { q = updatepoint(q, points[(i+k)%n], points[i%n], k); context.lineTo(q[0], q[1]); upoints.push([q[0], q[1]]); context.stroke(); } //context.fill(); var img = context.getImageData(0, 0, width, height); var boundaryPoints = []; var interior = PointsInside(img, [200, 150], width, height, boundaryPoints); var boundary = q => -Math.log(dist(q, [200, 150])); function callback(sol) { FillPoints(interior, sol, context); } var solver = MultiResLaplacianSolver(interior, boundary); var res = 64; while (!solver.SetResolution(res)) { res /= 2; } var solution, error, iterations = 0; async function tick() { [solution, error] = solver.Step(20); iterations += 20; console.log("Iteration " + iterations + ": error = " + error); if ((error > 0.001 || res > 1) && iterations < 10000 ) { if (error < 0.001) { res = res / 2; while (!solver.SetResolution(res)) { res /= 2; } } frame = requestAnimationFrame(tick); } else { // define the closed region composed of interior and boundary and extend the // solution to include the boundary conditions var closure = interior.concat(boundaryPoints); var uBar = solution.concat(boundaryPoints.map(p=>boundary(p))); viewof LaplaceData.setValue("closure", closure); viewof LaplaceData.setValue("interior", interior); viewof LaplaceData.setValue("boundary", boundary); viewof LaplaceData.setValue("boundaryPoints", boundaryPoints); viewof LaplaceData.setValue("points", upoints); viewof LaplaceData.setValue("u", solution); viewof LaplaceData.setValue("uBar", uBar); viewof LaplaceData.setValue("z0", [200, 150]); viewof LaplaceData.startEvent("e1"); } FillPoints(interior, solution, context); } tick(); active = false; points = []; } function mousedown() { active = true; reset(); var mousePos = d3.mouse(canvas.node()); context.moveTo(mousePos[0], mousePos[1]); points = [[mousePos[0], mousePos[1]]]; } function mousemove() { if (active) { var mousePos = d3.mouse(canvas.node()); context.lineTo(mousePos[0], mousePos[1]); context.stroke(); points.push([mousePos[0], mousePos[1]]); } } canvas .on("mousemove", mousemove) .on("mouseup", mouseup) .on("mousedown", mousedown); return context.canvas; }34:T61e,{ const height = 320; const width = 640+320; const context = DOM.context2d(width, height, 1); function update() { context.clearRect(0, 0, width, height); var points = viewof LaplaceData.value.points; var closure = viewof LaplaceData.value.closure; context.strokeStyle = "red"; var n = points.length context.beginPath(); context.moveTo(points[n-1][0], points[n-1][1]); for (var i = 0; i < n; ++i) { context.lineTo(points[i][0], points[i][1]); context.stroke(); } var u = viewof LaplaceData.value.uBar; var v = viewof LaplaceData.value.vBar; var z0 = viewof LaplaceData.value.z0; var k = 5; var n = closure.length var xc = Math.floor(BoundingBox(points)[0][1]+200); function CirclePos(w) { return [Math.floor(xc + 100 * w[0]), Math.floor(150 + 100 * w[1])]; } context.strokeStyle = "red" context.beginPath(); var c0 = CirclePos([1,0]); context.moveTo(c0[0], c0[1]); var T = 200; for (var t = 0; t <= T; ++t) { var w = Complex.exp([0,(t * Math.PI * 2) /T]); var wr = CirclePos(w); context.lineTo(wr[0], wr[1]); } context.stroke(); for (var t = 0; t < n; ++t) { var z = closure[t]; if (z[0]%k == 0 && z[1]%k == 0) { context.fillRect(z[0], z[1], 1, 1); var w = ConformalMap(z, z0, u[t], v[t]); var wr = CirclePos(w); context.fillRect(wr[0], wr[1], 1, 1); } } } viewof LaplaceData.setEvent("e2", update); return context.canvas; }35:Tfeb,{ const height = 320; const width = 640+320; const context = DOM.context2d(width, height, 1); function update() { var closure = viewof LaplaceData.value.closure; var m = closure.length var index = {} for (var t = 0; t < m; ++t) { index[closure[t]] = t; } context.clearRect(0, 0, width, height); var points = viewof LaplaceData.value.points; var closure = viewof LaplaceData.value.closure; context.strokeStyle = "red"; var n = points.length context.beginPath(); context.moveTo(points[n-1][0], points[n-1][1]); for (var i = 0; i < n; ++i) { context.lineTo(points[i][0], points[i][1]); context.stroke(); } var u = viewof LaplaceData.value.uBar; var v = viewof LaplaceData.value.vBar; var z0 = viewof LaplaceData.value.z0; var k = 5; function neighbors(p) { return [[p[0] - k, p[1]], [p[0] + k, p[1]], [p[0], p[1] - k], [p[0], p[1] + k]]; } var n = closure.length function AddSegment(p1, p2, color) { context.strokeStyle = color context.beginPath(); context.moveTo(p1[0], p1[1]); context.lineTo(p2[0], p2[1]); context.stroke(); } var xc = Math.floor(BoundingBox(points)[0][1]+200); function CirclePos(w) { return [Math.floor(xc + 100 * w[0]), Math.floor(150 + 100 * w[1])]; } context.strokeStyle = "red" context.beginPath(); var c0 = CirclePos([1,0]); context.moveTo(c0[0], c0[1]); var T = 200; for (var t = 0; t <= T; ++t) { var w = Complex.exp([0,(t * Math.PI * 2) /T]); var wr = CirclePos(w); context.lineTo(wr[0], wr[1]); } context.stroke(); function Cell(z, d1, d2) { function SearchBoundary(zi, di) { for (var ii = k; ii >= 0; --ii) { var z1 = [zi[0] + ii * di[0], zi[1] + ii * di[1]]; if (z1 in index) { return [z1, ii]; } } return [zi, 0]; } function Minus(d) { return [-d[0], -d[1]]; } var z1, i1, z2, i2, z11, z21; var zc = [z[0], z[1]]; z11 = null; z21 = null; [z1, i1] = SearchBoundary(z, d1); [z2, i2] = SearchBoundary(z, d2); var out = [zc, z1]; var it; if (i1 == k) { [z11, it] = SearchBoundary(z1, d2); out.push(z11); if (it == k && i2 < k) { [z21, it] = SearchBoundary(z11, Minus(d1)); out.push(z21); } } if (i2 == k) { [z21, it] = SearchBoundary(z2, d1); if (it == k && i1 < k) { [z11, it] = SearchBoundary(z21, Minus(d2)); out.push(z11); } out.push(z21); } out.push(z2); return out; } function PointsToSegments(pts) { var segs = []; for (var i = 0; i < pts.length; ++i) { segs.push([pts[i], pts[(i+1)%pts.length]]); } return segs; } function CellsAround(z) { var segs = []; var dirs = [[1,0], [0,1], [-1,0], [0,-1]]; for (var i = 0; i < 4; ++i) { segs = segs.concat(PointsToSegments(Cell(z, dirs[i], dirs[(i+1)%4]))); } return segs; } for (var t = 0; t < n; ++t) { var z = closure[t]; if (z[0]%k == 0 && z[1]%k == 0) { context.fillRect(z[0], z[1], 1, 1); var w = ConformalMap(z, z0, u[t], v[t]); var wr = CirclePos(w); //context.fillRect(wr[0], wr[1], 1, 1); CellsAround(z).forEach(function(s) { var z1 = s[0]; var t1 = index[z1]; var z2 = s[1]; var t2 = index[z2]; var w1 = ConformalMap(z1, z0, u[t1], v[t1]); var w2 = ConformalMap(z2, z0, u[t2], v[t2]); var wr1 = CirclePos(w1); var wr2 = CirclePos(w2); AddSegment(z1, z2, "MediumTurquoise"); AddSegment(wr1, wr2, "MediumTurquoise"); }); } } } viewof LaplaceData.setEvent("e3", update); return context.canvas; }36:Td4a,md` ### Solving Laplace's equation In the previous section we solved a [Dirichlet problem](https://en.wikipedia.org/wiki/Dirichlet_problem) which is the problem of finding a solution to Laplace's equation in a region when we fix the function on the boundary. There are various elegant analytic solutions for this problem, but above we solve it numerically. The way we do is that we discretize the space in a grid. It is convenient that each pixel in the screen is already a square, so our variables will be the values of the function ${tex`u(x)`} in a certain pixel ${tex`x`}. There are finitely many pixels, so this is a finite problem we will be solving. For the pixels in the boundary we set the value equal to the boundary condition. For the values inside, we add a constraint like that: ${tex.block`4 u(x,y) - u(x-1, y) - u(x+1, y) - u(x, y-1) - u(x, y+1) = 0`} We do that since a good approximation to the operator ${tex`\Delta u(x)`} is ${tex.block`\frac{1}{\delta^2} \left[4 u(x,y) - u(x - \delta, y) - u(x + \delta, y) - u(x, y - \delta) - u(x, y + \delta) \right]`} After we set those conditions, we end up with a linear system of the form: ${tex.block`Lu = b`} where ${tex`u`} is a column vector with one entry for every pixel inside the shape and ${tex`L`} is a matrix encoding the discrete version of the Laplace operator. The vector ${tex`b`} will collect the terms corresponding to values set in the boundary. Now, our problem is solving a linear system. The first thing that comes to mind is [Gaussian elimination](https://en.wikipedia.org/wiki/Gaussian_elimination) which can take ${tex`O(N^3)`} operations where ${tex`N`} is the number of pixels. This is a bit too much. Laplacian systems are special since they are sparse (most entries in ${tex`L`} are zero) and they have very nice properties. First, let's try to exploit the fact that ${tex`L`} is very sparse so for any vector ${tex`u \in \mathbb{R}^N`} we can evaluate ${tex`Lu`} in ${tex`O(N)`} operations. We will start with a vector ${tex`u_0`} and try to bring it closer to the solution. For that, we observe that the solution of ${tex`Lu=b`} is the minimizer of the function: ${tex.block`\phi(u) := \frac{1}{2} \langle u, Lu \rangle - \langle b, u \rangle`} so a natural idea is to do gradient descent on ${tex`\phi`}. We start with a certain solution, say ${tex`u_0 = 0`} and then we do a step in the direction of: ${tex.block`-\nabla \phi(u_t) := b - L u_t`} We usually call this the residue ${tex`r_t = b - L u_t`}. Now, we want to set ${tex.block`u_{t+1} = u_t + \eta_t r_t`} What value of ${tex`\eta_t`} to set ? We can look at the function: ${tex.block`\begin{aligned}\eta \mapsto \phi(u_t + \eta r_t) & = \frac{1}{2} \langle u_t + \eta r_t, L(u_t + \eta r_t) \rangle - \langle b, u_t + \eta r_t \rangle = \\ & = \phi(u_t) + \frac{\eta^2}{2} \langle r_t, L r_t \rangle - \eta \langle r_t, r_t \rangle \end{aligned}`} which a second-order algebraic equation in ${tex`\eta`}, so its minimizer is: ${tex.block`\eta_t = \frac{\langle r_t, r_t \rangle}{\langle r_t, L r_t \rangle}`} We repeat this iteration until the error is small enough, i.e. ${tex`\Vert L u_t - b \Vert_2 \leq \epsilon`}. This method is known as the [Richardson iteration](https://en.wikipedia.org/wiki/Modified_Richardson_iteration). If you click Play below, you will see the equation being solved applying the Richardson iteration. `37:T7de,{ await visibility(); const height = 230; const width = 640; const context = DOM.context2d(width, height, 1); var frame; function dist(q, r) { return Math.sqrt( (q[0] - r[0]) * (q[0] - r[0]) + (q[1] - r[1]) * (q[1] - r[1]) ); } var points = testPoints.map(p => [Math.floor(p[0]), Math.floor(p[1])]); var box = BoundingBox(points); var points = points.map(p => [40 + p[0] - box[0][0], 40 + p[1] - box[1][0]]); var point0 = [40 + testPoint0[0] - box[0][0], 40 + testPoint0[1] - box[1][0]]; context.clearRect(0, 0, width, height); function draw(xi) { context.strokeStyle = "red"; var pointsi = points.map(p => [p[0] + xi, p[1]]); var n = testPoints.length; context.beginPath(); context.moveTo(pointsi[n - 1][0], pointsi[n - 1][1]); for (var i = 0; i < n; ++i) { context.lineTo(pointsi[i][0], pointsi[i][1]); context.stroke(); } } draw(0); //draw(200); var img = context.getImageData(0, 0, width, height); var boundaryPoints = []; var interior = PointsInside(img, point0, width, height, boundaryPoints); console.log(interior); var boundary = q => -Math.log(dist(q, point0)); var solver = MultiResLaplacianSolver(interior, boundary); //var solverCG = MultiResLaplacianSolver(interior, boundary, true); solver.SetResolution(1); //solverCG.SetResolution(1); var solution, error = 1, iterations = 0; //var solutionCG, errorCG=1, iterationsCG = 0; async function tick() { if (error > 0.001) { [solution, error] = solver.Step(20); FillPoints(interior, solution, context); } /* if (errorCG > 0.001) { [solutionCG, errorCG] = solverCG.Step(20); FillPoints(interior, solutionCG, context, 200, 0); } */ iterations += 20; console.log("Iteration " + iterations + ": error = " + error); if (error > 0.001 && iterations < 10000) { frame = requestAnimationFrame(tick); } } if (button1) { tick(); } return context.canvas; }38:T7c1,{ await visibility(); const height = 230; const width = 840; const context = DOM.context2d(width, height, 1); var frame; function dist(q, r) { return Math.sqrt( (q[0] - r[0]) * (q[0] - r[0]) + (q[1] - r[1]) * (q[1] - r[1]) ); } var points = testPoints.map(p=>[Math.floor(p[0]), Math.floor(p[1])]); var box = BoundingBox(points); var points = points.map(p=>[40+p[0]-box[0][0], 40+p[1]-box[1][0]]); var point0 = [40+testPoint0[0]-box[0][0], 40+testPoint0[1]-box[1][0]]; context.clearRect(0, 0, width, height); function draw(xi, yi) { context.strokeStyle = "red"; var pointsi = points.map(p=>[p[0]+xi, p[1]+yi]); var n = testPoints.length context.beginPath(); context.moveTo(pointsi[n-1][0], pointsi[n-1][1]); for (var i = 0; i < n; ++i) { context.lineTo(pointsi[i][0], pointsi[i][1]); context.stroke(); } } var dp = [[0,0], [160,0], [320,0], [480,0]]; // [0,230], [200,230], [400,230]]; for (var i = 0; i < 4; ++i) { draw(dp[i][0], dp[i][1]); } var img = context.getImageData(0, 0, width, height); var boundaryPoints = []; var interior = PointsInside(img, point0, width, height, boundaryPoints); console.log(interior); var boundary = q => -Math.log(dist(q, point0)); var solver = MultiResLaplacianSolver(interior, boundary); var solution, error, iterations = 0; var xi = 0; var res = 32; solver.SetResolution(res); async function tick() { [solution, error] = solver.Step(20); iterations += 20; console.log("Iteration " + iterations + ": error = " + error); FillPoints(interior, solution, context, dp[xi][0], dp[xi][1]); if (error > 0.001 && iterations < 10000) { frame = requestAnimationFrame(tick); } else { if (res > 4) { res /= 2; xi += 1; solver.SetResolution(res); frame = requestAnimationFrame(tick); } } } if (button2) { tick(); } return context.canvas; }39:T1073,function MultiResLaplacianSolver(inside, boundaryfunc) { var n = inside.length; var index = {}; for (var i = 0; i < n; ++i) { index[inside[i]] = i; } // returns [A, b, pullback] function BuildSystem(res) { var adj = []; var b = []; function Neighbors(p, k) { return [[p[0] - k, p[1]], [p[0] + k, p[1]], [p[0], p[1] - k], [p[0], p[1] + k]]; } function FullSquare(p) { for (var i = 0; i < res; ++i) { for (var j = 0; j < res; ++j) { if (index[[p[0]+i, p[1]+j]] == undefined) { return false; } } } return true; } function HasOneNeighbor(p) { var hasNeighbor = false; Neighbors(p,1).forEach(function(q) { if (index[q] != undefined) { hasNeighbor = true; } }); return hasNeighbor; } function BoundaryValue(p) { var total = 0; var cells = 0; for (var i = 0; i < res; ++i) { for (var j = 0; j < res; ++j) { var q = [p[0]+i, p[1]+j]; if (index[q] == undefined && HasOneNeighbor(q)) { cells += 1; total += boundaryfunc(q); } } } return total / cells; } // maps (x,y) points to coordinates in the lower resolution vector var coord = {}; var boundaryVals = {}; var lowResPoints = []; var m = 0; for (var i = 0; i < n; ++i) { if ((inside[i][0] % res == 0) && (inside[i][1] % res == 0) && FullSquare(inside[i])) { coord[inside[i]] = m; lowResPoints.push(inside[i]); b.push(0); adj.push([]); m += 1; } } for (var i = 0; i < n; ++i) { if ((inside[i][0] % res == 0) && (inside[i][1] % res == 0) && FullSquare(inside[i])) { var j = coord[inside[i]]; Neighbors(inside[i], res).forEach(function(q) { if (coord[q] != undefined) { adj[j].push(coord[q]); } else { var f = BoundaryValue(q); b[j] += f; boundaryVals[q] = f; } }); } } function A(x) { var y = []; for (var i = 0; i < m; ++i) { y.push( 4 * x[i] ); adj[i].forEach(function(j) { y[i] -= x[j]; }); } return y; } // Maps a low resolution vector x \in \R^m to a high resolution vector x \in \R^n function PullBack(x) { var z = []; for (var i = 0; i < n; ++i) { var q = [inside[i][0] - (inside[i][0] % res), inside[i][1] - (inside[i][1] % res)]; if (coord[q] != undefined) { z.push(x[coord[q]]); } else { if (boundaryVals[q] == undefined) { boundaryVals[q] = BoundaryValue(q); } z.push(boundaryVals[q]); } } return z; } // Maps a high resolution vector z \in \R^n to a low resolution vector x \in \R^m function PushForward(z) { var x = []; for (var t = 0; t < m; ++t) { x.push(0); for (var i = 0; i < res; ++i) { for (var j = 0; j < res; ++j) { var p = [lowResPoints[t][0] + i, lowResPoints[t][1] + j]; x[t] += z[index[p]] / (res * res); } } } return x; } return [A, b, m, PullBack, PushForward]; } var A, b, PullBack, PushForward, m, solver, error; var sol_m = undefined; // low resolution solution var sol_n = undefined; // high resolution solution function SetResolution(res) { [A, b, m, PullBack, PushForward] = BuildSystem(res); solver = LinearSolver(A, b, m); if (m <= 1) return false; console.log("m = " + m); if (sol_n == undefined) { sol_m = solver.AvgInitializer(); sol_n = PullBack(sol_m); } else { sol_m = PushForward(sol_n); } return true; } function Step(k) { if (m == 0) { sol_n = PullBack([]); return [sol_n, 0]; } [sol_m, error] = solver.Step(sol_m, 100); sol_n = PullBack(sol_m); return [sol_n, error]; } return { Step: Step, SetResolution: SetResolution }; }3a:T5b0,// A is a linear function mapping n-vectors to n-vectors, // b is an n-vector and epsilon is the tolerance, our goal is to get x // such that || Ax - b ||_2 \leq \epsilon/. // The code below implements the Richardson iteration function LinearSolver(A, b, n) { // Some auxiliary functions function dot(u,v) { var prod = 0; for (var i = 0; i < u.length; ++i) { prod += u[i] * v[i]; } return prod; } function lin(a, u, b, v) { var comb = []; for (var i = 0; i < u.length; ++i) { comb.push(a * u[i] + b * v[i]); } return comb; } function norm(u) { var l = 0; for (var i = 0; i < u.length; ++i) { l += u[i] * u[i]; } return Math.sqrt(l); } function avg(u) { var sum = 0; for (var i = 0; i < u.length; ++i) { sum += u[i]; } return sum / u.length; } var mean = avg(b); var r, eta, error; function AvgInitializer() { return d3.range(n).map(u => mean); } function Step(x, k) { for (var t = 0; t < k; ++t) { // We take a gradient step in || Ax-b ||_2^2 and do exact optimization // in the direction of the gradient. r = lin(1, b, -1, A(x)); error = norm(r); if (error > 0.00001) { eta = dot(r,r) / dot(r, A(r)); x = lin(1, x, eta, r); } else { return [x, error]; } } return [x, error]; } return { Step : Step, AvgInitializer: AvgInitializer, }; }3b:T704,function ComplexComponent(u, points) { var n = points.length var index = {}; for (var i = 0; i < n; ++i) { index[points[i]] = i; } var v = d3.range(n).map(u => undefined); function neighbors(p) { return [[p[0] - 1, p[1]], [p[0] + 1, p[1]], [p[0], p[1] - 1], [p[0], p[1] + 1]]; } function grad_u(p) { var g = [0,0]; var dx = 0; var dt = 0; var ppdt = [p[0]+1, p[1]] if (ppdt in index) { dx += u[index[ppdt]]; dt += 1.0; } else { dx += u[index[p]]; dt += 0; } var pmdt = [p[0]-1, p[1]] if (pmdt in index) { dx -= u[index[pmdt]]; dt += 1.0; } else { dx -= u[index[p]]; dt += 0; } g[0] = dx / dt; var dy = 0; dt = 0; ppdt = [p[0], p[1]+1] if (ppdt in index) { dy += u[index[ppdt]]; dt += 1.0; } else { dy += u[index[p]]; dt += 0; } pmdt = [p[0], p[1]-1] if (pmdt in index) { dy -= u[index[pmdt]]; dt += 1.0; } else { dy -= u[index[p]]; dt += 0; } g[1] = dy/dt; return g; } function grad_v(p) { var g = grad_u(p); // The relation between u and v is given the the equations of Cauchy-Riemann // https://en.wikipedia.org/wiki/Cauchy%E2%80%93Riemann_equations return [-g[1], g[0]]; } // Returns v(p) assuming v(q) is an already computed neighbor function eval_v(p, q) { var vq = v[index[q]]; var dp = [p[0] - q[0], p[1] - q[1]]; var g = grad_v(q); return vq + dp[0] * g[0] + dp[1] * g[1]; } v[0] = 0; for (var i = 1; i < n; ++i) { neighbors(points[i]).forEach(function(q) { if (q in index && v[index[q]] != undefined) { v[i] = eval_v(points[i], q); } }); } return v; }3c:T184a,testPoints = [[152.49999999999997,169.53645833333334],[152.58333333333331,170.61979166666669],[152.66666666666666,171.53645833333334],[152.75,172.53645833333334],[152.83333333333334,173.61979166666669],[152.91666666666669,174.61979166666669],[152.91666666666669,175.53645833333334],[152.91666666666669,176.28645833333334],[152.91666666666669,177.03645833333334],[152.91666666666669,177.703125],[152.91666666666669,178.28645833333334],[152.91666666666669,178.78645833333334],[152.91666666666669,179.36979166666669],[152.91666666666669,179.95312500000003],[152.83333333333334,180.70312500000003],[152.66666666666669,181.3697916666667],[152.41666666666669,182.03645833333337],[152.08333333333334,182.70312500000006],[151.75,183.45312500000006],[151.33333333333334,184.2864583333334],[150.91666666666669,185.11979166666674],[150.50000000000003,185.95312500000009],[150.00000000000003,186.95312500000009],[149.50000000000003,188.11979166666674],[148.91666666666669,189.2864583333334],[148.25000000000003,190.70312500000006],[147.58333333333337,192.1197916666667],[147.00000000000003,193.45312500000003],[146.41666666666669,194.78645833333337],[145.83333333333337,196.20312500000003],[145.25000000000006,197.6197916666667],[144.6666666666667,199.03645833333337],[144.08333333333337,200.36979166666669],[143.41666666666669,201.86979166666669],[142.75000000000003,203.28645833333334],[142.00000000000003,204.61979166666669],[141.25000000000003,206.03645833333334],[140.58333333333334,207.28645833333334],[139.91666666666669,208.45312500000003],[139.25000000000003,209.70312500000003],[138.6666666666667,210.86979166666669],[138.08333333333337,211.95312500000003],[137.50000000000003,212.95312500000003],[136.91666666666669,213.95312500000003],[136.33333333333334,215.03645833333337],[135.83333333333334,216.03645833333337],[135.41666666666666,217.03645833333337],[135.08333333333331,218.03645833333337],[134.83333333333331,218.8697916666667],[134.58333333333331,219.78645833333337],[134.49999999999997,220.61979166666669],[134.49999999999997,221.45312500000003],[134.49999999999997,222.20312500000003],[134.74999999999997,222.95312500000003],[135.08333333333331,223.78645833333334],[135.66666666666666,224.53645833333334],[136.41666666666666,225.203125],[137.41666666666666,225.78645833333331],[138.66666666666666,226.28645833333331],[140.16666666666666,226.70312499999997],[142,227.11979166666663],[144.08333333333334,227.20312499999997],[146.33333333333334,227.20312499999997],[148.58333333333334,227.03645833333331],[151,226.86979166666666],[153.33333333333334,226.53645833333331],[155.75,225.953125],[158.08333333333334,225.28645833333334],[160.33333333333334,224.45312500000003],[162.50000000000003,223.45312500000003],[164.50000000000003,222.28645833333337],[166.4166666666667,221.03645833333337],[168.25000000000006,219.53645833333337],[170.00000000000006,218.03645833333337],[171.8333333333334,216.45312500000006],[173.8333333333334,214.70312500000006],[175.91666666666674,212.78645833333337],[178.0833333333334,210.78645833333337],[180.3333333333334,208.70312500000003],[182.75000000000006,206.53645833333337],[185.3333333333334,204.28645833333337],[188.00000000000006,202.03645833333337],[190.75000000000006,199.70312500000003],[193.58333333333337,197.28645833333337],[196.33333333333337,194.95312500000006],[199.1666666666667,192.6197916666667],[201.83333333333337,190.20312500000003],[204.50000000000003,187.86979166666669],[207.08333333333334,185.53645833333337],[209.75,183.1197916666667],[212.5,180.70312500000006],[215.33333333333334,178.20312500000006],[218.25,175.70312500000006],[221.33333333333334,173.11979166666674],[224.66666666666669,170.53645833333343],[228.25000000000003,167.95312500000009],[232.16666666666669,165.36979166666674],[236.33333333333334,162.78645833333343],[240.66666666666669,160.28645833333343],[245.25,157.78645833333343],[250.08333333333334,155.28645833333343],[255.08333333333334,152.78645833333343],[260.1666666666667,150.20312500000009],[265.1666666666667,147.70312500000009],[270.0833333333333,145.11979166666674],[274.75,142.5364583333334],[279.1666666666667,139.8697916666667],[283.1666666666667,137.20312500000006],[286.6666666666667,134.45312500000006],[289.6666666666667,131.78645833333337],[292.33333333333337,128.95312500000006],[294.33333333333337,125.95312500000006],[295.75000000000006,122.86979166666673],[296.50000000000006,119.70312500000007],[296.50000000000006,116.61979166666674],[295.91666666666674,113.45312500000009],[294.5833333333334,110.28645833333341],[292.8333333333334,107.20312500000009],[290.41666666666674,104.20312500000009],[287.66666666666674,101.20312500000009],[284.50000000000006,98.28645833333341],[280.91666666666674,95.28645833333341],[276.91666666666674,92.53645833333341],[272.7500000000001,89.95312500000009],[268.41666666666674,87.53645833333341],[263.91666666666674,85.36979166666674],[259.41666666666674,83.45312500000007],[255.00000000000006,81.7864583333334],[250.75000000000006,80.36979166666674],[246.50000000000006,79.11979166666674],[242.41666666666674,78.11979166666674],[238.25000000000006,77.36979166666674],[234.1666666666667,76.95312500000009],[230.0833333333334,76.86979166666674],[226.0833333333334,77.11979166666674],[222.00000000000009,77.7864583333334],[217.91666666666674,78.86979166666674],[213.50000000000009,80.70312500000007],[208.8333333333334,83.20312500000007],[203.6666666666667,86.5364583333334],[198.4166666666667,90.5364583333334],[193.00000000000006,95.0364583333334],[187.5833333333334,99.95312500000006],[182.25000000000006,105.11979166666673],[177.0833333333334,110.36979166666673],[172.1666666666667,115.70312500000007],[167.50000000000006,121.03645833333341],[163.1666666666667,126.2031250000001],[159.1666666666667,131.11979166666677],[156.33333333333337,136.36979166666677],[154.00000000000006,140.9531250000001],[152.50000000000006,144.78645833333343],[151.4166666666667,148.2031250000001],[150.75000000000003,151.2031250000001],[150.41666666666669,153.86979166666677],[150.33333333333334,156.36979166666677],[150.5,158.7031250000001],[150.83333333333334,160.9531250000001],[151.25,163.03645833333346],[151.75,165.1197916666668],[152.25,167.2031250000001],[152.33333333333334,168.36979166666677]];3d:[] 2e:{"id":"4c8175f65642a29d","slug":"laplace-equation","title":"Laplace Equation","thumbnail":"2ceac0d710d26f2baaa66372e1d0fbd3b9c3afbdf6221d21ed998353ef81127f","access_level":"public","create_time":"$D2019-04-21T22:57:38.109Z","update_time":"$D2019-09-09T02:05:49.876Z","view_time":null,"first_public_time":"$D2019-08-10T14:26:58.641Z","fork":false,"stars":25,"starred":false,"owner":{"type":"team","id":"cde07d8c7e702472","login":"renatoppl","name":"Renato Paes Leme","avatar_url":"https://avatars-next.observableusercontent.com/c22df5c297aca2188beaa018c7c11e6cd06d246c5c327df4305d6861533f225f","tier":"starter_2024"},"version":1082,"publish_version":1082,"role":null,"authored":false,"authors":[{"type":"individual","id":"d34a26730351ade8","login":"renatoppl","name":"Renato Paes Leme","avatar_url":"https://avatars-next.observableusercontent.com/c22df5c297aca2188beaa018c7c11e6cd06d246c5c327df4305d6861533f225f","tier":null}],"members":[{"type":"individual","id":"d34a26730351ade8","login":"renatoppl","name":"Renato Paes Leme","avatar_url":"https://avatars-next.observableusercontent.com/c22df5c297aca2188beaa018c7c11e6cd06d246c5c327df4305d6861533f225f","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"subscribed":null,"forks":6,"collection":null,"type":"notebook","body":{"title":"Laplace Equation","cells":[{"id":0,"value":"md`# Laplace Equation`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":805,"value":"md`[Laplace's equation](https://en.wikipedia.org/wiki/Laplace%27s_equation) appears everywhere in physics (fluid dynamics, gravitation, eletrostatics, heat conduction, ...). Here we will see a different application in Complex Analysis: the [Riemann conformal mapping theorem](https://en.wikipedia.org/wiki/Riemann_mapping_theorem). We will get to that in a minute. First, let me say what the equation is. We will be concerned with functions defined on a region of the plane, i.e.\n${tex.block`u : U \\subseteq \\mathbb{R}^2 \\rightarrow \\mathbb{R}`}\nWe will say that ${tex`u`} is _harmonic_ (i.e. a solution to Laplace's equation), when:\n${tex.block`\\Delta u(z) = 0 \\quad \\text{for } z \\in U`}\nwhere ${tex`\\Delta`} is the Laplacian operator:\n${tex.block`\\Delta u := \\frac{\\partial^2 u}{\\partial x^2} + \\frac{\\partial^2 u}{\\partial y^2}`}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":818,"value":"$30","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":842,"value":"$31","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":861,"value":"md`\n### Riemann Conformal Mapping Theorem\nIn the rest of this notebook, we will show how to (numerically) solve the Laplace equation and apply it to get a construction in the [Riemann conformal mapping theorem](https://en.wikipedia.org/wiki/Riemann_mapping_theorem). The Conformal Mapping Theorem says that given a simple compact region ${tex`U \\subset \\mathbb{C}`} there is always an holomorphic dipheomorphism from ${tex`U`} to the unit ball. In other words, there is a mapping:\n${tex.block`f:U \\rightarrow B := \\{z \\in \\mathbb{C}; \\vert z \\vert \\leq 1 \\}`}\nsuch that ${tex`f`} and ${tex`f^{-1}`} satisfy the equations of Cauchy-Riemann. The name _conformal_ refers to the fact that holomorphic functions preserve angles, i.e., given vectors ${tex`v_1, v_2 \\in \\mathbb{C}`} the angle between them is the same as the angle between ${tex`f'(z) v_1`} and ${tex`f'(z) v_2`}.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":626,"value":"$32","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":689,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":646,"value":"md`Now that we computed ${tex`u`}, we can compute ${tex`v`} using the Cauchy-Riemann equations. Since we know the gradient of ${tex`v`} since:\n${tex.block`\\nabla v(z) = \\begin{bmatrix} 0 & -1 \\\\ 1 & 0 \\end{bmatrix} \\nabla u(z) `}\nWe can just find ${tex`v`} by integrating its gradient. The solution is the function below. Now, we can construct an holomorphic function ${tex`g : \\mathbb{C} \\rightarrow \\mathbb{C}`} by setting:\n${tex.block` g(z) = u(z) + i v(z) `}\nThis is not yet the function we are looking for, since we want a bijection to the disk.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":561,"value":"{\n const height = 320;\n const width = 640;\n const context = DOM.context2d(width, height, 1);\n \n function update() {\n context.clearRect(0, 0, width, height);\n var points = viewof LaplaceData.value.points;\n var interior = viewof LaplaceData.value.interior;\n context.strokeStyle = \"red\";\n var n = points.length\n context.beginPath();\n context.moveTo(points[n-1][0], points[n-1][1]);\n for (var i = 0; i < n; ++i) {\n context.lineTo(points[i][0], points[i][1]);\n context.stroke();\n }\n var vBar = ComplexComponent(viewof LaplaceData.value.uBar, \n viewof LaplaceData.value.closure)\n viewof LaplaceData.setValue(\"vBar\", vBar);\n FillPoints(interior, vBar, context);\n viewof LaplaceData.startEvent(\"e2\");\n viewof LaplaceData.startEvent(\"e3\");\n }\n\n viewof LaplaceData.setEvent(\"e1\", update);\n return context.canvas;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":651,"value":"md`Constructing the bijection is easy, we just set:\n${tex.block`f(z) = (z-z_0) \\cdot \\exp(g(z))`}\nIt is simple to see it maps to the disk since:\n${tex.block`\\vert f(z) \\vert = \\vert z-z_0\\vert \\cdot \\exp(u(z)) = 1`}\non the boundary. By the [maximum principle](https://en.wikipedia.org/wiki/Maximum_principle) it should map to the disk. We still need to argue that it is a bijection and the inverse is holomorphic. See [here](http://www-users.math.umn.edu/~olver/ln_/cml.pdf) for the details. Below is a mapping from the points of ${tex`U`} to the points in the disk.\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":568,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":668,"value":"md`The next picture is a mapping from a grid on the shape to a grid in the circle.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":598,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":881,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined"},{"id":929,"value":"viewof button1 = {\n const button = html`