FD - http://mathworld.wolfram.com/FiniteDifference.htmlFE - despite being the most widely used PDE solver, very few good references http://www.fact-index.com/f/fi/finite_element_method.html
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]} { lappend pix \#ff0000
} elseif {[expr $tem > 40000]} {lappend pix \#dd6600
} elseif {$tem > 15000} {lappend pix \#cccc00
} elseif {$tem > 6000} {lappend pix \#88dd00
} elseif {$tem > 3200} {lappend pix \#00ff00
} elseif {$tem > 1000} {lappend pix \#00bb33
} elseif {$tem > 300} {lappend pix \#009966
} elseif {$tem > 50} {lappend pix \#0033aa
} elseif {$tem > 10} {lappend pix \#0000ff
} else {lappend pix \#aabbcc}
}
$tfd put [list $pix] -to 0 $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
} elseif {$i >$dim-15 && $i<$dim-12} { ;# add extra hot at edge of domain
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.GWM 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.
AM If we get serious about this, and why shouldn't we? :), then it is very useful to have a benchmark by which we can decide what is the fastest way to update the blocks of data that are typical for PDEs. I have set up the page Partial Differential Equations - performance benchmarks to record our experience.
AM (10 february 2005) I finally sat down and experimented with the NAP extension by Harvey Davies that is designed to work with multidimensional arrays. The code below is very immature as using NAP requires a slightly different mindset than plain Tcl. (Also, from a numerical point of view, I would rather use a finite-volumes approach than the finite-differences method below.) Still, I think it illustrates the power of NAP: deal with large amounts of data in a functional way instead of a procedural one.
# Experiment with NAP:
# A simple PDE solver for the diffusion equation in
# two dimensions
#
package require nap
proc make_grid {{n 20}} {
nap "reshape(0.0,2#n)"
}
proc initial_condition {conc value} {
foreach {n1 n2} [$conc shape] {break}
$conc set value $value [expr {$n1/2}],[expr {$n2/2}]
}
proc laplacian {cnc} {
upvar #0 laplace laplace
upvar #0 $cnc conc
foreach {n1 n2} [$cnc shape] {break}
nap "col = 0.0*(0..n2)"
nap "row = 0.0*(0..n1)"
nap "idxc = 1..(n1-2)"
nap "idxl = 0..(n1-3)"
nap "idxr = 2..(n1-1)"
nap "d2cdx2 = cnc(idxl,)+cnc(idxr,)-2*cnc(idxc,)"
nap "d2cdx2 = (col // d2cdx2) // col"
nap "idxc = 1..(n2-2)"
nap "idxd = 0..(n2-3)"
nap "idxu = 2..(n2-1)"
nap "d2cdy2 = cnc(,idxd)+cnc(,idxu)-2*cnc(,idxc)"
nap "d2cdy2 = transpose( (row // transpose(d2cdy2)) // row )"
nap "laplace = d2cdx2 + d2cdy2"
return $laplace
}
nap "conc = make_grid(20)"
initial_condition $conc 100.0
laplacian $conc
nap "conc_new = conc + 0.025*laplace"
puts [$conc_new value]TV I distinctly remember my 2d years EE course in differential equations, and the book about as thick as a phonebook with various kinds of them: I'm sure it is very hard to aim for some generality when staying with analytical methods, not even to speak of dealing with fundamental electrical circuits and their corresponding DE and solutions, I might life up to my promise at some point to make a tcl implementation about that problem.For iterative solutions, it might be good to include analytical symbolic) analysis, which is possible with maxima, which can be combined with tcl (e.g. Calling Maxima from Tcl).An maxima example:
ATVALUE(fe(t),t=0,1);
ATVALUE('diff(fe(t),t),t=0,1/2);
desolve([fe(t)=-'diff(fe(t),t,2)],[fe(t)]);
SIN(t)
--> fe(t) = ------ + COS(t)
2(see [1])EKB As near as I can tell from the Maxima documentation, Maxima only does ODEs, not PDEs. But the general idea of using an external library -- would that be useful? E.g., search for "PDE" at GAMS, the Guide to Available Mathematical Software. And there are probably some already existing libraries for finite element methods.TV The interesting point with maxima (like mathematica, I haven't linked that with tcl though) is that it can work symbolic, and with infinite accuracy if you like, which is stronger than a numerical solver.Lars H, 21 July 2005: Since symbolic metods only produce solutions to very tiny classes of differential equations (already for ODEs, and very definitely in the case of PDEs), that comparison is not a particularly enlightened one. Yes you may get an exact symbolic expression when the problem is easy, but most of the time you get nothing, which of course is far less than a numerical approximation. (Yet I'm going to spend seven weeks this autumn teaching symbolic ODE and PDE solution methods. Sigh. Well, it mostly boils down to linear algebra, and that's generally useful.)Even for those PDEs where one can produce symbolic "solutions", there is really a lot of approximation that has been swept under the carpet: eigenvalue equations that cannot be solved symbolically, Fourier coefficients that can only be computed via numeric integration, and Fourier series that have to be truncated just to name a few. The combined errors from these may well be larger than the approximation errors in a purely numeric solution, but since they're spread out all over they're much harder to analyse.TV My point was more in the analysis and the general idea of being able to construct solutions instead of leavingit all to some big named iteration scheme to solve a problem without insight. That most certainly holds for linear algebra based solutions, too. A major branch of science is supported by the theory of an symbolically computable set of DEs: electrical network theory.AM (20 february 2007) I started the development of a more-or-less general package for solving PDEs It uses NAP to do the heavey-duty computations.
