-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathindex.html
242 lines (229 loc) · 12.3 KB
/
index.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="UTF-8">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<link rel="preconnect" href="https://api.fonts.coollabs.io">
<title>Numerical Integration - Solving ODEs in JS</title>
<meta name="description" content="Solve ODEs using JavaScript and the Runge Kutta Method">
<meta property="og:image" content="https://jurasic.dev/images/ode.png">
<meta property="og:type" content="blog" />
<link rel="stylesheet" href="css/ode.css">
<link rel="stylesheet" href="https://api.fonts.coollabs.io/css2?family=Raleway:wght@300;400;600&display=swap">
<link rel="preload" href="css/uplot.min.css" as="style">
<link rel="preload" href="scripts/ode.js" as="script">
<link rel="preload" href="scripts/uplot.min.js" as="script">
<link rel="preload" href="css/prism.css" as="style">
<link rel="preload" href="scripts/prism.js" as="script">
</head>
<body>
<a class="homeLink" href="/"><svg width="40" height="40" viewBox="0 0 24 24" fill="none"
xmlns="http://www.w3.org/2000/svg">
<path
d="M11.9481 14.8285L10.5339 16.2427L6.29126 12L10.5339 7.7574L11.9481 9.17161L10.1197 11H17.6568V13H10.1197L11.9481 14.8285Z"
fill="currentColor" />
<path fill-rule="evenodd" clip-rule="evenodd"
d="M23 19C23 21.2091 21.2091 23 19 23H5C2.79086 23 1 21.2091 1 19V5C1 2.79086 2.79086 1 5 1H19C21.2091 1 23 2.79086 23 5V19ZM19 21H5C3.89543 21 3 20.1046 3 19V5C3 3.89543 3.89543 3 5 3H19C20.1046 3 21 3.89543 21 5V19C21 20.1046 20.1046 21 19 21Z"
fill="currentColor" />
</svg></a>
<h1>Numerical solutions to Ordinary differential equations</h1>
<h2>Basics</h2>
<p>
Ordinary differential equations are equations that involve ordinary derivatives of a function and the function
itself. The order of the differential equation is characterized by the highest order derivative it contains. ODEs
are used to model a variety of systems like the stock market, orbiting planets or predator prey systems. In this
article, we will explore a few numerical methods for solving the following first order equation in JavaScript.
</p>
<object style="width: 100%; max-height: 3.2em;" data="images/ode.svg" type="image/svg+xml">ordinary differential
equation</object>
<h2>Numerical methods</h2>
<p>
Since many differential equations are hard to solve mathematically, numerical methods have been developed to compute
aproximate solutions.
</p>
<h3 id="euler">Euler integration</h3>
<p>
Euler integration is the simplest method for numerical integration. It is the equivalent of following the slope of a
function in tiny steps. At each timestep, we calculate the current derivative, multiply it by the step size <b>h</b>
to get the change in y-direction, and add it to the previous y value to arrive at the next state.
<!-- y_{n+1} = y_n + h \cdot f(t,y_n) -->
</p>
<object style="width: 100%; max-height: 1.6em;" data="images/eulerStep.svg" type="image/svg+xml">Euler iteration
step</object>
<p>
In the following JS implementation of the Euler method <b>ts</b> is an array of
linearly spaced values ranging from the initial time <b>t_0</b> to the end time <b>t_1</b>. To store the resulting
values, we create an empty array <b>ys</b> and initialize the first entry to the initial ODE condition <b>y_0</b>.
</p>
<pre><code class="language-javascript" contenteditable spellcheck="false">const N = 100, t_0 = 0, t_1 = 1, y_0 = 2
const h = (t_1 - t_0) / N //time step size
var ts = Array.from(Array(N+1), (_, k) => k * h + t_0)
var ys = Array(N+1).fill(0) //empty array for the results
ys[0] = y_0 //initial conditions
for (let i = 0; i < N; i++) {
ys[i + 1] = ys[i] + f(ts[i], ys[i]) * h
}</code></pre>
<p>
The biggest drawback of the Euler method is, that the error scales with O(h²) and therefore oftentimes isn't
stable. An example of the emerging instability can be seen <a href="oscillator">in this demo</a> when the resolution
gets decreased.
</p>
<h3 id="midpoint">Midpoint method</h3>
<p>
The <a href="https://en.wikipedia.org/wiki/Midpoint_method">midpoint method</a> improves upon the Euler method, by
taking the "average slope" between the current location and the next one, instead of the slope at the current
location. As the name suggests, it determines the derivative at the midpoint between two steps. To do this, a
regular Euler iteration with half step size is performed and used to evaluate the slope. The next y value is
calculated by moving along the slope by a full step from the previous position, but this time using the "midpoint
slope".
</p>
<object style="width: 100%; max-height: 6em;" data="images/midpointStep.svg" type="image/svg+xml">Midpoint iteration
step</object>
<p>Using this method, the Order of the error decreases to O(h³) per step and O(h²)
for the total solution. The midpoint method increases stability significantly, while only requiring two evaluations
of the ODE per step. In realtime applications, this is often a good balance between accuracy and performance. Here's
the sample implementation in JavaScript:
</p>
<pre><code class="language-javascript" contenteditable spellcheck="false">for (let i = 0; i < resolution; i++) {
const k1 = f(ts[i], ys[i])
const k2 = f(ts[i] + h/2, ys[i] + k1 * h/2)
ys[i + 1] = ys[i] + k2 * h
}</code></pre>
<h3 id="rk4">Runge Kutta 4th Order </h3>
<p>
Now we're getting into the fun stuff.
The <a href="https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods">Runge Kutta family</a> of methods takes the
concept further, by introducing n evaluations, each more accurate than the previous one, to compute a single step.
For the vast majority of applications, the fourth order method also called the standard Runge Kutta Method (RK4),
is accurate enough. Even at large step sizes it retains a high stability, since the error per step only scales with
O(h^5)
</p>
<object style="width: 100%; max-height: 9em;" data="images/rk4.svg" type="image/svg+xml">runge kutta
equations</object>
<p>
The RK4 method really starts to shine when working with periodic solutions and nonlinear ODEs, since the lower order
methods oftentimes diverge in those situation, even for low step sizes.
<br>
A simple JS implementation for the RK4 method could look like this:
</p>
<pre><code class="language-javascript" contenteditable spellcheck="false">for (let i = 0; i < resolution; i++) {
const k1 = f(ts[i], ys[i])
const s1 = ys[i] + k1 * h/2
const k2 = f(ts[i] + h/2, s1)
const s2 = ys[i] + k2 * h/2
const k3 = f(ts[i] + h/2, s2)
const s3 = ys[i] + k3 * h
const k4 = f(ts[i] + h, s3) // f(t + h, y_n + k3*h)
ys[i + 1] = ys[i] + (k1/6 + k2/3 + k3/3 + k4/6) * h
}</code></pre>
<hr>
<p>The demo below demonstrates how the accuracy of all three methods scales with the step size. Use the slider below
to change their resolution!
<noscript>
<br>
Please enable JavaScript to see the demo. I promise it's not for tracking.
</noscript>
</p>
<label for="resolution_rng">resolution</label>
<input type="range" oninput="u.setData(calcOde(Math.pow(2,this.value)))" min="2" max="16" name="resolution slider"
id="resolution_rng">
<div id="plot"></div>
<h2>The Limits of Resolution</h2>
<p>As you can see, a smaller step size leads to more accurate results. This way of increasing accuracy unfortunately
has it's limits. At some point, floating point <a href="https://en.wikipedia.org/wiki/Round-off_error">roundoff
errors</a> from computers limited numeric resolution start outweighing the improvements a smaller step size gives
us. To get even higher accuracies, higher order methods have to be used.</p>
<a title="Berland, Public domain, via Wikimedia Commons"
href="https://commons.wikimedia.org/wiki/File:AbsoluteErrorNumericalDifferentiationExample.png">
<img width="1024" alt="AbsoluteErrorNumericalDifferentiationExample"
src="https://upload.wikimedia.org/wikipedia/commons/thumb/4/41/AbsoluteErrorNumericalDifferentiationExample.png/1024px-AbsoluteErrorNumericalDifferentiationExample.png"></a>
<h2>Moving into higher orders</h2>
<p>Until now, we have only looked at solving first order differential equations. To solve higher order ODEs like
the <a href="https://en.wikipedia.org/wiki/Electromagnetic_wave_equation">wave equation</a> we first have to turn
it into a system of first order ODEs. As an example, let's look at harmonic oscillators.
</p>
<object style="width: 100%; max-height: 1.8em;" data="images/oszi.svg" type="image/svg+xml">harmonic oscillator
equation</object>
<p>
To turn this into a system of first order ODEs, let's write the second order derivative <b>y''</b> as a derivative
of <b>y'</b>. Now we can represent <b>y</b> as a vector consisting of the value and it's first derivative like so
</p>
<object style="width: 100%; max-height: 3.2em;" data="images/firstorder.svg" type="image/svg+xml">system of first
order euqations (harmonic oscillator)</object>
<p>
This reqires some changes in the implementation, since <b>y</b> isn't a number anymore, but rather a vector.
To see one possible implementation, either inspect the site or take a look at the
<a href="https://github.com/missing-user/missing-user.github.io">GitHub repository</a> for all the demos. There's a
bunch of other improvements that could be done, like moving from a fixed resolution to an <a
href="https://en.wikipedia.org/wiki/Adaptive_step_size">adaptive step size</a> like I
did in this <a href="https://missing-user.github.io/solarSystem">solar system simulation</a> or implementing a
<a href="https://en.wikipedia.org/wiki/Linear_multistep_method#Adams%E2%80%93Bashforth_methods">multistep
method</a>, that takes more than one previous point into consideration for the calculation.
</p>
<h2>...and beyond</h2>
<p>Using these methods, you can now solve a variety of problems from physics, mathematics and even biology. Here's a
list of problems I find particularly interesting and an interactive implementation to go along with each one. If you
enjoy these kinds of projects, why don't you check out <a href="https://jurasic.dev/">my website</a>!</p>
<ul>
<li><a href="oscillator">Harmonic oscillators</a></li>
<li><a href="airy">Airy functions</a></li>
<li><a href="lotkavolterra">Lotka Volterra equations</a></li>
<li><a href="lorenz">Lorenz Attractors</a></li>
<li><a href="aizawa">Aizawa Attractors</a></li>
<li><a href="arnedo">Arnedo Attractors</a></li>
</ul>
</body>
<link href="css/prism.css" rel="stylesheet" type="text/css" />
<script src="scripts/prism.js" type="text/javascript"></script>
<link rel="stylesheet" href="css/uplot.min.css">
<script src="scripts/uplot.min.js"></script>
<script src="scripts/ode.js"></script>
<script async>
var plotElem = document.getElementById("plot")
const DefaultODE = {
f: (t, y) => [2 * Math.sin(3 * t) + y[0] / 4],
parameters: [[0.2], 0, 8]
}
var solver = new ODEsolver(DefaultODE.f, ...DefaultODE.parameters)
const opts = {
width: 900,
height: 900,
title: "Numerical ODE Solution",
scales: {
x: { time: false },
y: { range: [0, 7] }
},
series: [
{ label: "t" },
{
label: "Euler",
stroke: "#D32F2F",
},
{
label: "Midpoint",
stroke: "#ff8000",
},
{
label: "RK4",
stroke: "#12B02F",
},
],
}
function calcOde(resolution) {
const eres = solver.euler(~~resolution)
data = [
eres.ts,
eres.ys.map((t) => t[0]),
solver.midpoint(~~resolution).ys.map((t) => t[0]),
solver.rk4(~~resolution).ys.map((t) => t[0]),
]
return data
}
var u = new uPlot(opts, calcOde(100), plotElem)
u.setSize({ width: plotElem.clientWidth, height: plotElem.clientWidth });
window.addEventListener("resize", e => {
u.setSize({ width: plotElem.clientWidth, height: plotElem.clientWidth });
})
</script>
<script>if (!!localStorage.getItem("_dev") && !sessionStorage.getItem("_swa") && document.referrer.indexOf(location.protocol + "//" + location.host) !== 0) { fetch("https://counter.dev/track?" + new URLSearchParams({ referrer: document.referrer, screen: screen.width + "x" + screen.height, user: "missing-user", utcoffset: "1" })) }; sessionStorage.setItem("_swa", "1");</script>
</html>