Solving 2nd order differential equations in JavaScript using Finite Difference Method
The basic flow to solve the differential equations numerically is
- Find differential equations to solve
- Apply approximation such as centered difference method
- Set boundary conditions
- Solve system of linear equations (Matrix problem)
I was very fascinated in numerical method. When I first learned it, it was not that simple to understand how it works. I am hoping this article would bring down the barrier of numerical method of differential equations while I improve my JavaScript skills.
I also wrote the article doing the same but in two dimension.
Here is the 2nd order differential equation
The generic form of 2nd order differential equation is
Now there is a way to make this equation in approximated form. The centered difference method is the one, and at point xₙ the method yields the following form,
Please see references at the end of my post for more details!
The basic idea is to use the Taylor series expansion to approximate the equation, up to the second order in this case. I was always skeptical using approximation but I found that the approximation is quite good.
Laplace Equation
Now this is the generic form. Let’s make it simpler! The simplest form of the equation is called Laplace equation where p(x), q(x) and r(x) are all 0. So it will look
we will ignore the big O term where the approximation comes from, so the right hand side is 0. For the Laplace equation, this can be further simplified to
Poisson Equation
For the Poisson equation, r(x) is constant. So, multiplying rΔx² gives
Only right hand side is the difference between the two.
Let us start with Laplace equation in one dimension because it is the simplest of the simplest, followed up with Poisson equation in one dimension.
Building Matrix — Laplace Equation 1D
Set Boundary Condition
Since this is the 2nd order differential equations, there are numerous solutions exist unless boundary conditions are set. There are types of boundary conditions such as Dirichlet and Neumann boundary conditions. Let’s consider Dirichlet boundary condition because it is simpler. Let
which means at grid point 0, u is set to value a, and at grid point n, u is set to value b.
Build Matrix
hence, starting from u₀ = a, applying the above equation, there are n+1 equations to solve. Imagine there are grid points left to right, starting from point 0 to n. At point 0, it is u₀=a, at point 1, we get u₀–2u₁+u₂=0, and so on. Therefore
This can be written in the matrix form,
Or simply write in the Matrix form
This form is common and there are ways to solve. All I need is to build the matrix A and the vector b. math.js has an excellent function called lusolve which will solve the linear system.
For one dimensional case, I used the lusolve because it is easily implemented. Alternatively numerical methods such as Gauss-Seidel method may be used to solve the linear system, which I used for two dimensional case.
Before Begin, Getting Started with JavaScript
The html file should like like this below
<!DOCTYPE HTML>
<html>
<head>
<script src="https://cdnjs.cloudflare.com/ajax/libs/mathjs/5.0.0/math.min.js"></script>
<script src="https://cdn.plot.ly/plotly-2.0.0.min.js"></script>
</head>
<body>
<div id="myDiv" style="width:603px; height:500px"></div>
<script src="scratch.js"></script>
</body>
</html>I’m using two libraries, Math.js for matrix solver and Plotly.js for plotting. Inside the body tag, div tag with id named “myDiv” was prepared so Plotly can plot. We will write main code in scratch.js file.
Coding
Let’s consider the conditions such that:
- n = 10 (so that grid point from 0 to 10)
- u(x0 = 0) = b0 = 0
- u(xn = 10) = bn = 10
Hence,
let n= 10;
let x0 = 0;
let xn = 10;
let b0 = 0;
let bn = 10;x values of each node can be found from node number and Δx, or dx.
let dx = (xn - x0) / n;
let x = new Array(n + 1);
for (i = 0; i < x.length; i++) { x[i] = x0 + i * dx };The function lusolve requires parameters A and b where A is an [n x n] matrix and b is a [n] column vector. Then there are 11 equations to solve (n = 0 to n = 10 as explained above). The matrix A and the vector b can be created as such
let A = [];
A.push(new Array(n + 1).fill(0));
A[0][0] = 1;
for (i = 1; i < n; i++) {
A.push(new Array(n + 1).fill(0));
A[i][i - 1] = 1;
A[i][i] = -2;
A[i][i + 1] = 1;
}
A.push(new Array(n + 1).fill(0));
A[n][n] = 1; let b = new Array(n + 1).fill(0);
b[0] = b0;
b[b.length - 1] = bn;
At this moment, if you open the index.html from file explorer (in Windows), browser will run the script and should output white screen. Ctrl+Shift+I will open the developer tool for chrome browser, under the console tab, type console.log(A), and console.log(b). It will show the content of the Arrays.
console.log(A)
(11) [Array(11), Array(11), Array(11), Array(11), Array(11), Array(11), Array(11), Array(11), Array(11), Array(11), Array(11)]
0: (11) [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
1: (11) [1, -2, 1, 0, 0, 0, 0, 0, 0, 0, 0]
2: (11) [0, 1, -2, 1, 0, 0, 0, 0, 0, 0, 0]
3: (11) [0, 0, 1, -2, 1, 0, 0, 0, 0, 0, 0]
4: (11) [0, 0, 0, 1, -2, 1, 0, 0, 0, 0, 0]
5: (11) [0, 0, 0, 0, 1, -2, 1, 0, 0, 0, 0]
6: (11) [0, 0, 0, 0, 0, 1, -2, 1, 0, 0, 0]
7: (11) [0, 0, 0, 0, 0, 0, 1, -2, 1, 0, 0]
8: (11) [0, 0, 0, 0, 0, 0, 0, 1, -2, 1, 0]
9: (11) [0, 0, 0, 0, 0, 0, 0, 0, 1, -2, 1]
10: (11) [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1]
length: 11
__proto__: Array(0)console.log(b)
(11) [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 10]
It looks good. Matrix A and vector b are created as intended. Use lusolve to get the solution u
let u = math.lusolve(A, b);Let’s visualize the solution. I found the plotly.js is a good library for plotting. Now in order to plot, we need to create an object
let trace_u = { x: [], y: [] }then add the x values to x, the solution u to the y array. Only thing is the solution u is given in the form of size 1 array in size n array. So,
for (let i = 0; i < n+1; i++) {
trace_u.x.push(x[i]);
trace_u.y.push(u[i][0]);
}to plot, pointing to a div object with id=’myDiv’ and the array of trace_u just created,
Plotly.newPlot('myDiv', [trace_u])
This will give the straight line passing through (0, 0) and (10, 10) coordinates. The Laplace equation in 1D is a straight line as expected.
Poisson equation in 1D
It turned out coding portion of the Poisson equation is quite simple. But before begin, let’s see the difference.
Set Boundary condition and Build Matrix
Let’s set the boundary condition as
Considering the boundary condition, and applying the equation, there are n+1 equations to solve
we can write the system of equations in the matrix form,
Let’s implement this in code.
Coding
Let’s consider the conditions such that:
- n = 18
- r = 2
- u(x = 0) = b0 =4
- u(x = 4) = bn = 4
Modify the code from Laplace equation slightly
let n= 18;
let r = 2;
let x0 = 0;
let xn = 4;
let b0 = 4;
let bn = 4;Now remember for Laplace equation, b array was filled with 0. Instead, fill with rΔx²
let b = new Array(n + 1).fill(r * dx ** 2);This is it. Plot it and you will see
The analytical solution to the Poisson equation is
where C1 and C2 are constant. Since r = 2, and applying the boundary condition, you will have
The parabolic curve, at minimum of 0 at x=2, and it satisfies the boundary condition.