注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

阿弥陀佛

街树飘影未见尘 潭月潜水了无声 般若观照心空静...

 
 
 

日志

 
 
关于我

一直从事气象预报、服务建模实践应用。 注重气象物理场、实况场、地理信息、本体知识库、分布式气象内容管理系统建立。 对Barnes客观分析, 小波,计算神经网络、信任传播、贝叶斯推理、专家系统、网络本体语言有一定体会。 一直使用Java、Delphi、Prolog、SQL编程。

网易考拉推荐

Runge-Kutta  

2016-06-14 13:37:16|  分类: Scala |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
Runge-Kutta

Monday, June 10, 2013

Runge-Kutta ODE solver in Scala

Objective
The objective of this post 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
dydx=f(x,y)
The family of explicit Runge-Kutta numerical approximation methods is defined as
yn+1=yn+i=0s<nbikiwherekj=h.f(xn+cjh,yn+s=1j?1as,s?1ks?1)withh=xn+1?xnandΔ=dydx+s=1j?1as,s?1ks?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
yn+1=yn+hf(tn,yn)
and 4th order Runge-Kutta
yn+1=yn+h6(k1+2k2+2k3+k4);h=xn+1?xnk1=f(xn,yn)k2=f(xn+h2,yn+hk12)k3=f(xn+h2,yn+hk22)k4=f(xn+h,yn+hk3)

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:
c2a210.00.0...0.00.0c3a31a320.0...0.00.0c4a41a42a43...0.00.0ciai1ai2ai3......aii?1?b1b2b3......bi?1bi
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 omitted

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

trait RungeKuttaCoefs {
  type COEFS = Array[Array[Double]]
  
  private val EULER = Array(Array[Double](0.0, 1.0))
      // Coefficients for Runge-Kutta of order 3
  private 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 
  private 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
  private 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))
    
  protected 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 enumerator is at best not elegant and the design to hide the matrices of coefficients behind 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.

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 {}
}
....


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.

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 computation of the parameters to adjust the integration step is rather simple. A more elaborate implementation would include several alternative formulas implemented as sealed case class

case class AdjustParameters(
     maxDerivativeValue: Double = 0.01,
     minDerivativeValue: Double = 0.00001,
     gamma: Double = 1.0) {

  lazy val dx0 = 2.0*gamma/(maxDerivativeValue + minDerivativeValue)
}

The implementation of the generic algorithm over the interval [x + dx, x]. 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 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.

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).foldLeft(0.0)((sum, j) => { sum + 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 iterate through the entire integration interval. The accuracy of the solver depends on the value of the increment value dx.  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. Once again, the solve method is implemented as a tail recursion.

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.
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

No comments:

Post a Comment

Subscribe to: Post Comments (Atom)

Contact Form

Name
Email *
Message *
  评论这张
 
阅读(208)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2017