SimClass/BlogRungeKuttaMethods.html at master · blackencino/SimClass · GitHub

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
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
<!DOCTYPE html>
<html>
<head>
  <title>Encino Sim Class : Runge-Kutta Integration Methods</title>

  <style type="text/css" title="currentStyle">
    @import "css/SimClass.css";
  </style>

<script src='http://processingjs.org/js/processing.min.js' type='text/javascript'/></script>
<!Install the MathJax stuff so we can show LaTeX>
<script type="text/javascript"
  src="https://c328740.ssl.cf1.rackcdn.com/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>
</head>
<body>
<h2>Higher-Order Time Steps</h2>
In our
<a href="http://encinographic.blogspot.com/2013/05/simulation-class-midpoint-method-and.html">previous class</a>, we improved the accuracy and stability of our
simulation's time step by exploring a number of different integration techniques,
of which the best-performing were the Euler-Cromer and Velocity-Verlet
techniques. However, while stable, their solutions drifted from the
theoretically correct
answer over time.
<p>
In each of our previous examples, the system's equations were evaluated just
once per time-step, or in some cases twice, which means the solutions were
first- or second- order. In order to improve our estimation of the function
we're trying to solve for, we can make more than one or two evaluations per
time-step, resulting in a higher-order approximation.
<p>
<a href="http://en.wikipedia.org/wiki/Runge–Kutta_methods">
Runge-Kutta</a> methods are a general class of explicit and implicit
iterative methods for approximating functions from their derivatives. In order
to explore this class of techniques, let's quickly revisit the equations
which define our simple harmonic oscillator system.  We have a mass attached
to an ideal spring with a certain stiffness, without friction, resulting in
a force equal to a constant times its displacement from the spring's origin,
in the direction opposite the displacement. That system is defined like this:
<div class="LaTexEquation">
    \(
    x'_t = v_t \\
    v'_t = a_t \\
    a_t = -\frac{k}{m} x_t
    \)
</div>
The idea behind all of the explicit Runge-Kutta methods is that we will
estimate a velocity, a position, and then an acceleration, and then we will
improve our estimation by using the previous estimates to make better estimates.
We can do this any number of times, and the number of these estimation stages
is what we mean by the "order" of the approximation. After all the estimation
steps have been taken, the state of the system is finally updated using a
weighted average of all the estimates.  The only difference between different
variations of the explicit Runge-Kutta methods is how the different estimates
are weighted, and how each estimate is used to produce the next estimate.
<p>
<h2>Second-Order Runge-Kutta</h2>
To illustrate the concept simply, here is the traditional
<a href="http://en.wikipedia.org/wiki/Runge–Kutta_methods#Second-order_methods_with_two_stages">
second-order
Runge-Kutta method</a>, applied to our simple system. In this notation, we
use the \(\star\) superscript to indicate an estimate, and the numbered
subscripts to indicate which iteration of estimation we're in. Note that each
estimation iteration only references values with subscripts less than its own:
this is what makes the method explicit, it computes new values only from
previous values.
<p>
<div class="LaTexEquation">
    \(
    \begin{array}{l}

    v^\star_1 = v_t, &
    a^\star_1 = -\frac{k}{m} x_t,  \\

    v^\star_2 = v_t + \Delta t\ a^\star_1, &
    a^\star_2 = -\frac{k}{m}( x_t + \Delta t\ v^\star_1 ), \\

    x_{t+\Delta t} = x_t + \frac{\Delta t}{2}(v^\star_1 + v^\star_2), &
    v_{t+\Delta t} = v_t + \frac{\Delta t}{2}(a^\star_1 + a^\star_2), \\

    \end{array}
    \)
</div>
<p>
In processing code, we can implement the RK2 timestep as follows:
<pre><code>
// Acceleration from Position.
float A_from_X( float i_x )
{
    return -( Stiffness / BobMass ) * i_x;
}

// Time Step function.
void TimeStep( float i_dt )
{
    float vStar1 = State[StateVelocityX];
    float aStar1 = A_from_X( State[StatePositionX] );

    float vStar2 = State[StateVelocityX] + ( i_dt * aStar1 );
    float xTmp = State[StatePositionX] + ( i_dt * vStar1 );
    float aStar2 = A_from_X( xTmp );

    State[StatePositionX] += ( i_dt / 2.0 ) * ( vStar1 + vStar2 );
    State[StateVelocityX] += ( i_dt / 2.0 ) * ( aStar1 + aStar2 );

    // Update current time.
    State[StateCurrentTime] += i_dt;
}
</code></pre>
<p>
Note that we've extracted the calculation of acceleration from position into
a separate function. This is to emphasize the structure of the system and not
bog the implementation down by repeated instances of the spring and mass
constants. In larger systems, this is usually the most expensive part of
the evaluation, so this implementation helps us understand how much actual
work we're doing. We can see above that we're evaluating the acceleration
twice in this time step.  Here's what that looks like running:
<p>
<script type="application/processing" data-processing-target="pjsSpringRK2">
float Stiffness = 5.0;
float BobMass = 0.5;
int StateSize = 3;
float[] InitState = new float[StateSize];
float[] State = new float[StateSize];
int StateCurrentTime = 0;
int StatePositionX = 1;
int StateVelocityX = 2;

int WindowWidthHeight = 300;
float WorldSize = 2.0;
float PixelsPerMeter;
float OriginPixelsX;
float OriginPixelsY;

void setup()
{
    // Create initial state.
    InitState[StateCurrentTime] = 0.0;
    InitState[StatePositionX] = 0.65;
    InitState[StateVelocityX] = 0.0;

    // Copy initial state to current state.
    // notice that this does not need to know what the meaning of the
    // state elements is, and would work regardless of the state's size.
    for ( int i = 0; i < StateSize; ++i )
    {
        State[i] = InitState[i];
    }

    // Set up normalized colors.
    colorMode( RGB, 1.0 );

    // Set up the stroke color and width.
    stroke( 0.0 );
    //strokeWeight( 0.01 );

    // Create the window size, set up the transformation variables.
    size( WindowWidthHeight, WindowWidthHeight );
    PixelsPerMeter = (( float )WindowWidthHeight ) / WorldSize;
    OriginPixelsX = 0.5 * ( float )WindowWidthHeight;
    OriginPixelsY = 0.5 * ( float )WindowWidthHeight;

    textSize( 24 );
}

// Draw our State, with the unfortunate units conversion.
void DrawState()
{
    // Compute end of arm.
    float SpringEndX = PixelsPerMeter * State[StatePositionX];

    // Compute the CORRECT position.
    float sqrtKoverM = sqrt( Stiffness / BobMass );
    float x0 = InitState[StatePositionX];
    float v0 = InitState[StateVelocityX];
    float t = State[StateCurrentTime];
    float CorrectPositionX = ( x0 * cos( sqrtKoverM * t ) ) +
        ( ( v0 / sqrtKoverM ) * sin( sqrtKoverM + t ) );

    // Compute draw pos for "correct"
    float CorrectEndX = PixelsPerMeter * CorrectPositionX;

    // Draw the spring.
    strokeWeight( 1.0 );
    line( 0.0, 0.0, SpringEndX, 0.0 );

    // Draw the spring pivot
    fill( 0.0 );
    ellipse( 0.0, 0.0,
             PixelsPerMeter * 0.03,
             PixelsPerMeter * 0.03 );

    // Draw the spring bob
    fill( 1.0, 0.0, 0.0 );
    ellipse( SpringEndX, 0.0,
             PixelsPerMeter * 0.1,
             PixelsPerMeter * 0.1 );

    // Draw the correct bob in blue
    fill( 0.0, 0.0, 1.0 );
    ellipse( CorrectEndX, -PixelsPerMeter * 0.25,
             PixelsPerMeter * 0.1,
             PixelsPerMeter * 0.1 );
}

// Acceleration from Position.
float A_from_X( float i_x )
{
    return -( Stiffness / BobMass ) * i_x;
}

// Time Step function.
void TimeStep( float i_dt )
{
    float vStar1 = State[StateVelocityX];
    float aStar1 = A_from_X( State[StatePositionX] );

    float vStar2 = State[StateVelocityX] + ( i_dt * aStar1 );
    float xTmp = State[StatePositionX] + ( i_dt * vStar1 );
    float aStar2 = A_from_X( xTmp );

    State[StatePositionX] += ( i_dt / 2.0 ) * ( vStar1 + vStar2 );
    State[StateVelocityX] += ( i_dt / 2.0 ) * ( aStar1 + aStar2 );

    // Update current time.
    State[StateCurrentTime] += i_dt;
}

// Processing Draw function, called every time the screen refreshes.
void draw()
{
    // Time Step.
    TimeStep( 1.0/24.0 );

    // Clear the display to a constant color.
    background( 0.75 );

    // Label.
    fill( 1.0 );
    text( "RK2", 10, 30 );

    pushMatrix();

    // Translate to the origin.
    translate( OriginPixelsX, OriginPixelsY );

    // Draw the simulation
    DrawState();
}

// Reset function. If the key 'r' is released in the display,
// copy the initial state to the state.
void keyReleased()
{
    if ( key == 114 )
    {
        for ( int i = 0; i < StateSize; ++i )
        {
            State[i] = InitState[i];
        }
    }
}

</script>
<canvas id="pjsSpringRK2"> </canvas>
<p>
RK2 is basically the same as the midpoint method, applied to both the
position and velocity estimates.
<p>
It's not as bad as forward Euler, but it diverges from the solution quickly,
and if you watch it
long enough, you will see that it is not entirely stable, and will eventually
produce large amplitudes.  We could apply the same type of forward-backwards
approach that we applied to Euler to get Euler-Cromer to RK2, but the results
are still not a good fit, though then they would be at least stable again. We
are simply using the second-order formulation to illustrate the technique
behind Runge-Kutta methods, so let's explore what happens at higher orders.
<p>
<h2>Fourth-Order Runge-Kutta</h2>
As we increase the number of estimations, the solution improves
dramatically. The
<a href="http://en.wikipedia.org/wiki/Runge–Kutta_methods#Common_fourth-order_Runge.E2.80.93Kutta_method">
fourth-order Runge-Kutta</a> approximation is very popular, and with good
reason. It makes four sets of estimates, as you'd expect from a fourth-order
method, and then uses an unevenly weighted average of those estimates to
get the final result. The
derivation of exactly how RK4 chooses its estimation timesteps, and how it
weights the average of its estimates, is complex. If you're interested in
working through the
<a href="https://en.wikipedia.org/wiki/Taylor_series">Taylor-series</a>
expansion, it is detailed
<a href="http://en.wikipedia.org/wiki/Runge–Kutta_methods#Derivation_of_the_Runge.E2.80.93Kutta_fourth-order_method">here</a>. For now, we'll trust the coefficients
that are given and look at how it works.  Here's the formulation, as applied
to our system:
<p>
<div class="LaTexEquation">
    \(
    \begin{array}{l}

    v^\star_1 = v_t, &
    a^\star_1 = -\frac{k}{m} x_t, \\

    v^\star_2 = v_t + \frac{\Delta t}{2}\ a^\star_1, &
    a^\star_2 = -\frac{k}{m}( x_t + \frac{\Delta t}{2}\ v^\star_1 ), \\

    v^\star_3 = v_t + \frac{\Delta t}{2}\ a^\star_2, &
    a^\star_3 = -\frac{k}{m}( x_t + \frac{\Delta t}{2}\ v^\star_2 ), \\

    v^\star_4 = v_t + \Delta t\ a^\star_3, &
    a^\star_4 = -\frac{k}{m}( x_t + \Delta t\ v^\star_3 ), \\

    x_{t+\Delta t} = x_t + \frac{\Delta t}{6}
        (v^\star_1 + 2 v^\star_2 + 2 v^\star_3 + v^\star_4), &
    v_{t+\Delta t} = v_t + \frac{\Delta t}{6}
        (a^\star_1 + 2 a^\star_2 + 2 a^\star_3 + a^\star_4), \\

    \end{array}
    \)
</div>
Here's what that looks like in Processing code, computing temporary positions
at each sub-estimation, and again using the separate function for acceleration
from position:
<p>
<pre><code>
// Time Step function.
void TimeStep( float i_dt )
{
    float vStar1 = State[StateVelocityX];
    float aStar1 = A_from_X( State[StatePositionX] );

    float vStar2 = State[StateVelocityX] + ( ( i_dt / 2.0 ) * aStar1 );
    float xTmp2 = State[StatePositionX] + ( ( i_dt / 2.0 ) * vStar1 );
    float aStar2 = A_from_X( xTmp2 );

    float vStar3 = State[StateVelocityX] + ( ( i_dt / 2.0 ) * aStar2 );
    float xTmp3 = State[StatePositionX] + ( ( i_dt / 2.0 ) * vStar2 );
    float aStar3 = A_from_X( xTmp3 );

    float vStar4 = State[StateVelocityX] + ( i_dt * aStar3 );
    float xTmp4 = State[StatePositionX] + ( i_dt * vStar3 );
    float aStar4 = A_from_X( xTmp4 );

    State[StatePositionX] += ( i_dt / 6.0 ) *
        ( vStar1 + (2.0*vStar2) + (2.0*vStar3) + vStar4 );
    State[StateVelocityX] += ( i_dt / 6.0 ) *
        ( aStar1 + (2.0*aStar2) + (2.0*aStar3) + aStar4 );

    // Update current time.
    State[StateCurrentTime] += i_dt;
}
</code></pre>
Here's what RK4 looks like running in Processing:
<p>
<script type="application/processing" data-processing-target="pjsSpringRK4">
float Stiffness = 5.0;
float BobMass = 0.5;
int StateSize = 3;
float[] InitState = new float[StateSize];
float[] State = new float[StateSize];
int StateCurrentTime = 0;
int StatePositionX = 1;
int StateVelocityX = 2;

int WindowWidthHeight = 300;
float WorldSize = 2.0;
float PixelsPerMeter;
float OriginPixelsX;
float OriginPixelsY;

void setup()
{
    // Create initial state.
    InitState[StateCurrentTime] = 0.0;
    InitState[StatePositionX] = 0.65;
    InitState[StateVelocityX] = 0.0;

    // Copy initial state to current state.
    // notice that this does not need to know what the meaning of the
    // state elements is, and would work regardless of the state's size.
    for ( int i = 0; i < StateSize; ++i )
    {
        State[i] = InitState[i];
    }

    // Set up normalized colors.
    colorMode( RGB, 1.0 );

    // Set up the stroke color and width.
    stroke( 0.0 );
    //strokeWeight( 0.01 );

    // Create the window size, set up the transformation variables.
    size( WindowWidthHeight, WindowWidthHeight );
    PixelsPerMeter = (( float )WindowWidthHeight ) / WorldSize;
    OriginPixelsX = 0.5 * ( float )WindowWidthHeight;
    OriginPixelsY = 0.5 * ( float )WindowWidthHeight;

    textSize( 24 );
}

// Draw our State, with the unfortunate units conversion.
void DrawState()
{
    // Compute end of arm.
    float SpringEndX = PixelsPerMeter * State[StatePositionX];

    // Compute the CORRECT position.
    float sqrtKoverM = sqrt( Stiffness / BobMass );
    float x0 = InitState[StatePositionX];
    float v0 = InitState[StateVelocityX];
    float t = State[StateCurrentTime];
    float CorrectPositionX = ( x0 * cos( sqrtKoverM * t ) ) +
        ( ( v0 / sqrtKoverM ) * sin( sqrtKoverM + t ) );

    // Compute draw pos for "correct"
    float CorrectEndX = PixelsPerMeter * CorrectPositionX;

    // Draw the spring.
    strokeWeight( 1.0 );
    line( 0.0, 0.0, SpringEndX, 0.0 );

    // Draw the spring pivot
    fill( 0.0 );
    ellipse( 0.0, 0.0,
             PixelsPerMeter * 0.03,
             PixelsPerMeter * 0.03 );

    // Draw the spring bob
    fill( 1.0, 0.0, 0.0 );
    ellipse( SpringEndX, 0.0,
             PixelsPerMeter * 0.1,
             PixelsPerMeter * 0.1 );

    // Draw the correct bob in blue
    fill( 0.0, 0.0, 1.0 );
    ellipse( CorrectEndX, -PixelsPerMeter * 0.25,
             PixelsPerMeter * 0.1,
             PixelsPerMeter * 0.1 );
}

// Acceleration from Position.
float A_from_X( float i_x )
{
    return -( Stiffness / BobMass ) * i_x;
}

// Time Step function.
void TimeStep( float i_dt )
{
    float vStar1 = State[StateVelocityX];
    float aStar1 = A_from_X( State[StatePositionX] );

    float vStar2 = State[StateVelocityX] + ( ( i_dt / 2.0 ) * aStar1 );
    float xTmp2 = State[StatePositionX] + ( ( i_dt / 2.0 ) * vStar1 );
    float aStar2 = A_from_X( xTmp2 );

    float vStar3 = State[StateVelocityX] + ( ( i_dt / 2.0 ) * aStar2 );
    float xTmp3 = State[StatePositionX] + ( ( i_dt / 2.0 ) * vStar2 );
    float aStar3 = A_from_X( xTmp3 );

    float vStar4 = State[StateVelocityX] + ( i_dt * aStar3 );
    float xTmp4 = State[StatePositionX] + ( i_dt * vStar3 );
    float aStar4 = A_from_X( xTmp4 );

    State[StatePositionX] += ( i_dt / 6.0 ) *
        ( vStar1 + (2.0*vStar2) + (2.0*vStar3) + vStar4 );
    State[StateVelocityX] += ( i_dt / 6.0 ) *
        ( aStar1 + (2.0*aStar2) + (2.0*aStar3) + aStar4 );

    // Update current time.
    State[StateCurrentTime] += i_dt;
}

// Processing Draw function, called every time the screen refreshes.
void draw()
{
    // Time Step.
    TimeStep( 1.0/24.0 );

    // Clear the display to a constant color.
    background( 0.75 );

    // Label.
    fill( 1.0 );
    text( "RK4", 10, 30 );

    pushMatrix();

    // Translate to the origin.
    translate( OriginPixelsX, OriginPixelsY );

    // Draw the simulation
    DrawState();
}

// Reset function. If the key 'r' is released in the display,
// copy the initial state to the state.
void keyReleased()
{
    if ( key == 114 )
    {
        for ( int i = 0; i < StateSize; ++i )
        {
            State[i] = InitState[i];
        }
    }
}

</script>
<canvas id="pjsSpringRK4"> </canvas>
<p>
Well hot-damn. The popularity of RK4 is easy to understand! However, it comes
at a cost - four evaluations of acceleration in the timestep. Futhermore, there
is a temporary storage requirement which begins to be an issue in larger
simulations where a temporary position or velocity may amount to millions of
points. Also, RK4 is not unconditionally stable.
In many production cases, the stability of Velocity Verlet or Euler-Cromer may
be more attractive (especially with the fewer calculations), and the additional
accuracy provided by RK4 may not be visible.  These are the trade-offs
and design decisions made by simulation designers when implementing simulation
engines.
<p>
Interesting point - the "correct" solution to our simple system here is just a
cosine function of time. Our numerical integrator is achieving almost identical
results without ever evaluating anything other than a multiplication or an
addition. This illustrates how
<a href="https://en.wikipedia.org/wiki/Taylor_series">Taylor-series</a>
 expansions are used to provide
numerical approximations to mathematical functions.
<p>
So far, we've only been considering differential equations of a single variable,
time. Now that we've achieved a pretty good time step, the next class will
explore a system in which the equations feature derivatives in both time and
space.
</body>