[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: ======none total momentum = momentum ball 1 + momentum ball 2 ====== where: ======none momentum = mass * velocity ====== All we have to do is ''minimize'' the function: ======none 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: ======none 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: ======none 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! # ====== <> Mathematics | Numerical Analysis | Physics | Science