interfaces/xmaxima/Tkmaxima/Rk.tcl

   1 ############################################################
   2 # Rk.tcl                                                #
   3 #  Copyright William F. Schelter.  All rights reserved.    #
   4 #                                                          #
   5 #     Modified by Jaime E. Villate                         #
   6 #     Time-stamp: "2024-03-27 16:06:40 villate"            #
   7 ############################################################
   8
   9 proc rungeKutta { f g t0 x0 y0 h nsteps} {
  10     set n $nsteps
  11     set ans "$t0 $x0 $y0"
  12     set xn $x0
  13     set yn $y0
  14     set tn $t0
  15     set h2 [expr {$h / 2.0 }]
  16     set h6 [expr {$h / 6.0 }]
  17     catch {
  18         while { [incr nsteps -1] >= 0 } {
  19             set kn1 [$f $tn $xn $yn]
  20             set ln1 [$g $tn $xn $yn]
  21             set arg [list [expr {$tn+$h2}] [expr {$xn+$h2*$kn1}] \
  22                          [expr {$yn+$h2*$ln1}]]
  23             set kn2 [eval $f $arg]
  24             set ln2 [eval $g $arg]
  25             set arg [list [expr {$tn+$h2}] [expr {$xn+$h2*$kn2}] \
  26                          [expr {$yn+$h2*$ln2}]]
  27             set kn3 [eval $f $arg]
  28             set ln3 [eval $g $arg]
  29             set arg [list [expr {$tn+$h}] [expr {$xn+$h*$kn3}] \
  30                          [expr {$yn+$h*$ln3}]]
  31             set kn4 [eval $f $arg]
  32             set ln4 [eval $g $arg]
  33             set xn [expr {$xn+$h6*($kn1+2*$kn2+2*$kn3+$kn4)}]
  34             set yn [expr {$yn+$h6 * ($ln1+2*$ln2+2*$ln3+$ln4)}]
  35             set tn [expr {$tn+ $h}]
  36             lappend ans $tn $xn $yn}}
  37     return $ans}
  38
  39 proc pathLength {list} {
  40     set sum 0
  41     foreach {t x y} $list {
  42         set sum [expr {$sum + sqrt($x*$x+$y*$y)}]
  43     }
  44     return $sum}
  45
  46 proc rungeKuttaA {f g t0 x0 y0 h nsteps} {
  47     set ans [rungeKutta $f $g $t0 $x0 $y0 $h $nsteps]
  48     set count 0
  49     # puts "retrying([llength $ans]) .."
  50     while { [llength $ans] < $nsteps*.5  && $count < 7 } {
  51         incr count
  52         #set leng [pathLength $ans]
  53         #if { $leng == 0 } {set leng .001}
  54         set th [expr {$h/3.0}]
  55         if { $th  < $h }  { set h $th }
  56         set ans  [rungeKutta $f $g $t0 $x0 $y0 $h $nsteps]
  57         # puts -nonewline "..(h=[format "%.5f" $h],pts=[llength $ans])"
  58         # flush stdout
  59     }
  60     return $ans}
  61
  62 ## endsource rk.tcl