2d:null 2e:["$","title",null,{"children":"Numerical solution of 1D, linear BVPs | Mark McClure | Observable"}] 33:T61c,md`### Difference quotients In a first semester calculus course, you learn about the *difference quotient* that's used to define the derivative: ${tex.block`u'(x) \approx \frac{u(x+h)-u(x)}{h}.`} One fundamental technique to find approximate solutions to first order differential equations is to replace any occurence of ${tex`u'(x)`} with its difference quotient approximation ${tex`(u(x+h)-u(x))/h`}. If we know ${tex`u(x)`} (from, say, a boundary or initial condition) we can then solve for ${tex`u(x+h)`}, then use that to solve for ${tex`u(x+2h)`} and on down the line. This is called [Euler's method](https://en.wikipedia.org/wiki/Euler_method). In order to extend this idea to second order equations, we need a difference quotient for second derivatives. To find such a thing, we can use Taylor's theorem to expand ${tex`u(x+h)`} and ${tex`u(x-h)`} into power series: ${tex.block`u(x+h) = u(x) + u'(x)h + \frac{1}{2}u''(x)h^2 + \frac{1}{6}u'''(x)h^3 + O(h^4),`} ${tex.block`u(x-h) = u(x) - u'(x)h + \frac{1}{2}u''(x)h^2 - \frac{1}{6}u'''(x)h^3 + O(h^4).`} The big ${tex`O`} term indicates that the error of the order ${tex`h^4`}; that is, all subsequent terms involve powers of ${tex`h^n`}, where ${tex`n>4`}. Since ${tex`h`} is generally much less than one in our application, these terms will be much smaller than the previous terms. Now, we can add these two equations to get ${tex.block`u(x+h)+u(x-h) = 2u(x) + u''(x)h^2 + O(h^4)`} and then solve for ${tex`u''(x)`} to get ${tex.block`u''(x) = \frac{u(x+h) - 2u(x) + u(x-h)}{h^2} + O(h^2).`} `34:Tc64,partition_pic = { let w = 700; let h = 100; let xmin = -0.02; let xmax = 1.02; let ymin = -0.1; let ymax = 0.1; let pad = 15; let svg = d3 .create('svg') .attr('width', w) .attr('height', h); // .style('border', 'solid 1px black'); 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 pts_to_path = d3 .line() .x(function(d) { return x_scale(d[0]); }) .y(function(d) { return y_scale(d[1]); }); let n = 8; svg .append('line') .attr('x1', x_scale(0)) .attr('x2', x_scale(1)) .attr('y1', y_scale(0)) .attr('y2', y_scale(0)) .attr('stroke', 'black'); svg .selectAll('circle.interior') .data(d3.range(1 / n, 1, 1 / n)) .join('circle') .attr('class', 'interior') .attr('cx', d => x_scale(d)) .attr('cy', y_scale(0)) .attr('r', 3) .attr('fill', 'black'); svg .selectAll('circle.endpoint') .data([0, 1]) .join('circle') .attr('class', 'endpoint') .attr('cx', d => x_scale(d)) .attr('cy', y_scale(0)) .attr('r', 4) .attr('fill', '#ccc') .attr('stroke', 'black'); let xi = 0.375; let tex_group = svg.append('g'); tex_group .append('g') .attr('transform', `translate(${x_scale(xi - 0.01)} ${y_scale(-0.04)})`) .append(() => MathJax.tex2svg(String.raw`x_i`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(xi - 0.15)} ${y_scale(-0.04)})`) .append(() => MathJax.tex2svg(String.raw`x_{i-1}`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(xi + 0.1)} ${y_scale(-0.04)})`) .append(() => MathJax.tex2svg(String.raw`x_{i+1}`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(-0.04)} ${y_scale(-0.04)})`) .append(() => MathJax.tex2svg(String.raw`x_0=a`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(0.95)} ${y_scale(-0.04)})`) .append(() => MathJax.tex2svg(String.raw`x_n=b`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(xi - 0.05)} ${y_scale(0.085)})`) .append(() => MathJax.tex2svg(String.raw`u_i \approx u(x_i)`).querySelector("svg") ); tex_group .append('g') .attr('transform', `translate(${x_scale(xi - 0.15)} ${y_scale(0.07)})`) .append(() => MathJax.tex2svg(String.raw`u_{i-1}`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(xi + 0.1)} ${y_scale(0.07)})`) .append(() => MathJax.tex2svg(String.raw`u_{i+1}`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(-0.04)} ${y_scale(0.085)})`) .append(() => MathJax.tex2svg(String.raw`u_0=u(x_0)`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(0.9)} ${y_scale(0.085)})`) .append(() => MathJax.tex2svg(String.raw`u_n=u(x_n)`).querySelector("svg")); return svg.node(); }35:Taec,simple_partition_pic = { let w = 700; let h = 100; let xmin = -0.02; let xmax = 1.02; let ymin = -0.1; let ymax = 0.1; let pad = 15; let svg = d3 .create('svg') .attr('width', w) .attr('height', h); // .style('border', 'solid 1px black'); 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 pts_to_path = d3 .line() .x(function(d) { return x_scale(d[0]); }) .y(function(d) { return y_scale(d[1]); }); let n = 10; svg .append('line') .attr('x1', x_scale(0)) .attr('x2', x_scale(1)) .attr('y1', y_scale(0)) .attr('y2', y_scale(0)) .attr('stroke', 'black'); svg .selectAll('circle.interior') .data(d3.range(1 / n, 1, 1 / n)) .join('circle') .attr('class', 'interior') .attr('cx', d => x_scale(d)) .attr('cy', y_scale(0)) .attr('r', 3) .attr('fill', 'black'); svg .selectAll('circle.endpoint') .data([0, 1]) .join('circle') .attr('class', 'endpoint') .attr('cx', d => x_scale(d)) .attr('cy', y_scale(0)) .attr('r', 4) .attr('fill', '#ccc') .attr('stroke', 'black'); let xi = 0.4; let tex_group = svg.append('g'); tex_group .append('g') .attr('transform', `translate(${x_scale(xi - 0.045)} ${y_scale(-0.04)})`) .append(() => MathJax.tex2svg(String.raw`x_i=i/10`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(-0.04)} ${y_scale(-0.04)})`) .append(() => MathJax.tex2svg(String.raw`x_0=0`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(0.95)} ${y_scale(-0.04)})`) .append(() => MathJax.tex2svg(String.raw`x_{10}=1`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(xi - 0.04)} ${y_scale(0.085)})`) .append(() => MathJax.tex2svg(String.raw`u_i = \, ?`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(xi - 0.12)} ${y_scale(0.07)})`) .append(() => MathJax.tex2svg(String.raw`u_{i-1}`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(xi + 0.07)} ${y_scale(0.07)})`) .append(() => MathJax.tex2svg(String.raw`u_{i+1}`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(-0.04)} ${y_scale(0.085)})`) .append(() => MathJax.tex2svg(String.raw`u_0=0`).querySelector("svg")); tex_group .append('g') .attr('transform', `translate(${x_scale(0.95)} ${y_scale(0.085)})`) .append(() => MathJax.tex2svg(String.raw`u_{10}=0`).querySelector("svg")); return svg.node(); }36:T4dd,md`Again, we want values for the ${tex`u_i`}s that should be good approximations each ${tex`u(x_i)`}. We already know ${tex`u_0`} and ${tex`u_{10}`} from the boundary conditions. We can obtain equations relating the others by setting ${tex`h=0.1`} and ${tex`x=x_i`} in our [difference quotient for the second derivative (2)](https://observablehq.com/d/02363a70f5fffe82#diff_quotient_2). Doing so, we get ${tex.block`u''(x_i) \approx \frac{u(x_{i+1}) - 2u(x_i) + u(x_{i-1})}{0.1^2}`} or ${tex.block`100(u_{i+1} - 2u_i + u_{i-1}) = 2.`} Letting ${tex`i=1,2,\ldots,9`}, this yields nine equations in the relating ${tex`u_0,u_1,\ldots,u_{10}`}. Note that the equation for ${tex`u_1`} involves ${tex`u_0`} and the equation for ${tex`u_9`} involves ${tex`u_{10}`}, but ${tex`u_0`} and ${tex`u_{10}`} are known to have the value ${tex`0`}. Thus, we have nine equations in nine unkowns that we should be able to solve. #### Concrete formulation and solution There are loads of software packages for solving systems numerically. As input, they often require a matrix ${tex`A`} and a vector ${tex`\vec{b}`} so that the system is expressed as ${tex.block`A\vec{u}=\vec{b}.`} Thus, it will be convenient to express the system in that form, i.e. like so: `37:T645,{ let height = 0.6 * width; let container = d3 .create('div') .style('position', 'relative') .style('width', `${width}px`) .style('height', `${height}px`); let graph = container .append('div') .style('position', 'absolute') .style('width', `${width}px`) .style('height', `${height}px`); let target = graph.node(); let h = 1 / n; //let sol = [0, simple_solution.data.map(d => d[0]), 0].flat(); let pts = d3.zip(d3.range(0, 1 + h, h), simple_solution); console.log(['pts are', pts]); let ymin = d3.min(simple_solution); let ymax = d3.max(simple_solution); functionPlot({ target: target, width: width, height: height, xAxis: { domain: [-0.01, 1.01] }, yAxis: { domain: [-0.255, 0.01] }, disableZoom: true, data: [ { fn: 'x^2-x' }, { points: pts, fnType: 'points', graphType: 'polyline' }, { points: pts, fnType: 'points', graphType: 'scatter' } ] }); graph.selectAll('.origin').style('opacity', 0); graph .selectAll('circle') .attr('r', 4) .style('fill', 'black') .style('stroke', 'black'); graph .selectAll('.line-0') .style('stroke-width', '4px') .style('stroke', 'red'); graph .selectAll('.line-1') .style('stroke-width', '2px') .style('stroke', 'green'); container .append('div') .style('position', 'absolute') .style('width', `${width}px`) .style('height', `${height}px`) .style('top', '0px') .style('left', '0px'); return container.node(); }38:Tefb,md`## Extending the functionality The preceding simple example is fun and all but can clearly be generalized in a number of ways: 1. We'd like to include first and zeroth order terms, ${tex`p\,u'`} and ${tex`q\,u`}. 2. We'd like ${tex`f`} to be a function; while we're at it, we might as well allow ${tex`p`} and ${tex`q`} to be functions as well. 3. We'd like non-homogeneous (i.e., arbitary) boundary conditions. ### Adding first and zeroth order terms Just as the second order term was represented via a matrix, we'll need to represent the first and zeroth order terms with matrices. The key to the first order term is [equation (3)](#diff_quotient_1) that expresses an improved difference quotient for the first derivative. When expressed as a matrix, that yields ${tex.block` \frac{1}{2h} \begin{pmatrix} 0 & 1 & 0 & 0 & \cdots \\ -1 & 0 & 1 & 0 & \cdots \\ 0 & -1 & 0 & 1 & \cdots \\ 0 & 0 & -1 & 0 & \cdots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix}.`} The zeroth order term is easy, since the matrix representing ${tex`u`} itself is just the identity matrix. ### Allowing ${tex`f`} to be a function This part is easy. Recall that we discretize the interval ${tex`[a,b]`} into sub-intervals with endpoints ${tex`x_i=a+i\times h`}. When setting up the column vector representing the right side of our system, we simply choose the ${tex`i^{\text{th}}`} entry to be ${tex`f(x_i)`}. A similar idea works if ${tex`p`} and ${tex`q`} are functions. When we multiply the first or zeroth order matrix terms by ${tex`p`} or ${tex`q`}, we simply multiply the ${tex`i^{\text{th}}`} row by ${tex`p(x_i)`} or ${tex`q(x_i)`} ### Assembling the system If we account for what we've got so far, we should end up with a system that looks like so: ${tex.block` \left[\frac{1}{h^2} \begin{pmatrix} -2 & 1 & 0 & \cdots \\ 1 & -2 & 1 & \cdots \\ 0 & 1 & -2 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix} + \frac{1}{2h} \begin{pmatrix} 0 & p(x_1) & 0 & \cdots \\ -p(x_2) & 0 & p(x_2) & \cdots \\ 0 & -p(x_3) & 0 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix}\right. \\ + \left.\begin{pmatrix} q(x_1) & 0 & 0 & \cdots \\ 0 & q(x_2) & 0 & \cdots \\ 0 & 0 & q(x_3) & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{pmatrix}\right] \begin{pmatrix} u_1 \\ u_2 \\ u_3 \\ \vdots \end{pmatrix} = \begin{pmatrix} f(x_1) \\ f(x_2) \\ f(x_3) \\ \vdots \end{pmatrix} `} We can add those three matrices together to generate one big system of the form ${tex`M\vec{u}=\vec{b}`} and then pass ${tex`M`} and ${tex`\vec{b}`} to a linear solver to find ${tex`\vec{u}`}. ### Non-homogeneous boundary conditions At this point, we have a strategy to numerically solve equations of the [form (1)](#eq1), as long as ${tex`u_a=u_b=0`}. If ${tex`u_a`} and/or ${tex`u_b`} are non-zero, the strategy will be to translate the equation to an equivalent equation with homogeneous boundary conditions. We then solve the new equation and finally translate back. To be clear, suppose that ${tex`w`} solves ${tex.block`w'' + p\,w' + q\,w = F,`} where ${tex`F(x) = f(x) - \left[A\,p(x) + (A\,x+B)q(x)\right]`}. Then, ${tex.block` \begin{aligned} & (w+Ax+B)'' + p\,(w+Ax+B)' + q\,(w+Ax+B) = w'' + p\,w' + q\,w + A\,p(x) + (A\,x+B)q(x) \\ &= F(x) + A\,p(x) + (A\,x+B)q(x) = f(x) - \left[A\,p(x) + (A\,x+B)q(x)\right] + A\,p(x) + (A\,x+B)q(x) = f. \end{aligned} `} Thus, to solve the equation with arbitrary values of ${tex`u_a`} and ${tex`u_b`}, we simply choose ${tex`A`} and ${tex`B`} so that ${tex.block` \begin{aligned} Aa+B &= u_a \\ Ab+B &= u_b. \end{aligned} `} We then solve the system ${tex.block`w'' + p(x)\,w' + q(x)\,w = f(x) - \left[A\,p(x) + (A\,x+B)q(x)\right].`} Finally, we set ${tex`u(x) = w(x)+Ax+b`}. That's *all* there is to it!! 😕 `39:T6cc,// A function to plot the solutions that arise from this technique function plot(solution, opts = {}) { let { actual = false } = opts; let height = 0.6 * width; let container = d3 .create('div') .style('position', 'relative') .style('width', `${width}px`) .style('height', `${height}px`); let graph = container .append('div') .style('position', 'absolute') .style('width', `${width}px`) .style('height', `${height}px`) .style('background-color', 'white'); let target = graph.node(); let functionPlotData = [ { points: solution, fnType: 'points', graphType: 'polyline' }, { points: solution, fnType: 'points', graphType: 'scatter' } ]; if (actual) { functionPlotData.unshift({ fn: actual }); } console.log(functionPlotData); let ymin = d3.min(solution.map(pt => pt[1])); let ymax = d3.max(solution.map(pt => pt[1])); functionPlot({ target: target, width: width, height: height, xAxis: { domain: [solution[0][0], solution[solution.length - 1][0]] }, yAxis: { domain: [ymin, ymax] }, disableZoom: true, data: functionPlotData }); graph.selectAll('.origin').style('opacity', 0); graph .selectAll('circle') .attr('r', 4) .style('fill', 'black') .style('stroke', 'black'); graph .selectAll('.line-0') .attr('stroke-width', 4) .attr('stroke', 'red'); graph .selectAll('.line-1') .attr('stroke-width', 2) .attr('stroke', 'green'); container .append('div') .style('position', 'absolute') .style('width', `${width}px`) .style('height', `${height}px`) .style('top', '0px') .style('left', '0px'); return container.node(); }3a:[] 31:{"id":"02363a70f5fffe82","slug":"numerical-solution-of-1d-linear-bvps","title":"Numerical solution of 1D, linear BVPs","thumbnail":"6a31a0f8fec2ef5ec273ade279be6b39e4a2aff465f6ec08ff6ea50ea2ea24c3","access_level":"public","create_time":"$D2021-01-09T23:54:28.062Z","update_time":"$D2024-04-23T16:09:55.154Z","view_time":null,"first_public_time":"$D2021-03-17T18:58:13.936Z","has_unpublished_edits":null,"fork":false,"stars":2,"starred":false,"owner":{"type":"team","id":"d771c8ae812d5130","login":"mcmcclur","name":"Mark McClure","avatar_url":"https://avatars-next.observableusercontent.com/e615e819ce83d43faef607e48c53c5ff63943aa0ba9f5848499bd02fef7396f4","tier":"starter_2024"},"version":1368,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"c3301e515dd08a41","login":"mcmcclur","name":"Mark McClure","avatar_url":"https://avatars-next.observableusercontent.com/e615e819ce83d43faef607e48c53c5ff63943aa0ba9f5848499bd02fef7396f4","tier":null}],"members":[{"type":"individual","id":"c3301e515dd08a41","login":"mcmcclur","name":"Mark McClure","avatar_url":"https://avatars-next.observableusercontent.com/e615e819ce83d43faef607e48c53c5ff63943aa0ba9f5848499bd02fef7396f4","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":2,"collection":{"id":"524f747a860d90c5","slug":"teaching-pdes","title":"Teaching PDEs","owner":{"id":"d771c8ae812d5130","login":"mcmcclur","name":"Mark McClure"}},"type":"notebook","body":{"title":"Numerical solution of 1D, linear BVPs","cells":[{"id":0,"value":"md`# Numerical solution of 1D, linear BVPs\n## with the finite difference method`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":151,"value":"md`The objective here is to describe a technique called *the finite difference method* to solve linear boundary value problems of the form`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":281,"value":"eq1 = tex.block`u'' + p\\,u'+q\\,u = f(x), \\: u(a)=u_a, \\: u(b) = u_b. \\tag{1}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":286,"value":"md`These types of equations are of tremendous importance in applied mathematics. They arise, for example, via separation of variables in partial differential equations; as such, they are related to heat flow and wave phenomena.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":160,"value":"$33","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":181,"value":"md`### Notation for discretization\n\nUltimately, numerical approximations to boundary value problems of the [form (1)](#eq1) are obtained by *discretization*; that is, we replace the continuous differential equation over the interval ${tex`[a,b]`} with a collection of linear algebraic equations indexed by a discrete set of points chosen from that interval. More specifically, given a positive integer ${tex`n`}, we break the interval ${tex`[a,b]`} up evenly into ${tex`n`} sub-intervals to yield a partition that looks like so:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":185,"value":"$34","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":297,"value":"md`Since we break the interval up evenly, each of these sub-intervals has length ${tex`h=(b-a)/n`}. We denote the endpoints of the sub-intervals by ${tex`x_i`}; thus, each ${tex`x_i`} may be computed as\n${tex.block`x_i= a+i\\times h.`}\nThe objective of the finite difference method is to compute estimates for the value of the solution ${tex`u`} at each point ${tex`x_i`}. We'll denote this estimate by ${tex`u_i`}; thus,\n${tex.block`u_i \\approx u(x_i)`}\nfor each ${tex`i=1,\\ldots,n-1`}. Of course, we know the value of ${tex`u`} exactly at the endpoints from the boundary conditions.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":503,"value":"diff_quotient_2 = md`${tex.block`u''(x) = \\frac{u(x+h) - 2u(x) + u(x-h)}{h^2} + O(h^2). \\tag{2}`}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":506,"value":"md`\n\nOur error term is of the order ${tex`h^2`}, which is even better than the error inherent in the basic difference quotient for ${tex`u'`}. It's worth mentioning that we can use this same technique to find an improved difference quotient for the first derivative. If we subtract the second equation from the first, we get\n\n${tex.block`u(x+h)-u(x-h) = 2u'(x)h + O(h^3).`}\n\nWe can now solve for ${tex`u'`} to get`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":900,"value":"diff_quotient_1 = tex.block`u'(x) = \\frac{u(x+h) - u(x-h)}{2h} + O(h^2). \\tag{3}`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":315,"value":"simple_example_description = md`### Simple example\n\nLet's describe the technique with one of the simplest possible examples,\n${tex.block`u''=2,\\;u(0)=u(1)=0.`}\nNote that this equation has a simple closed form solution, namely ${tex`u(x)=x^2-x`}; thus, it will be easy to check the quality of our technique against this formula.\n\n#### Abstract formulation\nTo get started, we partion the interval from zero to one into 10 sub-intervals to yield a picture that looks like:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":470,"value":"$35","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":498,"value":"$36","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":540,"value":"{\n let n = 10;\n let h = 1 / n;\n let L = mlm.Matrix.zeros(n - 1, n - 1);\n for (let i = 0; i < n - 1; i++) {\n L.set(i, i, -2);\n if (i < n - 2) {\n L.set(i + 1, i, 1);\n L.set(i, i + 1, 1);\n }\n }\n L = L.mul((1 / h) ** 2);\n\n return tex`${matrix_to_TeX(L)}\n\\begin{pmatrix}\n u_1 \\\\ u_2 \\\\ u_3 \\\\ u_4 \\\\ u_5 \\\\ u_6 \\\\ u_7 \\\\ u_8 \\\\ u_9\n\\end{pmatrix} =\n\\begin{pmatrix}\n 2 \\\\ 2 \\\\ 2 \\\\ 2 \\\\ 2 \\\\ 2 \\\\ 2 \\\\ 2 \\\\ 2\n\\end{pmatrix}.\n`;\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":862,"value":"md`Now, the following bit of code sets up this matrix and solves the system using a Javascript library called [ml-matrix](https://mljs.github.io/matrix/).`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":555,"value":"simple_solution = {\n // Step size h defined in terms of parition size n\n // n is defined by a slider below.\n let h = 1 / n;\n\n // Initialize a matrix of zeros\n let L = mlm.Matrix.zeros(n - 1, n - 1);\n\n // Reset values on the diagonal and near diagonals\n for (let i = 0; i < n - 1; i++) {\n L.set(i, i, -2);\n if (i < n - 2) {\n L.set(i + 1, i, 1);\n L.set(i, i + 1, 1);\n }\n }\n // Multiply by 100\n L = L.mul((1 / h) ** 2);\n\n // Set up a the right hand side, a vector with value 2 in each entry.\n let B = mlm.Matrix.columnVector(d3.range(n - 1).map(() => 2));\n\n // Solve the system\n let solution = mlm.solve(L, B).getColumn(0);\n solution = [0].concat(solution.concat([0]));\n return solution;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":593,"value":"md`The solution is just a list of numbers representing approximations to each ${tex`u(x_i)`}. We can plot this computed solution, together with the graph of the analytic solution ${tex`y=x^2-x`}:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":654,"value":"viewof n = Range([2, 100], { label: tex`n`, step: 1, value: 6 })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":641,"value":"$37","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":868,"value":"md`Not bad! It should be noted that it's rare for the computed solution to be this close to exact at the actual points; this only happens when the actual solution is a polynomial of degree two or less.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":736,"value":"$38","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":976,"value":"md`## Examples\n\nThe code below defines a \\`solveBVP\\` function that implements exactly the strategy defined above. The signature is\n\n solveBVP(a, b, ua, ub, p, q, f, n)\n\nI guess the argments are self-explanatory, if you've gotten this far. We've also defined a \\`plot\\` function that's built to plot exactly the type of solution that \\`solveBVP\\` produces.\n\nAs an example, here the solution of\n${tex.block`u'' - x\\,u' + x^2 u = 2 x^5-x^4-10 x^3+2 x^2+16 x-2, \\: u(-1)=2, \\: u(2)=4.`}\nNote that this equation was built so that the solution would be \n${tex.block`u(x)=2x^3-x^2-4x.`}\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1035,"value":"plot(\n solveBVP(\n -1,\n 2,\n 1,\n 4,\n (x) => -x,\n (x) => x ** 2,\n (x) => -2 + 16 * x + 2 * x ** 2 - 10 * x ** 3 - x ** 4 + 2 * x ** 5,\n 15\n ),\n {\n actual: \"2*x^3-x^2-4*x\"\n }\n)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1336,"value":"//Signature:\n// solveBVP(a, b, ua, ub, p, q, f, n)\n// where we're trying to solve\n// u'' + p*u' + q*u = f(x)\n// with u(a) = u_a and u(b) = u_b\n\nsample = plot(\n solveBVP(\n 0,\n 13,\n 1,\n 0,\n (x) => 0,\n (x) => x,\n (x) => x,\n 400\n )\n)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1099,"value":"md`## Eigensystems\n\nThese same tools can be used to solve eigenfunction problems as well. Consider the matrix ${tex`L`} that represents the second derivative over an interval:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1227,"value":"tex.block`\nL=\\frac{1}{h^2}\n\\begin{pmatrix}\n -2 & 1 & 0 & 0 & \\cdots \\\\\n 1 & -2 & 1 & 0 & \\cdots \\\\\n 0 & 1 & -2 & 1 & \\cdots \\\\\n 0 & 0 & 1 & -2 & \\cdots \\\\\n \\vdots & \\vdots & \\vdots & \\vdots & \\ddots \n\\end{pmatrix}.`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1234,"value":"md`Now suppose that \n${tex.block`\\vec{u}=\\langle u_1,u_2,\\ldots,u_{n-1} \\rangle`}\nis an eigenvector for ${tex`L`}. If we plot the points\n${tex.block`\\{(x_1,u_1), (x_2,u_2),\\ldots,(x_{n-1},u_{n-1})\\},`}\nWe should get a decent approximation to the graph of an eigenfunction of ${tex`L`}.\n\nHere's an illustration that examines ${tex`u''=\\lambda u`} over the unit interval, so that ${tex`h=1/N`}, where ${tex`N`} is the partition size:\n`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1240,"value":"viewof N = Range([10, 100], {\n label: tex`N=\\text{partition size}`,\n step: 1,\n value: 50\n})","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1137,"value":"viewof m = Range([1, N / 3], { label: 'mode', step: 1, value: 1 })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1127,"value":"{\n let height = 0.6 * width;\n let container = d3\n .create('div')\n .style('position', 'relative')\n .style('width', `${width}px`)\n .style('height', `${height}px`);\n let graph = container\n .append('div')\n .style('position', 'absolute')\n .style('width', `${width}px`)\n .style('height', `${height}px`);\n\n let target = graph.node();\n\n let ymin = d3.min(pts.map(pt => pt[1]));\n let ymax = d3.max(pts.map(pt => pt[1]));\n functionPlot({\n target: target,\n width: width,\n height: height,\n xAxis: { domain: [pts[0][0], pts[pts.length - 1][0]] },\n yAxis: { domain: [ymin, ymax] },\n disableZoom: true,\n data: [\n {\n points: pts,\n fnType: 'points',\n graphType: 'polyline'\n },\n {\n points: pts,\n fnType: 'points',\n graphType: 'scatter'\n }\n ]\n });\n return container.node();\n}","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1251,"value":"md`Well the result sure looks like functions of the forum ${tex`\\sin(n\\pi x)`}!\n\nIt's also worth examining the eigenvalues. Here are the last few that are smallest in magnitude:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1248,"value":"ed.realEigenvalues.slice(-5)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1301,"value":"md`And here are minus the squares of ${tex`n\\pi`}:`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1275,"value":"d3\n .range(1, 6)\n .map(n => -1 * (n * Math.PI) ** 2)\n .reverse()","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1107,"value":"ed = new mlm.EigenvalueDecomposition(L)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1102,"value":"L = term2(0, 1, N)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1297,"value":"pts = [0]\n .concat(ed.eigenvectorMatrix.getColumn(N - (m + 1)))\n .concat([0])\n .map((y, i) => [(i + 1) / N, y])","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":986,"value":"md`## Code`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":998,"value":"md`### Computation`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":825,"value":"// The main function\nfunction solveBVP(a, b, ua, ub, p, q, f, n) {\n let h = (b - a) / n;\n\n let L = term2(a, b, n);\n let D = term1(a, b, n, p);\n let C = term0(a, b, n, q);\n let M = L.addM(D).addM(C);\n\n let A = (ub - ua) / (b - a);\n let B = (b * ua - a * ub) / (b - a);\n\n let F;\n if (typeof f == 'number') {\n F = x => f + (A * x + B);\n } else if (typeof f == 'function') {\n F = x => f(x) - q(x) * (A * x + B) - p(x) * A;\n }\n\n let RHS = mlm.Matrix.columnVector(d3.range(a + h, b - h / 2, h).map(F));\n let solution = mlm.solve(M, RHS).getColumn(0);\n solution = [0].concat(solution.concat([0]));\n solution = solution.map((u, i) => u + (A * (a + i * h) + B));\n solution = d3.zip(d3.range(a, b + h / 2, h), solution);\n return solution;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":746,"value":"// The second order term, u''\nfunction term2(a, b, n) {\n let h = (b - a) / n;\n let L = mlm.Matrix.zeros(n - 1, n - 1);\n for (let i = 0; i < n - 1; i++) {\n L.set(i, i, -2);\n if (i < n - 2) {\n L.set(i + 1, i, 1);\n L.set(i, i + 1, 1);\n }\n }\n L = L.mul((1 / h) ** 2);\n return L;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":755,"value":"// The first order term, p(x)*u'\nfunction term1(a, b, n, p) {\n let h = (b - a) / n;\n let n_over_2h = n / (2 * (b - a));\n let M = mlm.Matrix.zeros(n - 1, n - 1);\n if (typeof p == 'number') {\n for (let i = 0; i < n - 2; i++) {\n M.set(i + 1, i, -n_over_2h * p);\n M.set(i, i + 1, n_over_2h * p);\n }\n } else if (typeof p == 'function') {\n for (let i = 0; i < n - 2; i++) {\n M.set(i + 1, i, -n_over_2h * p(a + (i + 2) * h));\n M.set(i, i + 1, n_over_2h * p(a + (i + 1) * h));\n }\n }\n return M;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":768,"value":"// The zeroth order term, q(x)*u\nfunction term0(a, b, n, q) {\n let h = (b - a) / n;\n let M = mlm.Matrix.zeros(n - 1, n - 1);\n if (typeof q == 'number') {\n for (let i = 0; i < n - 1; i++) {\n M.set(i, i, q);\n }\n } else if (typeof q == 'function') {\n for (let i = 0; i < n - 1; i++) {\n M.set(i, i, q(a + (i + 1) * h));\n }\n }\n return M;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":733,"value":"md`### Utilities`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":990,"value":"$39","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":362,"value":"function matrix_to_TeX(M) {\n let ff = d3.format(\"0.2~f\");\n function format(x) {\n let formatted = ff(x);\n if (formatted == \"NaN\") {\n return '\\\\vdots';\n } else {\n return formatted;\n }\n }\n function makeRow(R) {\n let result = '';\n let cut;\n if (M.columns < 10) {\n cut = -1;\n } else {\n cut = 9;\n }\n R.slice(0, cut).forEach(function(x) {\n result = result + format(x) + ' & ';\n });\n if (cut == 9) {\n result = result + '\\\\cdots &';\n }\n result = result + format(R.slice(-1)[0]) + '\\\\\\\\';\n return result;\n }\n let MJSON = M.toJSON();\n let tex_string = '\\\\begin{pmatrix}';\n MJSON.slice(0, 9).forEach(function(r) {\n tex_string = tex_string + makeRow(r);\n });\n if (M.columns > 10) {\n tex_string = tex_string + makeRow(d3.range(M.columns).map(x => '.'));\n }\n if (M.columns >= 10) {\n tex_string = tex_string + makeRow(MJSON.slice(-1)[0]);\n }\n tex_string = tex_string + ' \\\\end{pmatrix}';\n return tex_string;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":730,"value":"md`## Imports`","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1325,"value":"sin = Math.sin","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1327,"value":"cos = Math.cos","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1329,"value":"exp = Math.exp","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1331,"value":"pi = Math.PI","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":1333,"value":"e = Math.E","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":652,"value":"import { Range } from \"@observablehq/inputs\"","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":724,"value":"MathJax = require('https://cdn.jsdelivr.net/npm/mathjax@3/es5/tex-svg.js').catch(\n () => window['MathJax']\n)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":118,"value":"functionPlot = require(\"function-plot@1/dist/function-plot\")","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":2,"value":"mlm = require('ml-matrix')","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"},{"id":48,"value":"d3 = require('d3@6')","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W3a","stdlib":"1","titling":"auto"}} 32:"$31:body"