[Arjen Markus] (5 march 2009) Yesterday I spoke with [Kevin Kenny] and [Cameron Laird] about the possibilities of developing an application a la XPPAUT [http://www.math.pitt.edu/~bard/xpp/xpp.html]. The idea is to use such packages as [Tclode] and [Plotchart] to solve systems of ordinary differential equations and plot the results. These equations should be definable in a very similar way as XPPAUT allows. The user-interface should be much more modern though. Anyway, here is a (second) shot at something like that. I have called it "Bernoulli" in honour of those great 17th-century mathematicians, Johan, Jacob and Daniel (if not either of the Nicolaus-s (Nikolas, Nicolas)). (It is also a bit of nostalgia for me, as some 15-20 years ago I developed a similar program, also called bernoulli). The version below defines a number of commands and can actually integrate systems of ODEs using the Tclode [extension]. No neat user-interface yet, just the basics. See the example at the end for how to use the commands. (The plotnow command simply plots one variable as function of the other after a computation has been done. Each time a new plot is created) [AM] Update: now with a [GUI]. Save the files with names as given in the comments. ---- [GWM] I have had the temerity to make my own version using the math::calculus to solve a system of ODEs [Bernoulli using math::calculus]. [AM] Tclode is on SourceForge: [http://tclode.sourceforge.net] [AM] I should create a [starkit] for this (and merge Geoffrey's version as it has some interesting extensions ;)) ---- ====== # bernoulli.tcl -- # Integrate systems of ordinary differential equations # # demo -- # Set this variable to 1, to run a simple demo # set demo 0 package require Tk package require Plotchart if {1} { set dir tclode1.0 source $dir/pkgIndex.tcl } package require Tclode # namespaces -- # Create a number of private namespaces # namespace eval ::time { variable steps 100 variable tmin 0.0 variable tmax 1.0 variable statevars {t} } namespace eval ::series {} namespace eval ::init {} namespace eval ::deriv {} # dummy routines -- # Dummy routines in case you want to use bernoulli.tcl # without a GUI or another GUI than the default # proc InitVariable {varname value} {} proc DerivVariable {varname expression} {} # plotf -- # Convenience procedure to plot a function of one variable # # Arguments: # range Range over which to plot the function (x-axis) # expression Expression by which the function is defined # (independent variable: t) # # Result: # None # proc plotf {range expression} { # Very simple for now: no parameters foreach {tmin tmax} $range {break} set nosteps 100 set tstep [expr {($tmax-$tmin)/double($nosteps)}] if { ![winfo exists .c] } { pack [canvas .c] } else { .c delete all } set tvalues {} set xvalues {} set xmin {} set xmax {} for {set i 0} {$i < $nosteps+1} {incr i} { set t [expr {$tmin + $i*$tstep}] set x [expr $expression] lappend tvalues $t lappend xvalues $x if { $xmin == {} || $x < $xmin } { set xmin $x } if { $xmax == {} || $x > $xmax } { set xmax $x } } set tscale [::Plotchart::determineScale $tmin $tmax] set xscale [::Plotchart::determineScale $xmin $xmax] set p [::Plotchart::createXYPlot .c $tscale $xscale] foreach t $tvalues x $xvalues { $p plot function $t $x } set ::time::plot $p } # plotnow -- # Convenience procedure to plot one variable as function of the other # # Arguments: # xvar Variable on the x-axis # yvar Variable on the y-axis # # Result: # None # proc plotnow {xvar yvar} { if { ![winfo exists .c] } { pack [canvas .c] } else { .c delete all } set xvalues [set ::series::$xvar] set yvalues [set ::series::$yvar] set xmin {} set xmax {} set ymin {} set ymax {} foreach x $xvalues y $yvalues { if { $xmin == {} || $x < $xmin } { set xmin $x } if { $xmax == {} || $x > $xmax } { set xmax $x } if { $ymin == {} || $y < $ymin } { set ymin $y } if { $ymax == {} || $y > $ymax } { set ymax $y } } set xscale [::Plotchart::determineScale $xmin $xmax] set yscale [::Plotchart::determineScale $ymin $ymax] set p [::Plotchart::createXYPlot .c $xscale $yscale] foreach x $xvalues y $yvalues { $p plot function $x $y } set ::time::plot $p } # timeframe -- # Define the time frame for the next computation # # Arguments: # tmin Minimum time # tmax Maximum time # # Result: # None # proc timeframe {tmin tmax} { set time::tmin $tmin set time::tmax $tmax } # init -- # Define the initial value of a variable # # Arguments: # varname Name of the (state) variable # value Initial value # # Result: # None # proc init {varname value} { set init::$varname $value InitVariable $varname $value } # deriv -- # Register the expression that determines the derivative # # Arguments: # varname Name of the (state) variable # expression Expression to be used # # Result: # None # proc deriv {varname expression} { if { [lsearch $::time::statevars $varname] < 0 } { lappend ::time::statevars $varname } set deriv::$varname $expression DerivVariable $varname $expression } # steps -- # Set the number of steps to report/plot # # Arguments: # nosteps Number of steps # # Result: # None # proc steps {nosteps} { set time::steps $nosteps } # go -- # Do the computation # # Arguments: # None # # Result: # None # proc go {} { # # Time administration and initial values # set t $::time::tmin set ::time::dt [expr {($::time::tmax-$t)/double($::time::steps)}] set args {} set ::series::t {} foreach v [lrange $::time::statevars 1 end] { lappend args $v [set ::deriv::$v] set $v [set ::init::$v] set ::series::$v {} } # # Construct the ODE solver # set ::time::solver [::tclode::odesolv -atol 1.0e-6 -rtol 1.0e-6 $args] # # Run the computation # $::time::solver run $t for {set ::time::i 0} {$::time::i <= $::time::steps} {incr ::time::i} { set ::time::tnext [expr {$::time::tmin + $::time::dt * $::time::i}] $::time::solver continue $::time::tnext foreach v $::time::statevars { lappend ::series::$v [set $v] } } } # reportResult -- # Write the result to an external file # # Arguments: # filename Name of the file to write to # # Result: # None # # Note: # The type of file is determined from the extension # proc reportResult {filename} { set type [lsearch {.txt .html .tsv .csv} [string tolower [file extension $filename]]] if { $type < 0 } { return -code error "Unknown file type: [file extension $filename]" } set outfile [open $filename w] set par "" set pre "" set nopre "" set break "" set sep " " switch -- $type { 0 { puts $outfile "Result of Bernoulli computation\n" } 1 { puts $outfile \ " Result of Bernoulli computation " set par

set break
set pre

            set nopre 
} 2 { puts $outfile "Result of Bernoulli computation\n" set sep \t } 3 { puts $outfile "Result of Bernoulli computation\n" set sep , } default { return -code error "Unknown file type: [file extension $filename]" } } puts $outfile "Differential equations:$par" foreach v [lrange $::time::statevars 1 end] { set text [list $v: [set ::deriv::$v]] puts $outfile [join $text $sep]$break } puts $outfile "Initial conditions:$par" foreach v [lrange $::time::statevars 1 end] { set text [list $v = [set ::init::$v]] puts $outfile [join $text $sep]$break } puts $outfile "Results:$par$pre" set text "" foreach v $::time::statevars { lappend text [format "%12.12s" $v] } puts $outfile [join $text $sep] set i 0 set n [llength $::series::t] for { set i 0} {$i < $n} {incr i} { set text {} foreach v $::time::statevars { lappend text [format "%12.4g" [lindex [set ::series::$v] $i]] } puts $outfile [join $text $sep] } puts $outfile $nopre if { $type == 1 } { puts $outfile "\n" } close $outfile } # main -- # Test the stuff # #plotf {0 10} {exp(-$t)*cos(2.0*$t)} #after 3000 { # plotf {0.0001 100} {sin($t)/$t} ;# Avoid t=0 of course #} if {$demo} { init s 0.0 init c 1.0 init x 0.0 deriv s {$c} deriv c {-$s} deriv x {1} timeframe 0.0 10.0 go catch { console show } foreach t $::series::t s $::series::s c $::series::c x $::series::x { puts [format "%10.4f %10.4f %10.4f %10.4f" $t $s $c $x] } plotnow t c } ====== The GUI is in this file: ====== # bern_gui.tcl -- # Integrate systems of ordinary differential equations: # GUI part # namespace eval ::gui { variable command "" variable version 0.1 variable date "March, 2009" } source [file join [file dirname [info script]] bernoulli.tcl] # ShowInfo -- # Show information on how to use Bernoulli # # Arguments: # None # # Result: # None # proc ShowInfo {} { if { [winfo exists .help] } { wm raise .help return } toplevel .help frame .help.f text .help.f.info -yscrollcommand {.help.f.y set} ::ttk::scrollbar .help.f.y -orient vertical -command {.help.f.info yview} button .help.ok -text OK -width 8 -command {destroy .help} grid .help.f.info .help.f.y -sticky news grid .help.f -sticky news grid .help.ok .help.f.info insert end \ "BERNOULLI - quick reference init varname value Defines a state variable and it's initial value deriv varname expression Defines the differential equation by which to compute the derivative. Use braces! steps number Set number of output steps timeframe begin end Set begin and end for the variable t (time) plotnow xvar yvar Plot the variable yvar versus xvar in the plot plotf {begin end} expression Plot an expression in t for t from 'begin' to 'end' go Do the computation (after this the results are available) " .help.f.info see end } # ShowAbout -- # Show an About box # # Arguments: # None # # Result: # None # proc ShowAbout {} { tk_dialog .about "About Bernoulli" \ "Bernoulli: Integration of systems of ordinary differential equations Version $::gui::version $::gui::date" "" 0 " OK " } # InitVariable -- # Store the variable and the initial value in the tree # # Arguments: # varname Name of the variable # value Initial value # # Result: # None # proc InitVariable {varname value} { if { ![info exists ::gui::var($varname,init)] } { set ::gui::var($varname,init) [.tree insert {} end -text $varname] } set node $::gui::var($varname,init) .tree item $node -values [list "Initial: $value"] } # DerivVariable -- # Store the variable and the expression for its derivative in the tree # # Arguments: # varname Name of the variable # expression Expression for derivative # # Result: # None # proc DerivVariable {varname expression} { if { ![info exists ::gui::var($varname,init)] } { set ::gui::var($varname,init) [.tree insert {} end -text $varname] } set node $::gui::var($varname,init) if { ![info exists ::gui::var($varname,deriv)] } { set ::gui::var($varname,deriv) [.tree insert $node end -text "Derivative:"] } set node $::gui::var($varname,deriv) .tree item $node -values [list "$expression"] } # AddMessage -- # Add a message to the history area # # Arguments: # msg Message to show # # Result: # None # proc AddMessage {msg} { .cmd.history.text insert end $msg\n message .cmd.history.text see end } # OpenDefFile -- # Open a system definition file # # Arguments: # None # # Result: # None # # Side effects: # System definition file loaded # proc OpenDefFile {} { set types { {{Definition Files} {.def}} {{All Files} * }} set defFile [tk_getOpenFile -filetypes $types] if { $defFile != "" } { foreach child [.tree children {}] { .tree delete $child } array unset ::gui::var source $defFile set ::gui::defFile $defFile AddMessage "Opened $defFile successfully" } } # ReportResult -- # Write the results to a file # # Arguments: # None # # Result: # None # # Side effects: # System definition file loaded # proc ReportResult {} { set types { {{Plain text Files} {.txt}} {{HTML Files} {.html}} {{TSV Files} {.tsv}} {{CSV Files} {.csv}}} set reportFile [tk_getSaveFile -filetypes $types -defaultextension ".txt"] if { $reportFile != "" } { reportResult $reportFile } AddMessage "Saved report in $reportFile" } # SavePlot -- # Save the current plot # # Arguments: # None # # Result: # None # proc SavePlot {} { set types { {{PostScript Files} {.eps}}} set plotFile [tk_getSaveFile -filetypes $types -defaultextension ".eps"] if { $plotFile != "" } { $::time::plot saveplot $plotFile } AddMessage "Saved plot in $plotFile" } # RunCommand -- # Run the command that was entered # # Arguments: # None # # Result: # None # # Side effects: # Whatever the command was # proc RunCommand {} { # TODO: Should probably use a safe interpreter ... eval $::gui::command .cmd.history.text insert end $::gui::command\n .cmd.history.text see end set ::gui::command {} } # CreateMenu -- # Set up the menu bar and the menus # # Arguments: # None # # Result: # None # # Side effects: # Unmanaged menubar created (.menu) # proc CreateMenu {} { if {[tk windowingsystem] eq "aqua"} { set modifier Command } elseif {$::tcl_platform(platform) == "windows"} { set modifier Control } else { set modifier Meta } set menu ".menu" # # Menu bar # menu $menu -tearoff 0 # # File menu # set f $menu.file menu $f -tearoff 0 $menu add cascade -label File -menu $f -underline 0 $f add command -label "Open..." -command OpenDefFile $f add separator $f add command -label "Save" -command {SaveDefFile 0} -accelerator $modifier+S $f add command -label "Save As ..." -command {SaveDefFile 1} $f add separator $f add command -label "Report ..." -command ReportResult $f add command -label "Save plot ..." -command SavePlot $f add separator $f add command -label "Exit" -command exit # # Help menu # set h $menu.help menu $h -tearoff 0 $menu add cascade -label Help -menu $h -underline 0 $h add command -label "Information" -command ShowInfo $h add separator $h add command -label "About ..." -command ShowAbout . configure -menu $menu } # CreateCommandArea -- # Create an area for entering commands # # Arguments: # None # # Result: # None # # Side effects: # Unmanaged command area (entry field + text widget) created (.cmd) # proc CreateCommandArea {} { frame .cmd frame .cmd.history text .cmd.history.text -height 5 \ -yscrollcommand {.cmd.history.y set} scrollbar .cmd.history.y -orient vertical \ -command {.cmd.history.text yview} grid .cmd.history.text .cmd.history.y -sticky news grid rowconfigure .cmd.history 0 -weight 1 grid columnconfigure .cmd.history 0 -weight 1 label .cmd.label -text "Run command:" entry .cmd.edit -textvariable ::gui::command button .cmd.go -text Run -command RunCommand grid .cmd.history - - -sticky news grid .cmd.label - - -sticky w grid .cmd.edit .cmd.go -sticky news grid columnconfigure .cmd 0 -weight 1 bind .cmd.edit RunCommand .cmd.history.text tag configure message -foreground blue } # MainWindow -- # Set up the main window # # Arguments: # None # # Result: # None # # Side effects: # Main window set up: # - Menu bar: load a system definition, save the results # - Left: tree view of the system definition # - Middle: Plot window # - Below: command area # proc MainWindow {} { CreateMenu CreateCommandArea canvas .c -width 400 -height 400 ::ttk::treeview .tree -columns {value} .tree insert {} end -text "No model definition ..." .tree heading #0 -text Variable .tree heading value -text Value grid .c .tree -sticky news grid .cmd ^ -sticky news grid rowconfigure . 0 -weight 1 grid columnconfigure . 1 -weight 1 .c create text 10 40 -text "Load a model definition or define one interactively" -anchor w } # main -- # Bring up the main window and get started # MainWindow ====== And here is a simple example: ====== # circle.def -- # Simple example of a system definition file for Bernoulli # init s 0.0 init c 1.0 init x 0.0 deriv s {$c} deriv c {-$s} deriv x {1} timeframe 0.0 10.0 go plotnow t s ====== ---- !!!!!! %| [Category Mathematics] |% !!!!!!