## Bernoulli using math::calculus

GeoffM My derivation from Bernoulli.

I modified the basic plotf to handle multiple graphs, add a title, set colours and other options. Then added a call to RungeKutta4 to solve the damped harmonic oscillator equation (as an example).

```  package require Tk
package require Plotchart
console show

# plotf --
#     Convenience procedure to plot a function of one variable or
# call a user defined proc with arguments:
#
# Arguments:
#     range       Range over which to plot the function (x-axis)
#     expression  Expression by which the function is defined
#                 (independent variable: t)
#
proc plotf {range expression args} {
# intial is the initial conditions
# nosteps - number of time steps to progress.
# args are arguments to be interpreted; currently
#  -title "XXX"
# -nosteps Number of steps to use

# Very simple for now: no parameters
set nosteps 100
foreach {name value} \$args {
switch -- \$name {
"-nosteps" {set nosteps \$value}
"-initial" {set x \$value}
default {}
}
}

foreach {tmin tmax} \$range {break}
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}]
if {[info command \$expression]!={}} {
# Note that RungeKutta expects start of time step as time
set x [\$expression [expr {\$t-\$tstep}] \$x \$tstep]
} else {
set x [expr \$expression]
}
lappend tvalues \$t
lappend xvalues \$x
if {\$xmin=={}} {set nvalues [llength \$x]}
foreach xvalu \$x { ;# handle vector of results
if { \$xmin == {} || \$xvalu < \$xmin } {
set xmin \$xvalu
}
if { \$xmax == {} || \$xvalu > \$xmax } {
set xmax \$xvalu
}
}
}

set tscale [::Plotchart::determineScale \$tmin \$tmax]
set xscale [::Plotchart::determineScale \$xmin \$xmax]

set p [::Plotchart::createXYPlot .c \$tscale \$xscale]
set i 0
foreach {name value} \$args {
switch -- \$name {
"-title" { .c create text 200 50 -text \$value }
"-colours" -
"-colors" {
set i 0
foreach col \$value {
\$p dataconfig function\$i -colour \$col
incr i
}
}
default {}
}
}
foreach t \$tvalues x \$xvalues {
set i 0
while {\$i<\$nvalues} { ;# handle vector of results
\$p plot function\$i \$t [lindex \$x \$i]
incr i
}
}
}

# now define a runge-kutta- solver
package require math::calculus  ;# supplies RK4 algorithm

proc dampened_oscillator {t xvec} {
# called by RungeKuttaStep. Returns list of derivatives of xvec.
global damp
global k
set x  [lindex \$xvec 0] ;# pos
set speed [lindex \$xvec 1] ;# speed
# Solves the damped oscillator equation: x'' = -k*x-damp*speed
# where x is position, and x'' the second time derivative.
#  This equation is broken into 2 first order equations:
# dx/dt=v
# dv/dt= k.x-damp*speed
#   k is the 'spring constant' of the oscillator;
# more negative gives higher oscillation frequency.
# a list of derivatives at time t.
return [list \$speed [expr {\$k*\$x-\$damp*\$speed}]]
}

proc osc {t x tstep} { ;# return the value at the time t+dt for step tstep
# x is a list of values (position, speed for the moving part)
return [::math::calculus::rungeKuttaStep \$t \$tstep \$x dampened_oscillator]
}
set damp -0.2
set k -9