interfaces/xmaxima/Tkmaxima/Adams.tcl

   1 ############################################################
   2 # Adams.tcl                                                #
   3 #  Copyright William F. Schelter.  All rights reserved.    #
   4 #                                                          #
   5 #     Modified by Jaime E. Villate                         #
   6 #     Time-stamp: "2024-03-16 20:45:39 villate"            #
   7 ############################################################
   8
   9 proc adamsMoulton {f g t0 x0 y0 h nsteps} {
  10     set ans [rungeKutta $f $g $t0 $x0 $y0 $h 3]
  11     catch {
  12         set i 0
  13         set h24 [expr {$h/24.0}]
  14         foreach {t x y } $ans {
  15             lappend listXff [xff [expr {$t0+$i*$h}] $x $y]
  16             lappend listYff [yff [expr {$t0+$i*$h}] $x $y]
  17             incr i
  18             set xn $x
  19             set yn $y}
  20         set n [expr $nsteps -3]
  21         while { [incr n -1] >= 0 } {
  22             #puts "listXff = $listXff"
  23             #puts "listYff = $listYff"
  24             # adams - bashford formula:
  25             set xp [expr {$xn + ($h24)*(55 *[lindex $listXff 3]-59*[lindex $listXff 2]+37*[lindex $listXff 1]-9*[lindex $listXff 0]) }]
  26             set yp [expr {$yn + ($h24)*(55 *[lindex $listYff 3]-59*[lindex $listYff 2]+37*[lindex $listYff 1]-9*[lindex $listYff 0]) }]
  27             #puts "i=$i,xp=$xp,yp=$yp"
  28             # adams-moulton corrector-predictor:
  29             # compute the yp = yn+1 value..
  30             set t [expr {$t0+$i*$h}]
  31             incr i
  32             if { 1 } {
  33                 set xap [expr { $xn+($h24)*(9*[xff $t $xp $yp]+19*[lindex $listXff 3]-5*[lindex $listXff 2]+[lindex $listXff 1]) }]
  34                 set yap [expr { $yn+($h24)*(9*[yff $t $xp $yp]+19*[lindex $listYff 3]-5*[lindex $listYff 2]+[lindex $listYff 1]) }]
  35
  36                 set xn $xap
  37                 set yn $yap
  38                 # puts "after correct:i=[expr $i -1],xn=$xn,yn=$yn"
  39                 # could repeat it, or check against previous to see if changes too much.
  40             }
  41             set listXff [lrange $listXff 1 end]
  42             set listYff [lrange $listYff 1 end]
  43
  44             lappend listXff [xff $t $xn $yn]
  45             lappend listYff [yff $t $xn $yn]
  46
  47             lappend ans $t $xn $yn
  48             # puts "ans=$ans"}
  49         #puts "adams:t=$t"
  50     }
  51     return $ans}
  52
  53 ## endsource adams.tcl