Bernoulli

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