[GWM] PDE's express a function of at least 2 variables so that the value of the function has some sort of relationship to values at nearby points. 3 common methods of solution are Finite Element, Finite Volume & Finite Difference methods. Try a google search for these names. This page (will) shows how a simple PDE can be solved numerically. The next commonest method is Boundary Element methods which are reputed to be fast but need a good understanding of why they work at all. [GWM] 8 Sep 2004 minor change to code to remove dead variable and allow user to set relaxation coefficient. Added new reference. FD - http://mathworld.wolfram.com/FiniteDifference.html FE - despite being the most widely used PDE solver, very few good references http://www.fact-index.com/f/fi/finite_element_method.html FV - http://mathworld.wolfram.com/FiniteVolumeMethod.html BEM - http://www.boundary-element-method.com/ ---- The heat conduction equation is: Tt = Txx +Tyy (where Txx implies differentiation with respect to x twice, Tt diferentiate temperature with time.). Very few of the best solution methods are suitable for Tcl - they could be implemented but the sheer volume of data and storage plus solution time makes a long wait. {An old English custom for new starts in an engineering company was to send the new employee to the stores for a long weight; they would ask (say it aloud)for a long wait, and after a few hours the new start might notice that nothing had happened - clearly management material.} One of the few solution methods that can be implemented easily is a relaxation method for the Heat equation on a 2D rectangular grid. We can approximate (FD!) the above equations using [Euler Methods] forward marching in time with: Tnew(x,y) = Told(x,y) + k.(T(x+1,y)-T(x,y)-T(x,y)+T(x-1,y) +T(x,y+1)-T(x,y)-T(x,y)+T(x,y-1)) = Told(x,y) + k.(T(x+1,y)-4 T(x,y)+ T(x-1,y) +T(x,y+1) + T(x,y-1)) Why is that equation right? - dT/dt implies (Euler) that t(new)=told +dt*(RHS of equation). Txrt (differential of T to right of point) = (T(x+dx)-T(x))/dx Txlt (differential of T to left of point) = (T(x)-T(x-dx))/dx, SO: Txx = (Txrt - Txlt)/dx = T(x,y+1)-2*T(x,y)+T(x-1,y). Add a similar term for Tyy and there you are. k is a constant relating speed of transmission of heat to the time step of the solution (rate of transmission of heat per unit of space dimension). In the code below k is 0.25 (argument kdt to proc nexttime) - if you increase k to >0.25 (try 0.3) the equations become numerically unstable, and after a few steps the solver will die as one value will exceed the largest storage (you could amend this solver sot hat kdt is changed after step 100 or so). So remember: you can write a discrete version of an equation but it might not work! There are long theoretical discussions about why some values work and others dont, and the reasons are good but a little complex for this space. (See Crank-Nicholson, Implicit methods and so on for more numerically stable results.) A good description of PDE solvers for this equation is at http://www.upscale.utoronto.ca/GeneralInterest/Harrison/HeatEqn.pdf Here is a simple algorithm to demonstrate how a hot point dissipates. I set a block of pixels in an image to a very high intensity, then apply the FD equation to observe the spread of the 'heat' = seen as a spreading and cooling area of color. If you are using a machine slower than about 2GHz, try reducing dim in 'proc heateq' from 100 to 50. And the number of time steps from 350 to a manageable number. This is not going to proceed like lightning I'm afraid! But 10 minutes would have been fast 20 years ago on a mainframe, and this takes only 4 on a home PC. proc setcolour { tfd data } { ;# map 0-10000 to red->blue # push new pixel colours into the rectangular solution domain (or photo) upvar 1 $tfd f ;# field (image) values in the calling routine upvar 1 $data t ;# calculated values in the calling routine set nh [llength $t] ;# No of rows set nx [llength [lindex $t 0]] ;# No of cols for {set i 0} {$i<$nh} { incr i} { set pix [list] set row [lindex $t $i] for {set j 0} {$j<$nx} { incr j} { set tem [lindex $row $j] if {[expr $tem > 300000]} then { lappend pix \#ff0000 } \ elseif {[expr $tem > 40000]} then {lappend pix \#dd6600} \ elseif {$tem > 15000} then {lappend pix \#cccc00} \ elseif {$tem > 6000} then {lappend pix \#88dd00} \ elseif {$tem > 3200} then {lappend pix \#00ff00} \ elseif {$tem > 1000} then {lappend pix \#00bb33} \ elseif {$tem > 300} then {lappend pix \#009966} \ elseif {$tem > 50} then {lappend pix \#0033aa} \ elseif {$tem > 10} then {lappend pix \#0000ff} \ else {lappend pix \#aabbcc} } $tfd put [list $pix] -to 1 $i } } proc nexttime {data kdt} { ;# advect heat # data is the array of temperatures on a rectangular grid, stored as a list of rows, each row a list of temperatures # kdt is heat conductivity.dtime - value>0.25 is numerically unstable upvar 1 $data t ;# values in the calling routine set nh [llength $t] ;# No of rows set tnew {} lappend tnew [lindex $t 0] ;# copy old first line - boundary condition set nx [llength [lindex $t 0] ] ;# No of cols for {set i 1} {$i<$nh-1} { incr i} { set pix {} set row [lindex $t $i] set rowa [lindex $t [expr $i-1] ] set rowb [lindex $t [expr $i+1] ] lappend pix [lindex $row 0] for {set j 1} {$j<$nx-1} { incr j} { set txy [lindex $row $j] set txym1 [lindex $rowa $j] set txyp1 [lindex $rowb $j] set txm1y [lindex $row [expr $j-1]] set txp1y [lindex $row [expr $j+1]] set tn [expr int($txy+$kdt*($txym1 + $txyp1 + $txm1y + $txp1y -4*$txy))] lappend pix $tn } lappend pix [lindex $row [expr $nx-1]] lappend tnew $pix } lappend tnew [lindex $t [expr $nh-1]] return $tnew } proc heateq {} { catch [destroy .temp] {} ;# delete it. Same as clear. catch [destroy .lob] {} ;# delete it. Same as clear. set dim 100 ;# make smaller for smaller domain (and slower computer) set mid [expr $dim/2] set hot {} set row {} for {set j 0} { $j<$dim } {incr j} { lappend row 0 lappend hot [expr ($j>$mid-2 && $j<$mid+2) ? 999999 : 0] } for {set i 0} { $i<$dim } {incr i} { if {$i >$mid-2&& $i<$mid+2} { lappend tdata $hot } else { lappend tdata $row } } image create photo tempfield -height $dim -width $dim setcolour tempfield tdata label .lob -text Temperature label .temp -image tempfield -bd 3 -relief raised pack .lob .temp -side top update for {set t 0 } {$t<350} {incr t} { set tdata [nexttime tdata .25] setcolour tempfield tdata update puts "Time $t " } } time { heateq } Close inspection of this code shows a similarity to cellular automata such as [Life], at each time step the 'colour' of each cell is updated dependent on neighbour cells. ---- [AM] Could we join forces? I would love a PDE solver in Tcl - even it is not that particularly fast or even remotely useful :) As you have not signed this page with your name nor any of the others on similar subject, I have to take this public route. [GWM] Sorry - forgot to sign it, just like my letters to the taxman. I think a general PDE solver would be very hard to implement, and yes it is less than useful but educational is still useful. [AM] Well, I know of no completely general PDE solver, most concentrate on particular classical equations and variations thereof. But: a generic solver for typical classical PDEs should be doable. [GM] OK. I am currently working for a company called CHAM who sell a general CFD & stress equation solver, but if you look inside the hood there are a lot of different solution methods which are semi-automatically chosen by the software. That code uses about 1 million lines of FORTRAN and C. Certainly a range of 'matrix' like routines to handle the nearest neighbour values and produce the correction terms for d/dx, d2/dx2 under different basic solvers would be a start. Also a general 'contour' plotting program (as in the code above but better) could be useful.