interfaces/xmaxima/Tkmaxima/rk4.tcl

   1 #  $Id: rk4.tcl,v 1.5 2009-11-15 23:00:37 villate Exp $
   2 #
   3 # Copyright (C) 2008 Jaime E. Villate <villate@fe.up.pt>
   4 #
   5 # Fourth order Runge-Kutta method, with fixed-lenght steps in phase space.
   6 # Based on Rk.tcl by William F. Schelter
   7 # (the license information can be found in COPYING.tcl)
   8
   9 proc fieldlines { f g t0 x0 y0 sx sy nsteps dir} {
  10     set n $nsteps
  11     set ans "$t0 $x0 $y0"
  12     set xn $x0
  13     set yn $y0
  14     set tn $t0
  15     catch {
  16         while { [incr nsteps -1] >= 0 } {
  17             set kn1 [$f $tn $xn $yn]
  18             set ln1 [$g $tn $xn $yn]
  19
  20             set h [expr {$dir/[vectorlength [expr {$kn1/$sx}] [expr {$ln1/$sy}]]}]
  21             set h2 [expr {$h / 2.0 }]
  22             set h6 [expr {$h / 6.0 }]
  23
  24             set arg [list [expr {$tn + $h2}] [expr {$xn + $h2 * $kn1}] [expr {$yn + $h2*$ln1}]]
  25             set kn2 [eval $f $arg]
  26             set ln2 [eval $g $arg]
  27
  28             set arg [list [expr {$tn + $h2}] [expr {$xn + $h2 * $kn2}] [expr {$yn +$h2*$ln2}]]
  29             set kn3 [eval $f $arg]
  30             set ln3 [eval $g $arg]
  31
  32             set arg [list [expr {$tn + $h}] [expr {$xn + $h * $kn3}] [expr {$yn + $h*$ln3}]]
  33             set kn4 [eval $f $arg]
  34             set ln4 [eval $g $arg]
  35
  36             set dx [expr {$h6 * ($kn1+2*$kn2+2*$kn3+$kn4)}]
  37             set dy [expr {$h6 * ($ln1+2*$ln2+2*$ln3+$ln4)}]
  38
  39             if { [vectorlength $dx $dy] > 5*[vectorlength $sx $sy] } { return $ans }
  40             set xn [expr {$xn + $dx}]
  41             set yn [expr {$yn + $dy}]
  42             set tn [expr {$tn + $h}]
  43
  44             lappend ans $tn $xn $yn
  45         }
  46     }
  47     return $ans
  48 }
  49
  50 proc curves { f g t0 x0 y0 sx sy nsteps dir} {
  51     set n $nsteps
  52     set ans "$t0 $x0 $y0"
  53     set xn $x0
  54     set yn $y0
  55     set tn $t0
  56     catch {
  57         while { [incr nsteps -1] >= 0 } {
  58             set kn1 [expr -1*[$g $tn $xn $yn] ]
  59             set ln1 [$f $tn $xn $yn]
  60
  61             set h [expr {$dir/[vectorlength [expr {$kn1/$sx}] [expr {$ln1/$sy}]]}]
  62             set h2 [expr {$h / 2.0 }]
  63             set h6 [expr {$h / 6.0 }]
  64
  65             set arg [list [expr {$tn + $h2}] [expr {$xn + $h2 * $kn1}] [expr {$yn + $h2*$ln1}]]
  66             set kn2 [expr -1*[eval $g $arg] ]
  67             set ln2 [eval $f $arg]
  68
  69             set arg [list [expr {$tn + $h2}] [expr {$xn + $h2 * $kn2}] [expr {$yn +$h2*$ln2}]]
  70             set kn3 [expr -1*[eval $g $arg] ]
  71             set ln3 [eval $f $arg]
  72
  73             set arg [list [expr {$tn + $h}] [expr {$xn + $h * $kn3}] [expr {$yn + $h*$ln3}]]
  74             set kn4 [expr -1*[eval $g $arg] ]
  75             set ln4 [eval $f $arg]
  76
  77             set dx [expr {$h6 * ($kn1+2*$kn2+2*$kn3+$kn4)}]
  78             set dy [expr {$h6 * ($ln1+2*$ln2+2*$ln3+$ln4)}]
  79
  80             if { [vectorlength $dx $dy] > 5*[vectorlength $sx $sy] } { return $ans }
  81             set xn [expr {$xn + $dx}]
  82             set yn [expr {$yn + $dy}]
  83             set tn [expr {$tn + $h}]
  84
  85             lappend ans $tn $xn $yn
  86         }
  87     }
  88     return $ans
  89 }