[Tcllib-bugs] [ tcllib-Bugs-1078215 ] Possible bug in Runge-Kutta method

[SourceForge.net <no...@so...>]

Bugs item #1078215, was opened at 2004-12-03 11:21
Message generated for change (Tracker Item Submitted) made by Item Submitter
You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=112883&aid=1078215&group_id=12883

Category: math :: calculus
Group: None
Status: Open
Resolution: None
Priority: 5
Submitted By: Arjen Markus (arjenmarkus)
Assigned to: Arjen Markus (arjenmarkus)
Summary: Possible bug in Runge-Kutta method

Initial Comment:
Mark Stucky reported this in c.l.t.:

I think that I've found a small bug / typo within
tcllib 
(tcllib::math::calculus)
rungeKuttaStep routine.

Possibly I am just missing something here (it wouldn't
be the first time), 
but
here is the routine (from tcllib-1.7) with my suggested
fix :

 proc ::math::calculus::rungeKuttaStep { t tstep xvec
func } {

   set funcq    [uplevel 1 namespace which -command
$func]

   #
   # Four steps:
   # - k1 = tstep*func(t,x0)
   # - k2 = tstep*func(t+0.5*tstep,x0+0.5*k1)
   # - k3 = tstep*func(t+0.5*tstep,x0+0.5*k2)
   # - k4 = tstep*func(t+    tstep,x0+    k3)
   # - x1 = x0 + (k1+2*k2+2*k3+k4)/6
   #
   set tstep2   [expr {$tstep/2.0}]
   set tstep6   [expr {$tstep/6.0}]

   set xk1      [$funcq $t $xvec]
   set xvec2    {}
   foreach x1 $xvec xv $xk1 {
      lappend xvec2 [expr {$x1+$tstep2*$xv}]
   }

   set xk2      [$funcq [expr {$t+$tstep2}] $xvec2]
   set xvec3    {}
   foreach x1 $xvec xv $xk2 {
      lappend xvec3 [expr {$x1+$tstep2*$xv}]
   }

   set xk3      [$funcq [expr {$t+$tstep2}] $xvec3]
   set xvec4    {}
   foreach x1 $xvec xv $xk3 {
      lappend xvec4 [expr {$x1+$tstep2*$xv}]

>********************************************************************************
>  Here is the "bug" : $tstep2 in the above line should
actually be $tstep
>  From the description above
>           # - k4 = tstep*func(t+    tstep,x0+    k3)
>  tstep is used rather than tstep*0.5 (ie tstep2)
>********************************************************************************

   }

   set xk4      [$funcq [expr {$t+$tstep}] $xvec4]
   set result   {}
   foreach x0 $xvec k1 $xk1 k2 $xk2 k3 $xk3 k4 $xk4 {
      set dx [expr {$k1+2.0*$k2+2.0*$k3+$k4}]
      lappend result [expr {$x0+$dx*$tstep6}]
   }

   return $result
}


--Mark 



----------------------------------------------------------------------

You can respond by visiting: 
https://sourceforge.net/tracker/?func=detail&atid=112883&aid=1078215&group_id=12883