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. 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 }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.