-
Notifications
You must be signed in to change notification settings - Fork 7
Expand file tree
/
Copy pathtime_integration.html
More file actions
491 lines (421 loc) · 18.9 KB
/
Copy pathtime_integration.html
File metadata and controls
491 lines (421 loc) · 18.9 KB
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
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
<!doctype html>
<html class="no-js" lang="en">
<head>
<meta charset="utf-8">
<meta name="viewport" content="width=device-width, initial-scale=1">
<link rel="stylesheet" href="https://interactivecomputergraphics.github.io/physics-simulation/examples/style.css">
<script type="text/x-mathjax-config">
MathJax.Hub.Config({
extensions: ["tex2jax.js"],
jax: ["input/TeX", "output/HTML-CSS"],
tex2jax: {
inlineMath: [ ['$','$'], ["\\(","\\)"] ],
displayMath: [ ['$$','$$'], ["\\[","\\]"] ],
processEscapes: true
},
"HTML-CSS": { fonts: ["TeX"] }
});
</script>
<script type="text/javascript" async src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.4/MathJax.js"></script>
<title>Time integration</title>
</head>
<body>
<header class="page-header">
<h1>Time integration</h1>
</header>
<main>
<!-- Simulation panel: canvas + controls -->
<table style="align_center;border-radius: 20px;padding: 20px;margin:auto">
<col width="70%">
<col width="30%">
<tr>
<td>
<div class="card sim-panel">
<div class="sim-canvas-wrap">
<canvas id="simCanvas" width="1024" height="960" style="border:2px solid #000000;border-radius: 20px;background-color:#EEEEEE">Your browser does not support the HTML5 canvas tag.</canvas>
</div>
</div>
</td>
<td>
<div class="controls-panel">
<h3>Controls</h3>
<div class="controls-grid">
<label>Current time</label>
<span class="stat-value"><span id="time">0.00</span> s</span>
<label for="timeStepSizeInput">Time step size</label>
<input onchange="gui.restart()" id="timeStepSizeInput" type="number" value="0.05" step="0.01">
<label for="timeIntegrationInput">Time integration method</label>
<select onchange="gui.restart()" id="timeIntegrationInput">
<option>Explicit Euler</option>
<option selected="selected">Symplectic Euler</option>
<option>Runge-Kutta 2</option>
</select>
<label for="gravityInput">Gravity</label>
<input onchange="gui.restart()" id="gravityInput" type="number" value="-9.81" step="0.01">
<label for="massInput">Mass</label>
<input onchange="gui.restart()" id="massInput" type="number" value="1.0" step="0.01">
<div class="full-width">
<button onclick="gui.restart()" id="restart">▶ Restart</button>
</div>
<div class="full-width">
<button onclick="gui.doPause()" id="Pause" class="btn-secondary">▮▮ Pause</button>
</div>
</div>
</div>
</td>
</tr>
</table>
<!-- Theory section -->
<div class="card theory">
<h2>Time integration methods:</h2>
Simulating the motion of a physical system means solving its equations of motion, a system of ordinary differential equations (ODEs) that relate positions, velocities and accelerations. Since these equations can generally not be solved analytically, we approximate the solution step by step: given the state (position and velocity) at time $t$, a <b>time integration method</b> computes an approximation of the state at time $t + \Delta t$ using the currently acting forces. Repeating this update for a sequence of small time steps $\Delta t$ advances the simulation over time.
<p>
This example shows the motion of a single particle that has an initial velocity and is subject to a constant gravitational acceleration. Because the acceleration is constant here, this special case is one of the rare ones that admits a closed-form solution:
$$\begin{align*}
\mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t) + \frac12 \Delta t^2 \mathbf a^\text{grav} \\
\mathbf v(t + \Delta t) &= \mathbf v(t) + \Delta t \mathbf a^\text{grav} .
\end{align*}$$
This exact trajectory is plotted in green and serves as ground truth against which the numerical solutions (blue), computed with the methods introduced below, can be compared. As soon as forces depend on the state itself in a more complex way — e.g. spring, contact or SPH pressure forces — no analytic solution is available anymore and numerical time integration becomes essential.
</p>
<p>
Three properties are commonly used to characterize and compare time integration methods:
<ul>
<li><b>Order of accuracy:</b> A method of order $p$ has a local error (introduced in a single step) of $O(\Delta t^{p+1})$ and, consequently, a global error (accumulated over a fixed time interval) of $O(\Delta t^{p})$. Halving $\Delta t$ for a method of order $p$ roughly reduces the global error by a factor of $2^p$.</li>
<li><b>Computational cost:</b> given by the number of evaluations of the acceleration/force function required per time step. This is important in practice since evaluating the forces (e.g. for cloth, deformable solids or SPH fluids) is usually the most expensive part of a simulation step.</li>
<li><b>Stability and energy behavior:</b> describes how the numerical solution behaves over long simulation times, in particular whether it stays close to the true trajectory or whether energy is spuriously gained, lost, or oscillates around the correct value.</li>
</ul>
Try the different methods in the simulation with a large time step (e.g. $\Delta t = 0.1$) and compare how far the blue (numerical) trajectory deviates from the green (exact) one.
</p>
<h3>Explicit Euler</h3>
The explicit Euler method (also called forward Euler) is the simplest time integration scheme for ODEs. It uses only quantities known at the current time $t$ to extrapolate the state to $t + \Delta t$:
$$\begin{align*}
\mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t) \\
\mathbf v(t + \Delta t) &= \mathbf v(t) + \Delta t \mathbf a(t).
\end{align*}$$
<p><b>Properties:</b></p>
<ul>
<li>First-order accurate ($O(\Delta t)$ global error); requires only a single force evaluation per step, so it is the cheapest of the three methods.</li>
<li><b>Not symplectic:</b> it systematically injects energy into the system because the position update uses the <em>old</em> velocity $\mathbf v(t)$ instead of the newly computed one. For an oscillating system (e.g. a mass on a spring) this causes the amplitude to grow over time, even though the exact solution conserves energy.</li>
<li>Only conditionally stable: for stiff forces (e.g. stiff springs) $\Delta t$ must be chosen small enough, otherwise the simulation diverges numerically. For this reason explicit Euler is rarely used in practice and mainly serves as a baseline for comparison.</li>
</ul>
<h3>Symplectic Euler</h3>
The symplectic Euler method is also known as semi-implicit Euler: it first updates the velocity using the current acceleration, and then uses this <em>already updated</em> velocity to advance the position:
$$\begin{align*}
\mathbf v(t + \Delta t) &= \mathbf v(t) + \Delta t \mathbf a(t) \\
\mathbf x(t + \Delta t) &= \mathbf x(t) + \Delta t \mathbf v(t + \Delta t).
\end{align*}$$
<p><b>Properties:</b></p>
<ul>
<li>Still only first-order accurate, just like explicit Euler, and has the same cost of a single force evaluation per step — switching from explicit to symplectic Euler is essentially "free" in terms of performance.</li>
<li><b>Symplectic:</b> the method (approximately) preserves phase-space volume, which implies that the energy error stays bounded and oscillates around the true value instead of drifting away over time. This yields qualitatively much better long-term behavior than explicit Euler, even though the two methods have the same order and the same error at any single instant.</li>
<li>Because it is cheap, simple to implement and behaves well over long simulations, symplectic Euler is the de-facto standard in real-time physics: it is used for particle systems, cloth, rigid bodies, SPH fluids and position-based dynamics alike.</li>
</ul>
<h3>Runge-Kutta 2</h3>
Higher-order integration methods are usually formulated for a general first-order ODE $\dot{\mathbf s}(t) = \mathbf f(t, \mathbf s(t))$, where $\mathbf s$ is the state of the system at time $t$. In our example the state consists of the particle's position and velocity, and the ODE is defined by:
$$\begin{equation*}
\mathbf f(t, \mathbf s(t)) = \begin{pmatrix} \dot{\mathbf x} \\ \dot{\mathbf v} \end{pmatrix} = \begin{pmatrix} \mathbf v(t) \\ \mathbf a(t) \end{pmatrix}, \quad\quad \mathbf s(t) = \begin{pmatrix} \mathbf x(t) \\ \mathbf v(t) \end{pmatrix}.
\end{equation*}$$
Using this notation, the second-order Runge-Kutta method (also known as the midpoint method) is defined as:
$$\begin{align*}
\mathbf k_1 &= \Delta t \, \mathbf f(t, \mathbf s(t)) \\
\mathbf k_2 &= \Delta t \, \mathbf f\!\left(t + \tfrac12 \Delta t,\; \mathbf s(t) + \tfrac12 \mathbf k_1\right) \\
\mathbf s(t + \Delta t) &= \mathbf s(t) + \mathbf k_2.
\end{align*}$$
Intuitively, $\mathbf k_1$ is a trial (explicit Euler) step used only to estimate the state at the midpoint $t + \frac12 \Delta t$; the actual step $\mathbf k_2$ is then taken using the derivative evaluated at that midpoint, which cancels out the leading-order error term of the Euler methods.
<p><b>Properties:</b></p>
<ul>
<li>Second-order accurate ($O(\Delta t^2)$ global error): halving $\Delta t$ roughly quarters the error, so for smooth, non-stiff problems RK2 reaches a given accuracy with a much larger time step than either Euler method.</li>
<li>Has two stages, i.e. the function $\mathbf f$ must be evaluated twice per time step — roughly twice the cost of an Euler step. In general, an order-$p$ Runge-Kutta method typically needs at least $p$ stages, so higher accuracy is traded against more force evaluations per step.</li>
<li><b>Not symplectic:</b> despite its higher accuracy, RK2 can still exhibit long-term energy drift on oscillatory problems; the higher order just makes the drift smaller and less visible over short simulation times.</li>
</ul>
<h3>Summary</h3>
<p>The following table summarizes the properties discussed above:</p>
<table>
<tr><th>Method</th><th>Order</th><th>Force evals / step</th><th>Symplectic</th><th>Typical use</th></tr>
<tr><td>Explicit Euler</td><td>1</td><td>1</td><td>No</td><td>Baseline / teaching example</td></tr>
<tr><td>Symplectic Euler</td><td>1</td><td>1</td><td>Yes</td><td>Default for real-time simulation (games, cloth, rigid bodies, SPH, PBD)</td></tr>
<tr><td>Runge-Kutta 2</td><td>2</td><td>2</td><td>No</td><td>Higher short-term accuracy for smooth, non-stiff problems</td></tr>
</table>
There is no single "best" method: the right choice depends on whether short-term accuracy, long-term (energy) stability, or per-step cost matters most for the application at hand. Higher-order methods such as Runge-Kutta 4, or implicit methods for stiff problems, extend these same trade-offs further and are covered in later examples.
<h3>References</h3>
<ul>
<li>[HLW06] Ernst Hairer, Christian Lubich, Gerhard Wanner. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations, 2nd edition. Springer, 2006.</li>
<li>[BW98] David Baraff, Andrew Witkin. Large Steps in Cloth Simulation. In Proceedings of SIGGRAPH, 1998.</li>
</ul>
</div>
</main>
<script id="simulation_code" type="text/javascript">
class Particle
{
constructor (x, y)
{
this.x = x;
this.y = y;
this.fx = 0.0;
this.fy = 0.0;
this.vx = 0.0;
this.vy = 0.0;
}
}
class Simulation
{
constructor()
{
this.particles = [];
this.mass = 1.0;
this.gravity = -9.81;
this.timeStepSize = 0.05;
this.time = 0;
this.endTime = 3.0;
this.solution = [];
this.reference_solution = [];
this.timeIntegrationMethod = "Symplectic Euler";
this.pause = true;
this.init();
}
init()
{
this.computeReferenceSolution();
// create particle
this.particles = [];
this.particles.push(new Particle(0.0, 0.0));
let p = this.particles[0];
p.vx = 1;
p.vy = 4.5;
this.time = 0.0;
this.solution = [];
}
computeReferenceSolution()
{
// create particle
this.particles = [];
this.particles.push(new Particle(0.0, 0.0));
let p = this.particles[0];
p.vx = 1;
p.vy = 4.5;
this.time = 0.0;
this.reference_solution = [];
while (this.time < this.endTime)
{
let p = this.particles[0];
this.reference_solution.push([p.x, p.y]);
this.exactIntegration();
this.time = this.time + this.timeStepSize;
}
}
exactIntegration()
{
let dt = this.timeStepSize;
// exact time integration for constant force
for (let i = 0; i < this.particles.length; i++)
{
let p = this.particles[i];
// integrate position
p.x = p.x + dt * p.vx + 0.5 * dt*dt* p.fx / this.mass;
p.y = p.y + dt * p.vy + 0.5 * dt*dt* (p.fy / this.mass + this.gravity);
// integrate velocity considering gravitational acceleration and forces
p.vx = p.vx + dt * p.fx / this.mass;
p.vy = p.vy + dt * (p.fy / this.mass + this.gravity);
}
}
explicitEuler()
{
let dt = this.timeStepSize;
// symplectic Euler step
for (let i = 0; i < this.particles.length; i++)
{
let p = this.particles[i];
// integrate position
p.x = p.x + dt * p.vx;
p.y = p.y + dt * p.vy;
// integrate velocity considering gravitational acceleration and forces
p.vx = p.vx + dt * p.fx / this.mass;
p.vy = p.vy + dt * (p.fy / this.mass + this.gravity);
}
}
symplecticEuler()
{
let dt = this.timeStepSize;
// symplectic Euler step
for (let i = 0; i < this.particles.length; i++)
{
let p = this.particles[i];
// integrate velocity considering gravitational acceleration and forces
p.vx = p.vx + dt * p.fx / this.mass;
p.vy = p.vy + dt * (p.fy / this.mass + this.gravity);
// integrate position
p.x = p.x + dt * p.vx;
p.y = p.y + dt * p.vy;
}
}
rungeKutta2()
{
let dt = this.timeStepSize;
// Runge-Kutta 2 step
for (let i = 0; i < this.particles.length; i++)
{
let p = this.particles[i];
let k1_x = dt * p.vx;
let k1_y = dt * p.vy;
let k1_vx = dt/this.mass * p.fx;
let k1_vy = dt*(p.fy / this.mass + this.gravity);
let k2_x = dt * (p.vx + 0.5*k1_vx);
let k2_y = dt * (p.vy + 0.5*k1_vy);
let k2_vx = dt/this.mass * p.fx;
let k2_vy = dt*(p.fy / this.mass + this.gravity);
p.x = p.x + k2_x;
p.y = p.y + k2_y;
p.vx = p.vx + k2_vx;
p.vy = p.vy + k2_vy;
}
}
// simulation step
simulationStep()
{
if (this.pause)
return;
let dt = this.timeStepSize;
let p = this.particles[0];
if (this.time < this.endTime)
this.solution.push([p.x, p.y]);
else
this.pause = true;
// time integration step
if (this.timeIntegrationMethod == "Explicit Euler")
this.explicitEuler();
else if (this.timeIntegrationMethod == "Symplectic Euler")
this.symplecticEuler();
else if (this.timeIntegrationMethod == "Runge-Kutta 2")
this.rungeKutta2();
// update simulation time
this.time = this.time + dt;
}
}
class GUI
{
constructor()
{
this.canvas = document.getElementById("simCanvas");
this.c = this.canvas.getContext("2d");
this.requestID = -1;
this.origin = { x : this.canvas.width / 2, y : this.canvas.height/2};
this.zoom = 150;
this.particleRadius = 0.025;
// register mouse event listeners (zoom/selection)
this.canvas.addEventListener("wheel", this.wheel.bind(this), false);
}
// set simulation parameters from GUI and start mainLoop
restart()
{
window.cancelAnimationFrame(this.requestID);
if (this.sim == undefined)
this.sim = new Simulation();
this.sim.gravity = parseFloat(document.getElementById('gravityInput').value);
this.sim.timeStepSize = parseFloat(document.getElementById('timeStepSizeInput').value);
this.sim.mass = parseFloat(document.getElementById('massInput').value);
this.sim.timeIntegrationMethod = document.getElementById('timeIntegrationInput').value;
this.sim.init();
this.mainLoop();
}
drawCoordinateSystem()
{
// draw x-axis
this.c.strokeStyle = "#FF0000";
this.c.beginPath();
this.c.moveTo(this.origin.x, this.origin.y);
this.c.lineTo(this.origin.x+1*this.zoom, this.origin.y);
this.c.stroke();
// draw y-axis
this.c.strokeStyle = "#00FF00";
this.c.beginPath();
this.c.moveTo(this.origin.x, this.origin.y);
this.c.lineTo(this.origin.x, this.origin.y-1*this.zoom);
this.c.stroke();
}
draw()
{
this.c.clearRect(0, 0, this.canvas.width, this.canvas.height);
this.drawCoordinateSystem();
// draw particles as circles
let r = this.particleRadius;
for (let i = 0; i < this.sim.particles.length; i++)
{
let p = this.sim.particles[i];
let r = this.particleRadius;
this.c.fillStyle = "#0000FF";
let px = this.origin.x + p.x * this.zoom;
let py = this.origin.y - p.y * this.zoom;
this.c.beginPath();
this.c.arc(px, py, r * this.zoom, 0, Math.PI*2, true);
this.c.closePath();
this.c.fill();
}
this.c.strokeStyle = "#00BB00";
this.c.fillStyle = "#00BB00";
for (let i = 0; i < this.sim.reference_solution.length; i++)
{
let x = this.sim.reference_solution[i];
let px = this.origin.x + x[0] * this.zoom;
let py = this.origin.y - x[1] * this.zoom;
this.c.beginPath();
this.c.arc(px, py, r * this.zoom, 0, Math.PI*2, true);
this.c.closePath();
this.c.fill();
if (i > 0)
{
let x2 = this.sim.reference_solution[i-1];
this.c.beginPath();
this.c.moveTo(this.origin.x + x[0]*this.zoom, this.origin.y - x[1]*this.zoom);
this.c.lineTo(this.origin.x + x2[0]*this.zoom, this.origin.y - x2[1]*this.zoom);
this.c.stroke();
}
}
this.c.strokeStyle = "#0011AA";
this.c.fillStyle = "#0011AA";
for (let i = 0; i < this.sim.solution.length; i++)
{
let x = this.sim.solution[i];
let px = this.origin.x + x[0] * this.zoom;
let py = this.origin.y - x[1] * this.zoom;
this.c.beginPath();
this.c.arc(px, py, r * this.zoom, 0, Math.PI*2, true);
this.c.closePath();
this.c.fill();
if (i > 0)
{
let x2 = this.sim.solution[i-1];
this.c.beginPath();
this.c.moveTo(this.origin.x + x[0]*this.zoom, this.origin.y - x[1]*this.zoom);
this.c.lineTo(this.origin.x + x2[0]*this.zoom, this.origin.y - x2[1]*this.zoom);
this.c.stroke();
}
}
}
mainLoop()
{
this.sim.simulationStep();
document.getElementById("time").innerHTML = this.sim.time.toFixed(2);
this.draw();
if (!this.sim.pause)
this.requestID = window.requestAnimationFrame(this.mainLoop.bind(this));
}
doPause()
{
this.sim.pause = !this.sim.pause;
if (!this.sim.pause)
this.mainLoop();
}
wheel(event)
{
event.preventDefault();
this.zoom += event.deltaY * -0.05;
if (this.zoom < 1)
this.zoom = 1;
this.draw();
}
}
gui = new GUI();
gui.restart();
</script>
</body>
</html>