Version 5 of Partial Differential Equations

Updated 2004-10-07 21:44:59

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.

 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) + a.(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) + a.(T(x+1,y)-4 T(x,y)+ T(x-1,y) +T(x,y+1) + T(x,y-1))

a 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 a is 0.25 (see line beginning 'set tn') - if you increase a 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.)

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, and this takes only 2.

 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} { ;# advect heat
        upvar  $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+.25*($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 pxdata {}
        set row {}
        for {set j 0} { $j<$dim } {incr j} {
                lappend row 0
                lappend pxdata #004488
                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]
                setcolour tempfield tdata
                update
                puts "Time $t "
        }
 }

 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.