Target audience: Advanced

Estimated reading time: 30'

This post describes the implementation of the different Runge-Kutta method to solve differential equations in Scala.

The objective is to leverage the functional programming components of the Scala programming language to create a generic solver of ordinary differential equations (ODE) using

The objective is to leverage the functional programming components of the Scala programming language to create a generic solver of ordinary differential equations (ODE) using

**Runge-Kutta**family of approximation algorithms.Overview

Most of ordinary differential equations cannot be solved analytically. In this case, a numeric approximation to the solution is
often good enough to solve an engineering problem. Oddly enough most of commonly used algorithm to compute such an approximation have been establish a century ago. Let's consider the differential equation \[\frac{\mathrm{d}y }{\mathrm{d} x} = f(x,y)\] The family of explicit Runge-Kutta numerical approximation methods is defined as \[y_{n+1} = y_{n} + \sum_{i=0}^{s<n}b_{i}k_{i}\\where\,\,k_{j}=h.f(x_{n} + c_{j}h, y_{n} + \sum_{s=1}^{j-1} a_{s,s-1}k_{s-1} )\\with\,\,h=x_{n+1}-x_{n}\,\,and\,\,\Delta = \frac{dy}{dx}+ \sum_{s=1}^{j-1} a_{s,s-1}k_{s-1}\] k(j) is the increment based on the slope at the midpoint of the interval [x(n),x(n+1)] using delta. The Euler method defined as \[y_{n+1} = y_{n} + hf(t_{n},y_{n})\] and 4th order Runge-Kutta \[y_{n+1} = y_{n} + \frac{h}{6}(k_{1} + 2 k_{2}+2 k_{3}+ k_{4})\,\,\,;h = x_{n+1}-x_{n}\\k_{1}=f(x_{n},y_{n})\\k_{2}=f(x_{n}+\frac{h}{2}, y_{n}+\frac{hk_{1}}{2})\\k_{3}=f(x_{n}+\frac{h}{2},y_{n}+\frac{hk_{2}}{2})\\k_{4}=f(x_{n}+h,y_{n}+hk_{3})\]

The implementation relies on the functional aspect of the Scala language and should be flexible enough to support any new future algorithm. The generic Runge-Kutta coefficients a(i), b(i) and c(i) are represented as a matrix: \[\begin{vmatrix} c_{2}\,\,a_{21}\,\,0.0\,\,0.0\,\,...\,\,0.0\,\,0.0\\ c_{3}\,\,a_{31}\,\,a_{32}\,\,0.0\,\,...\,\,0.0\,\,0.0 \\ c_{4}\,\,a_{41}\,\,a_{42}\,\,a_{43}\,\,...\,\,0.0\,\,0.0\\ \\ c_{i}\,\,a_{i1}\,\,a_{i2}\,\,a_{i3}\,\,...\,\,\,\,...\,\,a_{ii-1}\\*\,\,b_{1}\,\,b_{2}\,\,b_{3}\,\,...\,\,...\,\,b_{i-1}\,\,b_{i} \end{vmatrix}\] In order to illustrate the flexibility of this implementation using Scala, I encapsulate the matrix of coefficients of the Euler, 3th order Runge-Kutta, 4th order Runge-Kutta and Felhberg methods using enumeration and case classes.

**Note**: For the sake of readability of the implementation of algorithms, all non-essential code such as error checking, comments, exception, validation of class and method arguments, scoping qualifiers or import is omittedCoefficients: Enumerators and case classes

Java developers, are familiar with enumerators as a data structure to list values without the need to instantiate an iterable collection.

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 | trait RungeKuttaCoefs { type COEFS = Array[Array[Double]] val EULER = Array( Array[Double](0.0, 1.0) ) // Coefficients for Runge-Kutta of order 3 val RK3 = Array( Array[Double](0.0, 0.0, 1/3, 0.0, 0,0), Array[Double](0.5, 0.5, 0.0, 2/3, 0.0), Array[Double](1.0, 0.0, -1.0, 0.0, 1/3) ) // Coefficients for Runge-Kutta of order 4 val RK4 = Array( Array[Double](0.0, 0.0, 1/6, 0.0, 0,0, 0.0), Array[Double](0.5, 0.5, 0.0, 1/3, 0.0, 0.0 ), Array[Double](0.5, 0.0, 0.5, 0.0, 1/3, 0.0), Array[Double](1.0, 0.0, 0.0, 1.0, 0.0, 1/6) ) // Coefficients for Runge-Kutta/Felberg of order 5 val FELBERG = Array( Array[Double](0.0, 0.0, 25/216, 0.0, 0.0, 0.0, 0.0, 0.0), Array[Double](0.25, 0.25, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 ), Array[Double](3/8, 3/32, 0.0, 0.0, 1408/2565, 0.0, 0.0, 0.0), Array[Double](12/13,1932/2197,-7200/2197, 7296/2197, 0.0,2197/4101,0.0,0.0), Array[Double](1.0, 439/216, -8.0, 3680/513, -845/4104, 0.0, -1/5, 0.0), Array[Double](0.5, -8/27, 2.0, -3544/2565, 1859/4104, -11/40, 0.0, 0.0) ) val rk = List[COEFS](EULER, RK3, RK4, FELBERG) } object RungeKuttaForms extends Enumeration with RungeKuttaCoefs{ type RungeKuttaForms = Value val Euler, Rk3, Rk4, Fehlberg = Value @inline final def getRk(value: Value): COEFS = rk(value.id) } |

The first step is to encapsulate the coefficients of the various versions of the Runge-Kutta formuler (lines 4 & 5), Runge-Kutta order 3 (lines 9 - 12), Runge-Kutta order 4 (lines 16 - 20) and Runge-Kutta-Felberg order 5 (lines 24 - 30) in an Scala enumerator.

The enumerator is at best not elegant. As you browse throught the code snippet above, it is clear that the design to wrap the matrices of coefficients with the enumerator is quite cumbersome. There is a better way:

The enumerator is at best not elegant. As you browse throught the code snippet above, it is clear that the design to wrap the matrices of coefficients with the enumerator is quite cumbersome. There is a better way:

**pattern matching**Case classes could be used instead of the singleton enumeration. Setters or getters can optionally be added as in the example below.*The validation of the arguments of methods, exceptions and auxiliary method or variables are omitted for the sake of simplicity.*1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 | trait RKMethods { type COEFS = Array[Array[Double]] def getRk(i: Int, j: Int): COEFS object class Euler extends RKMethods{ Array(Array[Double](0.0, 1.0)) override def getRk(i: Int, j: Int): COEFS {} } object class RK3 extends RKMethods { Array( Array[Double](0.0, 0.0, 1/3, 0.0, 0,0), Array[Double](0.5, 0.5, 0.0, 2/3, 0.0), Array[Double](1.0, 0.0, -1.0, 0.0, 1/3) ) override def getRk(i: Int, j: Int): COEFS {} } .... |

The first step is to encapsulate the coefficients of the various versions of the Runge-Kutta formuler (lines 4 & 5), Runge-Kutta order 3 (lines 9 - 12), Runge-Kutta order 4 (lines 16 - 20) and Runge-Kutta-Felberg order 5 (lines 24 - 30) in an Scala enumerator.

In this second approach, the values of the enumerator are replaced y

In this second approach, the values of the enumerator are replaced y

*Euler*object (line 5),*RK3*- Runge-Kutta order 3 (line 10).Integration

The main class

*RungeKutta*implements all the components necessary to resolve the ordinary differential equation. This simplified implementation relies on adjustable step for the integration.1 2 3 4 5 6 7 8 9 10 11 12 | class RungeKutta( rungeKutta: RungeKuttaForm, adjustStep: (Double, AdjustParameters) => (Double), adjustParameters: AdjustParameters) { final class StepIntegration(val coefs: Array[Array[Double]]) {} def solve( xBegin: Double, xEnd: Double, derivative: (Double, Double) => Double): Double } |

The class

**RungeKutta**has three arguments**rungeKutta**form or type of the Runge-Kutta formula (line 2)**adjustStep**Metric function to adjust the integration step, dx (line 3)<.li>**adjustParameters**Parameters used to compute the derivative (line 4)

1 2 3 4 5 6 7 8 | case class AdjustParameters( maxDerivativeValue: Double = 0.01, minDerivativeValue: Double = 0.00001, gamma: Double = 1.0) { lazy val dx0 = 2.0*gamma/(maxDerivativeValue + minDerivativeValue) } |

The sum of the previous Ks value is computed through an inner loop. The outer loop computes all the values for k using the Runge-Kutta matrix coefficients for this particular method. The integration step is implemented as a tail recursion (lines 14 - 22). but an iterative methods using

*foldLeft*can also be used. The tail recursion may not be as effective in this case because it is implemented as a closure: the method has to close on**ks**.

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 | final class StepIntegration(val coefs : Array[Array[Double]] ) { // Main routine def compute( x: Double, y: Double, dx: Double, derivative: (Double, Double) => Double): Double = { val ks = new Array[Double](coefs.length) // Tail recursion closure @scala.annotation.tailrec def compute(i: Int, k: Double, sum: Double): Double= { ks(i) = k val sumKs= (0 until i)./:(0.0)((s, j) => s + ks(j)*coefs(i)(j+1)) val newK = derivative(x + coefs(i)(0)*dx, y + sumKs*dx) if( i >= coefs.size) sum + newK*coefs(i)(i+2) else compute(i+1, newK, sum + newK*coefs(i)(i+2)) } dx*compute(0, 0.0, 0.0) } |

The next method implements the generic solver that iterates through the entire integration interval. As a matter of fact the solver is indeed implemented as a tail recursion (lines 8 - 20).

## Solver

The accuracy of the solver depends on the value of the increment value**dx**as computed on line 17. We need to weight the accuracy provided by infinitesimal increment with its computation cost. Ideally an adaptive algorithm that compute the value

**dx**according the value

**dy/dx**or

**delta**would provide a good compromise between accuracy and cost. The recursion ends when the recursion value

**x**reaches the end of the integration interval (line 14)..

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 | def solve( xBegin: Double, xEnd: Double, derivative: (Double, Double) => Double): Double ={ val rungeKutta = new StepIntegration(rungeKuttaForm) @scala.annotation.tailrec def solve( x: Double, y: Double, dx: Double, sum: Double): Double = { val z = rungeKutta.compute(x, y, dx, derivative) if( x >= xEnd) sum + z else { val dx = adjustStep(z - y, adjustParameters) solve(x + dx, z, dx, sum+z) } } solve(xBegin, 0.0, adjustParameters.initial, 0.0) } |

The invocation of the solver is very straight forward and can be verified against the analytical solution.

The first step is to define the function that adjusts the integration step (lines 1, 10). This implementation uses the default adjust parameters (line 19) in the initialization of the solver (lines 16 - 19).

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 | val adjustingStep = (diff: Double, adjustParams: AdjustParameters) => { val dx = Math.abs(diff)*adjustParams.dx0/adjustParams.gamma if( dx < adjustParams.minDerivativeValue) adjustParams.minDerivativeValue else if ( dx > adjustParams.maxDerivativeValue) adjustParams.maxDerivativeValue else dx } final val x0 = 0.0 final val xEnd = 2.0 val solver = new RungeKutta( Rk4, adjustingStep, AdjustParameters()) solver.solve( x0, xEnd, (x: Double, y: Double) => Math.exp(-x*x)) |

The family of

*explicit Runge-Kutta methods*provides a good starting point to resolve ordinary differential equations. The current implementation could and possibly should be extended to support adaptive**dx**managed by a control loop using a reinforcement learning algorithm of a Kalman filter of just simple exponential moving average.References

*The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods - E. Hairer, C Lubich, M. Roche - Springer - 1989**Programming in Scala - M. Odersky, L. Spoon, B. Venners - Artima Press 2008**https://github.com/prnicolas*