[Arjen Markus] (22 december 2014) The program below implements a
user-interface to set up a system of ordinary differential equations and
solve it. Via the GUI you solve the equations with the initial
conditions you set and the parameter values you selected and
subsequently present them as a simple graph.

The program is not perfect - there are a few quirks that I need to
iron out, mostly to do with the way it manages the state variables and
the contents of the various tab windows: if you switch from one system
to another, it leaves too much of the old stuff. But that is mostly
laziness from me and the desire to write it from start to finish in a
single weekend. I say this, not to boast about my own capabilities to
program fast, but to praise Tcl and Tk and the packages used for making
this kind of productivity possible!

Here is a screenshot of the calculation window.

[Screenshot GUI defining and solving ODEs]

======
# solve_ods.tcl --
#     Program to specify a system of ordinary differential equations and to
#     solve them. Provides a graphical display of the results
#
#     TODO:
#     - Scrollbars
#     - Use more vectcl - for instance: storing the results
#     - Use initialisation code
#     - Help text
#     - Keep track of changes in the definition and parameters since the last save action
#     - Erase old parameters, table and plot when starting with a different ODE file/anew
#
#     Enhancements:
#     - Monte Carlo simulation
#     - Allow customised plotting
#     - Different integration options
#     - More examples
#
#     Note:
#     Issue with tablelist:
#     If you edit an entry and then switch the tab window, the new value is not saved!
#

lappend auto_path c:/tcl/lib  ;# Hack!

package require vectcl
package require Plotchart
package require tablelist

# GUI --
#     Provide a convenient namespace
#
namespace eval GUI {
    variable documentation
    variable variables
    variable filename
    variable plot
    variable helpText

    set documentation(t)         "Time (or independent variable in general)"
    set documentation(tbegin)    "Start time of simulation"
    set documentation(tend)      "End time of simulation"
    set documentation(tstep)     "Time step"
    set documentation(tstepout)  "Time step for output"

    set filename   ""

    set ::t        0.0
    set ::tbegin   0.0
    set ::tend     1.0
    set ::tstep    0.025
    set ::tstepout 0.1

    set helpText "
*Solving systems of ordinary differential equations (ODE)
This program solves systems of ordinary differential equations. You define
the equations in the tab window \"Definition\". Here, a prefix \"d_\" indicates
that the variable is a state variable. All other variables are parameters or
auxiliary variables.

For example:

d_x = -k*y - r*x
d_y = k*x

defines a system of two equations with the state variables x and y. The system
is parameterised with the parameters k, the angular velocity, and r, the
friction coefficient.

*Tab window: Calculation
You set the values for the calculation in the \"Calculation\" tab. You
can also change the description of the variables, so that they are better
recognisable.

Note: there is a quirk in the editing of values - you should not change the
tab window or press \"Run\" without clicking on another value, otherwise the
value is lost.

*Tab window: Result
When you press the \"Run\" button, the calcuation is started and the results
are shown in the \"Result\" tab window. You can store these results in an
external file via the \"File\" menu.

*Tab window: Display
The program allows you to plot the results - either a state variable as
function of t (time) or a state variable as function of another state variable.

*Caveats
There are plenty of quirks to correct still in this version, especially
when it comes to switching from one ODE file to another."
}

# setupMainWindow --
#     Set up the tab windows that constitute the GUI
#
# Arguments:
#     None
#
# Returns:
#     Nothing
#
# Side effects:
#     Create the tab windows "Definition", "Calculation", "Result", "Plot" and "Help"
#
proc ::GUI::setupMainWindow {} {
    variable plot
    #
    # Dedicated toplevel
    #
    toplevel .ode

    #
    # Set up the menu
    #
    menu .ode.menu -tearoff 0
    menu .ode.menu.file -tearoff 0
    .ode.menu add cascade -label "File" -menu .ode.menu.file -underline 0
    .ode.menu.file add command -label "New" -command {::GUI::newOdeFile}
    .ode.menu.file add command -label "Open..." -command {::GUI::openOdeFile 1}
    .ode.menu.file add separator
    .ode.menu.file add command -label "Save" -command {::GUI::saveOdeFile 0}
    .ode.menu.file add command -label "Save as..." -command {::GUI::saveOdeFile 1}
    .ode.menu.file add separator
    .ode.menu.file add command -label "Save result" -command {::GUI::saveResult}
    .ode.menu.file add separator
    .ode.menu.file add command -label "Exit" -command {::GUI::exitProgram}

    .ode configure -menu .ode.menu

    #
    # Set up the tab windows
    #
    ::ttk::notebook .ode.nb

    #
    # Definition tab
    #

    # TODO: scrollbars

    ::ttk::frame .ode.nb.definition
    text .ode.nb.definition.def  -height 10 \
              -xscrollcommand {.ode.nb.definition.xscrolldef set} \
              -yscrollcommand {.ode.nb.definition.yscrolldef set}
    text .ode.nb.definition.init -height 10 \
              -xscrollcommand {.ode.nb.definition.xscrollinit set} \
              -yscrollcommand {.ode.nb.definition.yscrollinit set}
    ::ttk::scrollbar .ode.nb.definition.xscrolldef  -orient horizontal -command {.ode.nb.definition.def xview}
    ::ttk::scrollbar .ode.nb.definition.yscrolldef  -orient vertical   -command {.ode.nb.definition.def yview}
    ::ttk::scrollbar .ode.nb.definition.xscrollinit -orient horizontal -command {.ode.nb.definition.init xview}
    ::ttk::scrollbar .ode.nb.definition.yscrollinit -orient vertical   -command {.ode.nb.definition.init yview}

    grid .ode.nb.definition.def  .ode.nb.definition.yscrolldef  -pady 3 -sticky news
    grid .ode.nb.definition.xscrolldef                          -pady 3 -sticky news
    grid .ode.nb.definition.init .ode.nb.definition.yscrollinit -pady 3 -sticky news

    grid columnconfigure .ode.nb.definition 0 -weight 1
    grid rowconfigure    .ode.nb.definition 0 -weight 1
    grid rowconfigure    .ode.nb.definition 2 -weight 1

    .ode.nb.definition.def  insert end "# This text pane is for defining the differential equations\n# See \"Help\" tab for more information\n"
    .ode.nb.definition.init insert end "# This text pane is for defining auxiliary functions\n\n"

    .ode.nb add .ode.nb.definition -text "Definition" -padding 3 -underline 0

    grid .ode.nb -sticky news
    grid columnconfigure .ode 0 -weight 1
    grid rowconfigure    .ode 0 -weight 1

    #
    # Calculation
    #
    ::ttk::frame .ode.nb.calculation
    ::tablelist::tablelist .ode.nb.calculation.table -width 80 -height 10 -stretch all \
        -columns {0 Variable 0 Value 0 Description} \
        -editendcommand ::GUI::applyValue \
        -yscrollcommand {.ode.nb.calculation.yscroll set}

    ::ttk::scrollbar .ode.nb.calculation.yscroll -orient vertical -command {.ode.nb.calculation.table yview}

    .ode.nb.calculation.table columnconfigure 1 -editable 1 -editwindow ttk::entry
    .ode.nb.calculation.table columnconfigure 2 -editable 1 -editwindow ttk::entry

    ::ttk::button .ode.nb.calculation.run -command {::GUI::runCalculation 1} -text Run
    grid .ode.nb.calculation.table .ode.nb.calculation.yscroll -sticky news
    grid .ode.nb.calculation.run
    grid columnconfigure .ode.nb.calculation 0 -weight 1
    grid rowconfigure    .ode.nb.calculation 0 -weight 1
    .ode.nb add .ode.nb.calculation -text "Calculation" -padding 3 -underline 0

    bind .ode.nb <ButtonPress> {::GUI::checkDefinition [.ode.nb identify tab %x %y]}
    bind .ode.nb.definition.def  <Alt-KeyPress> {::GUI::checkDefinition %K}
    bind .ode.nb.definition.init <Alt-KeyPress> {::GUI::checkDefinition %K}

    #
    # Tabular output
    #
    ::ttk::frame  .ode.nb.result
    ::tablelist::tablelist .ode.nb.result.table -width 80 -height 10 -stretch all -columns {0 t 0 ...} \
         -xscrollcommand {.ode.nb.result.xscroll set} -yscrollcommand {.ode.nb.result.yscroll set}

    ::ttk::scrollbar .ode.nb.result.xscroll -orient horizontal -command {.ode.nb.result.table xview}
    ::ttk::scrollbar .ode.nb.result.yscroll -orient vertical   -command {.ode.nb.result.table yview}

    grid   .ode.nb.result.table   .ode.nb.result.yscroll -sticky news
    grid   .ode.nb.result.xscroll                        -sticky news

    grid columnconfigure .ode.nb.result 0 -weight 1
    grid rowconfigure    .ode.nb.result 0 -weight 1

    .ode.nb add .ode.nb.result -text "Result" -padding 3 -underline 0

    #
    # Graphical output
    #
    ::ttk::frame  .ode.nb.plot
    canvas .ode.nb.plot.cnv -bg white
    ::ttk::label    .ode.nb.plot.xaxislabel -text "X-axis:" -width 10
    ::ttk::combobox .ode.nb.plot.xaxis -state readonly
    ::ttk::label    .ode.nb.plot.empty      -text " "
    ::ttk::label    .ode.nb.plot.yaxislabel -text "Y-axis:" -width 10
    ::ttk::combobox .ode.nb.plot.yaxis -state readonly
    ::ttk::button   .ode.nb.plot.display -text Display -command {::GUI::displayResult}
    grid   .ode.nb.plot.cnv - - - - - -sticky news
    grid   .ode.nb.plot.xaxislabel .ode.nb.plot.xaxis \
           .ode.nb.plot.empty                         \
           .ode.nb.plot.yaxislabel .ode.nb.plot.yaxis \
           .ode.nb.plot.display -sticky w -pady 3

    grid columnconfigure .ode.nb.plot 0 -weight 1
    grid rowconfigure    .ode.nb.plot 0 -weight 1

    .ode.nb add .ode.nb.plot -text "Plot" -padding 3 -underline 0

    set plot [::Plotchart::createXYPlot .ode.nb.plot.cnv {0 1 0.25} {0 1 0.25}]

    #
    # Help text
    #
    ::ttk::frame  .ode.nb.help
    text .ode.nb.help.text -xscrollcommand {.ode.nb.help.xscroll set} \
                           -yscrollcommand {.ode.nb.help.yscroll set}
    ::ttk::scrollbar .ode.nb.help.xscroll -orient horizontal -command {.ode.nb.help.text xview}
    ::ttk::scrollbar .ode.nb.help.yscroll -orient vertical   -command {.ode.nb.help.text yview}
    grid .ode.nb.help.text    .ode.nb.help.yscroll -sticky news
    grid .ode.nb.help.xscroll -                    -sticky news
    .ode.nb add .ode.nb.help -text "Help" -padding 3 -underline 0
    displayHelpText

    grid columnconfigure .ode.nb.help 0 -weight 1
    grid rowconfigure    .ode.nb.help 0 -weight 1
    grid rowconfigure    .ode.nb.help 1 -weight 0

    ttk::notebook::enableTraversal .ode.nb
}

# applyValue --
#     Apply the new value to the variable in the "Calculation" tab
#
# Arguments:
#     widget       Name of the tablelist widget
#     row          Row of the variable that was edited
#     column       Column of the variable that was edited
#     text         New value
#
# Returns:
#     Nothing
#
# Side effects:
#     Updated value of the underlying variable
#
proc ::GUI::applyValue {widget row column text} {
    variable variables
    variable documentation

    set v [lindex $variables $row]
    if { $column == 1 } {
        set ::$v $text
    } else {
        set documentation($v) $text
    }
    return $text
}

# displayHelpText --
#     Display the help text
#
# Arguments:
#     None
#
# Returns:
#     Nothing
#
# Side effects:
#     Help window filled with text - simple formatting
#
proc ::GUI::displayHelpText {} {
    variable helpText

    .ode.nb.help.text tag configure Section -font "Helvetica 12 bold"

    foreach line [split $helpText \n] {
        if { [string index $line 0] == "*" } {
            .ode.nb.help.text insert end "[string range $line 1 end]\n" Section
        } else {
            .ode.nb.help.text insert end "$line\n"
        }
    }

    .ode.nb.help.text configure -state disabled
}

# checkDefinition --
#     Analyse the text in the definition pane for variables and parameters
#
# Arguments:
#     tab          The index of the tab that was pressed
#
# Returns:
#     Nothing
#
# Side effects:
#     Updated list of variables and parameters in the "calculation" tab
#
proc ::GUI::checkDefinition {tab} {
    variable documentation
    variable variables
    variable stateVariables

    if { $tab != 1 && $tab != "c" && $tab != "C" } {
        return
    }

    #
    # Get the text of the definition pane, remove comments and
    # analyse it
    #
    set definition [.ode.nb.definition.def get 1.0 end]
    set cleanedup ""
    foreach line [split $definition \n] {
        regsub {#.*} $line {} line
        append cleanedup "$line\n"
    }
    puts $cleanedup

    set rawVarList    [regexp -all -inline {[_[:alnum:]]+} $cleanedup]
    set deriveVarList [lsearch -all -inline -regexp $cleanedup {^d_.+}]

    set stateVariables {}
    foreach v $deriveVarList {
        lappend stateVariables [string range $v 2 end]
    }

    set parameterList {}
    foreach v $rawVarList {
        if { $v in $stateVariables || $v in $deriveVarList || $v in $parameterList } {
            continue
        }
        if { [string is integer [string index $v 0]] } {
            continue
        }
        lappend parameterList $v
    }

    #
    # Display it
    #
    .ode.nb.calculation.table delete 0 end
    set variables {}
    foreach v {tbegin tend tstep tstepout} {
        lappend variables $v
        .ode.nb.calculation.table insert end [list $v [set ::$v] $documentation($v)]
    }

    foreach v $stateVariables {
        if { ![info exists documentation($v)] } {
            set documentation($v) "State variable $v"
        }
        if { ![info exists ::$v] } {
            set ::$v 0.0
        }
        lappend variables $v
        .ode.nb.calculation.table insert end [list $v [set ::$v] $documentation($v)]
    }

    foreach v $parameterList {
        if { ![info exists documentation($v)] } {
            set documentation($v) "Parameter $v"
        }
        if { ![info exists ::$v] } {
            set ::$v 0.0
        }
        lappend variables $v
        .ode.nb.calculation.table insert end [list $v [set ::$v] $documentation($v)]
    }

    .ode.nb.plot.xaxis configure -values [concat t $stateVariables]
    .ode.nb.plot.yaxis configure -values [concat t $stateVariables]
    .ode.nb.plot.xaxis current 0
    #.ode.nb.plot.yaxis current 1 ;# Not applicable
}

# newOdeFile --
#     Reset everything
#
# Arguments:
#     None
#
# Returns:
#     Nothing
#
# Side effects:
#     Stops the program
#
proc ::GUI::newOdeFile {} {
    variable filename
    variable variables
    variable documentation

    #
    # Force an update of the system
    #
    checkDefinition 1

    if { [llength $variables] > 4 } {
        if { [tk_messageBox -icon question -type yesno -message "Reset the system?"] == "yes" } {
            foreach v [lrange $variables 4 end] {
                unset documentation($v)
            }
            set filename  {}
            set variables {}
            .ode.nb.definition.def  delete 1.0 end
            .ode.nb.definition.init delete 1.0 end
            .ode.nb.calculation.table delete 0 end
        }
    }
}

# saveResult --
#     Save the current results
#
# Arguments:
#     None
#
# Returns:
#     Nothing
#
# Side effects:
#     Saves the results in external file
#
proc ::GUI::saveResult {} {
    variable values
    variable stateVariables
    variable documentation

    #
    # Ask for the file name
    #
    set filename [tk_getSaveFile -defaultextension ".csv" -filetypes {{.csv "CSV files"} {*.* "All files"}} \
                      -title "Save the results" -parent .ode]
    if { $filename != "" } {
        set outfile [open $filename w]
        puts $outfile [join [concat t $stateVariables] \t]
        foreach record $values {
            puts $outfile [join $record \t]
        }
        close $outfile
    }
}

# saveOdeFile --
#     Save the definition and the initial values
#
# Arguments:
#     saveas               Save the data under a new name or not
#
# Returns:
#     Nothing
#
# Side effects:
#     Saves definition and initial values in external file
#
proc ::GUI::saveOdeFile {saveas} {
    variable filename
    variable variables
    variable documentation

    #
    # Force an update of the system
    #
    checkDefinition 1

    #
    # Ask for the file name
    #
    if { $saveas || $filename == "" } {
        set filename [tk_getSaveFile -defaultextension ".ode" -filetypes {{.ode "ODE files"} {*.* "All files"}} \
                          -title "Save the system definition" -initialfile $filename -parent .ode]
    }
    if { $filename != "" } {
        set outfile [open $filename w]
        puts $outfile "set definition     {[.ode.nb.definition.def  get 1.0 end]}"
        puts $outfile "set initialisation {[.ode.nb.definition.init get 1.0 end]}"
        foreach v $variables {
            puts $outfile "set documentation($v) {$documentation($v)}"
            puts $outfile "set ::$v [set ::$v]"
        }
        close $outfile
    }
}

# openOdeFile --
#     Open an existing ODE file and load the definition and the initial values
#
# Arguments:
#     ask                  Ask for the name or not
#
# Returns:
#     Nothing
#
# Side effects:
#     Restores definition and initial values from external file
#
proc ::GUI::openOdeFile {ask} {
    variable filename
    variable variables
    variable documentation

    #
    # Ask for the file name
    #
    if { $ask || $filename == "" } {
        set filename [tk_getOpenFile -defaultextension ".ode" -filetypes {{.ode "ODE files"}} \
                          -title "Open the system definition" -parent .ode]
    }
    if { $filename != "" } {
        source $filename
        .ode.nb.definition.def  delete 1.0 end
        .ode.nb.definition.init delete 1.0 end
        .ode.nb.definition.def  insert end $definition
        .ode.nb.definition.init insert end $initialisation
    }
}

# exitProgram --
#     Exit the program after asking
#
# Arguments:
#     None
#
# Returns:
#     Nothing
#
# Side effects:
#     Stops the program
#
proc ::GUI::exitProgram {} {
    if { [tk_messageBox -icon question -type yesno -message "Exit the program?"] == "yes" } {
        exit
    }
}

# runCalculation --
#     Run the actual calculation
#
# Arguments:
#     show                 Show the results in the Result tab or not
#
# Returns:
#     Nothing
#
# Side effects:
#     Stores the result in the variable values
#
proc ::GUI::runCalculation {show} {
    variable values
    variable variables
    variable stateVariables

    global t
    foreach v $variables {
        global $v
    }

    set initialValues {}
    foreach v $stateVariables {
        global d_$v
        lappend initialValues [set $v]
    }

    #
    # Construct the body of the routine
    #
    set body ""
    foreach v [lrange $variables 4 end] {
        append body "global $v\n"
    }
    foreach v $stateVariables {
        append body "global d_$v\n"
    }
    append body "vexpr {[.ode.nb.definition.def get 1.0 end]}"

    proc ::GUI::Calculate {} $body

    #
    # Run the calculation
    #
    set values {}
    set step   0
    set t      $tbegin
    set tout   $tbegin
    while { $t <= $tend + 0.5 * $tstep } {
        if { abs($t-$tout) < 0.5 * $tstep } {
            set record $t
            foreach v $stateVariables {
                lappend record [set $v]
            }
            set tout [expr {$tout + $tstepout}]
            lappend values $record
        }

        Calculate

        foreach v $stateVariables {
            set current [set $v]
            set deriv   [set d_$v]
            set $v      [expr {$current + $tstep * $deriv}]
        }

        incr step
        set  t [expr {$tbegin + $step * $tstep}]
    }

    #
    # Restore the initial values
    #
    foreach v $stateVariables initv $initialValues {
        set $v $initv
    }

    if { $show } {
        .ode.nb.result.table delete 0 end
        .ode.nb.result.table deletecolumns 0 end

        set columns {0 t}
        foreach v $stateVariables {
            lappend columns 0 $v
        }
        .ode.nb.result.table configure -columns $columns

        foreach record $values {
            .ode.nb.result.table insert end $record
        }
        .ode.nb select 2
    }
}

# displayResult --
#     Display a graph of the chosen variables
#
# Arguments:
#     None
#
# Returns:
#     Nothing
#
# Side effects:
#     Draw a graph, appropriately scaled, of the results
#
proc ::GUI::displayResult {} {
    variable plot
    variable values
    variable stateVariables

    .ode.nb.plot.cnv delete all ;# Why is this needed?
    ::Plotchart::eraseplot $plot

    set xcol [.ode.nb.plot.xaxis current]
    set ycol [.ode.nb.plot.yaxis current]

    set xvalues {}
    set yvalues {}
    foreach record $values {
        lappend xvalues [lindex $record $xcol]
        lappend yvalues [lindex $record $ycol]
    }

    set xextremes [::Plotchart::determineScaleFromList $xvalues]
    set yextremes [::Plotchart::determineScaleFromList $yvalues]

    set plot [::Plotchart::createXYPlot .ode.nb.plot.cnv $xextremes $yextremes]

    foreach x $xvalues y $yvalues {
        $plot plot data $x $y
    }
}

# main --
#     Get it all started
#
catch {
    console show
}
::GUI::setupMainWindow

if { [llength $argv] > 0 } {
    set ::GUI::filename [lindex $argv 0]
    ::GUI::openOdeFile 0
}

======


<<categories>>Category Mathematics|Category Application