Arjen Markus (5 march 2009) Yesterday I spoke with Kevin Kenny and Cameron Laird about the possibilities of developing an application a la XPPAUT [L1 ]. 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: [L2 ]
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 \
"<html>
<head>
<title>Result of Bernoulli computation</title>
</head>
<body>"
set par <p>
set break <br>
set pre <pre>
set nopre </pre>
}
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 "</body>\n</html>"
}
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 <Key-Return> 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