Arjen Markus (2 june 2004) I have been looking for ways to approach
physical problems using Tcl, but most of the interesting problems I could come with, are governed by complicated partial differential equations, such as the Navier-Stokes equations for fluid dynamics and Schroedinger's equation for quantum mechanics (the latter is not my area of expertise anyway :)).
The solution of such equations requires much more work than I can do in, say, one or two evenings.
Other simpler problems, from mechanics to thermodynamics, often involve setting up a momentum or heat balance and then solving that for a particular situation. And that is not easily captured in a static script either.
Well, that might not be entirely true, as shown below for the problem of two point masses (small balls, if you like) that collide head on. The equations that will give us the velocities
after the collision are simple:
- The total momentum must remain the same
- The total kinetic energy must remain the same
Nevertheless, we get a set of equations with two unknowns and one of the equations is quadratic.
So how do we set up these equations and solve them in the laziest way possible? Let us look at the momentum balance:
total momentum = momentum ball 1 + momentum ball 2
where:
momentum = mass * velocity
All we have to do is
minimize the function:
func = (total momentum after - total momentum before) ** 2
for the two velocities to make sure we "conserve" momentum. This is achieved at the minimum we are looking for. Then the function attains a value 0 and the total momentum before and after is the same.
What we do not want to do - at least I do not - is fill in all the formulae to make one big thing. It is so easy to make mistakes ...
And here is where Tcl's
traces come in very handy:
If you need the momentum to compute the above function,
then let the trace mechanism invoke a small procedure to
calculate that momentum.
To make this idea a bit clearer, suppose we have a block of a certain length, width and height. Then its volume is:
volume = length * width * height
If you change the width, you will need to recalculate the volume when next you need it. Now, if you put a
trace on the variable "volume", so that it gets computed when you ask for its value, like:
set width 2.0
puts $volume
you have automated the process!
That is done in the script below. Some remarks:
- You will need the code from Differentiation and steepest-descent to run it
- The solve procedure assumes that any minimum is a true solution
- It is tacitly assumed that all contributions to the function to be minimized are of the same order of magnitude - large differences may result in inaccurate solutions.
# phys.tcl --
# Solve equations without too much fuss, so that
# certain physical problems become a piece of cake
#
namespace eval ::phys {
namespace export const var solve
}
# const --
# Define a variable that does not depend on others
# Arguments:
# vname Name of the variable
# value (Optional) value
# Result:
# New variable (the one in the caller's scope is linked to the private one)
#
proc ::phys::const {vname args} {
variable $vname
uplevel [list upvar 0 ::phys::$vname $vname]
if { $args != {} } {
set $vname [lindex $args 0]
return [set $vname]
}
}
# var --
# Define a variable that _does_ depend on others
# Arguments:
# vname Name of the variable
# vexpr Expression to use when the value is required
# Result:
# New variable (the one in the caller's scope is linked to the private one)
#
proc ::phys::var {vname vexpr} {
variable $vname
uplevel [list upvar 0 ::phys::$vname $vname]
set body {
variable VNAME
set VNAME [namespace eval ::phys expr {VEXPR}]
}
proc ::phys::__set_$vname {vname1 vname2 op} \
[string map [list VNAME $vname VEXPR $vexpr] $body]
trace variable $vname r ::phys::__set_$vname
}
# solve --
# Solve a set of equations as defined by one or more variables
# Arguments:
# equations Names of the variables that define the equations
# variables Names of the variables for which to solve them
# startvalues (Optional) the initial values of these variables
# Result:
# The variables are set to their new values, approximating a solution
# the solution is printed on the screen
#
proc ::phys::solve { equations variables {startvalues {}} } {
#
# Set the initial values
#
if { $startvalues != {} } {
foreach v $variables val $startvalues {
set ::phys::$v $val
}
} else {
foreach v $variables {
lappend startvalues [set ::phys::$v]
}
}
#
# Construct a temporary procedure that defines the
# function we will minimize:
#
# 1. The variable commands
#
set ns_vars ""
foreach v $equations {
append ns_vars "variable $v\n"
}
#
# 2. The function: the sum of squares of each equation
#
set eqs 0.0
foreach e $equations {
append eqs "+\$$e*\$$e"
}
#
# 3. The procedure
#
proc ::phys::__minfunc args "$ns_vars; \
foreach v [list $variables] val \$args {set ::phys::\$v \$val} ; \
expr {$eqs}"
#
# Now, minimize the function and transfer the result to the variables
#
set result [::math::optimize::minimumSteepestDescent ::phys::__minfunc $startvalues]
foreach v $variables val $result {
set ::phys::$v $val
puts "$v = $val"
}
}
#
# A simple problem: the volume of a block
#
namespace import ::phys::*
source differentiate.tcl
const width 10.0
const length 3.0
const height 7.0
var volume {$width*$height*$length}
puts "Volume: $volume"
set width 4.0
puts "Volume: $volume"
#
# Now something less trivial:
# Calculate the width of a block with given volume, height and length
#
var dvol {$volume-100.0}
solve dvol width 1
#
# Just for the fun of it:
# Compute the velocities after a frontal collision of
# two point bodies with equal mass (to make the equations simpler)
# The trick:
# - the momentum before and after the collision must be the same
# - the kinetic energy must be the same
#
const start_vel1 1.0
const start_vel2 0.1
const vel1
const vel2
var start_momentum {$start_vel1+$start_vel2}
var start_kin_energy {$start_vel1*$start_vel1+$start_vel2*$start_vel2}
var end_momentum {$vel1+$vel2}
var end_kin_energy {$vel1*$vel1+$vel2*$vel2}
var dmom {$end_momentum-$start_momentum}
var dkin {$end_kin_energy-$start_kin_energy}
solve {dmom dkin} {vel1 vel2} {0.0 10.0}
#
# What remains to be done: check that we have indeed
# found an acceptable solution!
#
