[Arjen Markus] (20 february 2007) I would like to claim that the package presented below is a general package for
solving partial differential equations, but that would be a gross exaggeration:

   * It is not ready for any practical use yet
   * When/if it is finished it will still only deal with a particular class of PDEs
   * It will use finite-volume methods on rectangular grids, so the PDE must be fit for such a treatment
   * It will be limited to one and two spatial dimensions for practical reasons

But on the other hand:

   * You will be able to deal with systems of PDEs
   * Both Dirichlet and Neumann boundary conditions are supported
   * It will easy to define the particular PDEs 

Well, the type of PDE it can solve can be expressed like this:

      du                d2u              d2u               du              du
     ----  = f(x,y,u) * --- + g(x,y,u) * --- + h(x,y,u) * ---- + j(x,y,u) ---- + k(x,y,u)
      dt                dx2              dy2               dx              dy

where u is the dependent variable or a set of such variables.

Without further ado, here is the code I have so far (limited to 1D, no first-order derivatives yet). 
It relies on [NAP] for the heavy-duty computations. 

----

 # pde.tcl --
 #     Implementation of a basic solver for partial differential
 #     equations (PDE), using NAP
 #
 #     This solver is meant for quasilinear parabolic PDEs, but
 #     it can handle systems of PDes. Via a simple transformation
 #     a hyperbolic equation, like the wave equation, can be
 #     written as a system of two or more equations that are first-order
 #     in the time.
 #
 #     The current implementation is limited to equidistant grids in
 #     one spatial dimension. The types of boundary conditions are
 #     - Dirichlet (value given at the boundary)
 #     - Neumann (value of the flux given at the boundary)
 #     The same type is assumed for all equations, but the boundaries
 #     are independent.
 #
 #     Because of the amount of information involved in solving PDEs
 #     we store everything in a separate namespace (this is partly
 #     required because of the properties of NAP). This is one of the
 #     tasks taken care of in the initialisation phase.
 #
 #     TODO:
 #     nextTime
 #     getValues
 #
 package require nap

 namespace eval PDE {
     namespace import ::NAP::*

     variable pdeCount 0

     #
     # Public procedures
     #
     namespace export initPde1D boundary timestep values initial
 }

 # initPde1D --
 #     Create a namespace to solve a one-dimensional PDE
 # Arguments:
 #     None
 # Result:
 #     The name of the namespace dedicated to the PDE
 # Note:
 #     The result is required for accessing the PDE solution system.
 #     Pass it as the first argument.
 #
 proc ::PDE::initPde1D {} {
     variable pdeCount

     namespace eval $pdeCount [list     \
         variable _dimension 1          \
         variable _pde ::PDE::$pdeCount \
     ]

     set ns "::PDE::$pdeCount"
     incr pdeCount

     return $ns
 }

 # timestep --
 #     Set the value of the time step for integration
 # Arguments:
 #     pde              Name of the PDE solution system
 #     deltt            Time step to be used
 # Result:
 #     None
 #
 proc ::PDE::timestep {pde deltt} {

     set [set pde]::_deltt $deltt
 }

 # printPde --
 #     Print the PDE information array
 # Arguments:
 #     pde              Name of the array holding all information
 # Result:
 #     None
 # Side effects:
 #     Prints the information contained in the array on stdout
 #
 proc ::PDE::printPde {pde} {

     set [set pde]::_pde $pde
     namespace eval $pde {
         puts "PDE: $_pde"
         puts "    -- grid coordinates (cell boundaries) --"
         puts "[$_grid value]"
         puts "    -- volumes of cells --"
         puts "[$_volume value]"
         puts "    -- grid distances --"
         puts "[$_dx value]"

         if { [info exists _statevars] } {
             puts "    -- state variables --"
             foreach v $_statevars {
                 puts "$v: [[set $v] value]"
             }
         }
     }
 }

 # values --
 #     Return the values of a state variables as a list
 # Arguments:
 #     pde              Name of the array holding all information
 #     name             Name of the state variable
 # Result:
 #     None
 #
 proc ::PDE::values {pde name} {

     ::PDE::HasParameter $_pde $name

     set [set pde]::_name $name
     namespace eval $pde {
         return [[set $_name] value]
     }
 }

 # HasParameter --
 #     Check that a correct parameter name is given
 # Arguments:
 #     pde          Name of the PDE system
 #     name         Parameter name (empty: any parameter)
 #     hasvalue     Does it have a value? (optional)
 # Result:
 #     1 if there is a parameter by that name, error otherwise
 #
 proc ::PDE::HasParameter {pde name {hasvalue 0}} {

     set [set pde]::_name     $name
     set [set pde]::_hasvalue $hasvalue

     namespace eval $pde {
         if { ! [info exists _statevars] } {
             return -code error "No state variables defined yet"
         }
         if { $_name != {} } {
             if { [lsearch $_statevars $_name] < 0 } {
                 return -code error "No such state variable: $_name"
             }
             if { $_hasvalue && ! [info exists $_name] } {
                 return -code error "State variable $_name has no values yet"
             }
         }
     }
 }

 # HasGrid --
 #     Check that the PDE system has a grid
 # Arguments:
 #     pde          Name of the PDE system
 # Result:
 #     1 if there is a grid, error otherwise
 #
 proc ::PDE::HasGrid {pde} {
     namespace eval $pde {
         if { ! [info exists _grid] } {
             return -code error "No grid defined yet"
         }
     }
 }

 # boundary --
 #     Set the type and value for the boundaries
 # Arguments:
 #     pde              Name of the array holding all information
 #     name             Which parameter
 #     bound            Which boundary (left, right, top or bottom)
 #     type             Type of condition (dirichlet or neumann)
 #     value            Value of the variable or the derivative
 # Result:
 #     None
 #
 proc ::PDE::boundary {pde name bound type value} {
     switch -re $bound {
         "l|left"   { set b left   }
         "r|right"  { set b right  }
         "t|top"    { set b top    }
         "b|bottom" { set b bottom }
         default {
             return -code error "Boundary must be left, right, top or bottom"
         }
     }

     switch -re $type {
         "d|dirichlet" { set t dirichlet }
         "n|neumann"   { set n neumann   }
         default {
             return -code error "Type must be dirichlet or neumann"
         }
     }

     set [set pde]::_btype($b,$name)  $t
     set [set pde]::_bvalue($b,$name) $value
 }

 # SetBoundaryValues --
 #     Set the correct values for the boundary elements
 # Arguments:
 #     pde          Name of the PDE system
 # Result:
 #     None
 #
 proc ::PDE::SetBoundaryValues {pde} {

     #
     # TODO: Neumann boundaries!
     #

     namespace eval $pde {
         set _maxcellp1 [expr {$_maxcell + 1}]
         foreach _v $_statevars {
             foreach _b {left right} _idx [list 0 $_maxcellp1] {
                 if { $_btype($_b,$_v) == "dirichlet" } {
                     [set $_v] set value $_bvalue($_b,$_v) $_idx
                 }
             }
         }
     }
 }

 # initial --
 #     Set the initial value for a particular state variable
 # Arguments:
 #     pde              Name of the array holding all information
 #     name             Which variable
 #     value            Initial value
 # Result:
 #     None
 # Notes:
 #     Useful only when the grid has been defined and the state variables
 #
 proc ::PDE::initial {pde name value} {

     set [set pde]::_name  $name
     set [set pde]::_value $value

     HasParameter $pde $name
     HasGrid $pde

     namespace eval $pde {
         set $_name $_value
         nap "$_name = $_value + 0.0 * $_volume"
     }
 }

 # zeroDerivs --
 #     Set the derivatives for all state variables to zero
 # Arguments:
 #     pde              Name of the array holding all information
 # Result:
 #     None
 #
 proc ::PDE::zeroDerivs {pde} {

     HasParameter $pde {}
     HasGrid $pde

     namespace eval $pde {
         foreach name $_statevars {
             nap "_d_$name = 0.0 * $_volume"
         }
     }
 }

 # addTerm --
 #     Add a term to the derivative of a particular state variable
 # Arguments:
 #     pde              Name of the array holding all information
 #     name             Which variable
 #     term             The value in each cell to be added
 # Result:
 #     None
 #
 proc ::PDE::addTerm {pde name term} {

     HasParameter $pde $name 1
     HasGrid $pde

     set [set pde]::_name $name
     set [set pde]::_term $term

     namespace eval $pde {
         [set _d_$_name] \
             set value [string map [list Name $_name] "_d_Name(1 .. _maxcell) + _term"] \
                 "(1 .. _maxcell)"
     }
 }

 # addSecondDeriv --
 #     Compute the second order derivative of a particular state variable,
 #     multiplied by a factor (for convenience)
 # Arguments:
 #     pde              Name of the array holding all information
 #     name             Which variable
 #     factor           Expression for multiplication factor
 #     dname            Which variable to take the derivative of
 # Result:
 #     None
 #
 proc ::PDE::addSecondDeriv {pde name factor dname} {

     set [set pde]::_name   $name
     set [set pde]::_dname  $dname
     set [set pde]::_factor $factor

     HasParameter $pde $name 1
     HasParameter $pde $dname 1
     HasGrid $pde

     SetBoundaryValues $pde

     namespace eval $pde {
         set _maxcellm1 [expr {$_maxcell-1}]
         nap [string map [list Name $_dname] \
                "_term1 = (Name(1 .. _maxcell) - Name(0 .. _maxcellm1)) /
                          _dx(1 .. _maxcell)"]
         nap [string map [list Name $_dname] \
                "_term2 = (Name(2 .. _nocells) - Name(1 .. _maxcell)) /
                          _dx(2 .. _nocells)"]
         nap "_term  = $_factor * (_term2 - _term1)"

         [set _d_$_name] \
             set value [string map [list Name $_name] "_d_Name(1 .. _maxcell) + _term"] \
                 "(1 .. _maxcell)"
     }
 }

 # nextStep --
 #     Compute for all state variables the values in the next time step
 # Arguments:
 #     pde              Name of the array holding all information
 # Result:
 #     None
 #
 proc ::PDE::nextStep {pde} {

     HasParameter $pde {} 0
     HasGrid $pde

     namespace eval $pde {
         foreach _v $_statevars {
             nap "$_v = $_v + $_deltt * _d_$_v"
         }
     }
 }

 # statevars --
 #     Set the names of the state variables
 # Arguments:
 #     pde              Name of the array holding all information
 #     args             Names of the variables
 # Result:
 #     None
 #
 proc ::PDE::statevars {pde args} {

     foreach name $args {
         if { [string index $name 0] == "_" } {
             return -code error "Names starting with an underscore are reserved - do not use them for state variables"
         }
     }

     set [set pde]::_statevars $args
 }

 # grid1D --
 #     Create a one-dimensional grid
 # Arguments:
 #     pde              Name of the array holding all information
 #     nocells          Number of grid cells
 #     xmin             Minimum x-coordinate (left boundary)
 #     xmax             Maximum x-coordinate (right boundary)
 # Result:
 #     None
 #
 proc ::PDE::grid1D {pde nocells xmin xmax} {

     #
     # Create a namespace to hold the NAP variables
     #
     namespace eval $pde {
         namespace import ::NAP::nap
     }

     #
     # Check the consistency
     #

     if { $nocells < 1 } {
         return -code error "Number of cells must be at least 1"
     }
     if { $xmin >= $xmax } {
         return -code error "Minimum x-coordinate must be smaller than maximum"
     }

     set [set pde]::_nocells $nocells
     set [set pde]::_maxcell [expr {$nocells-1}]
     set [set pde]::_xmax    $xmax
     set [set pde]::_xmin    $xmin

     namespace eval $pde {
         set _deltx [expr {($_xmax-$_xmin)/double($_nocells)}]
         nap "_grid   = _xmin + _deltx * (0 .. _nocells)"
         nap "_dx     = _deltx + 0.0 * (0 .. _maxcell)"
         nap "_volume = _deltx + 0.0 * (0 .. _maxcell)"
     }
 }

 # main --
 #     Testing:
 #     - Generation of the grid must work properly
 #
 set g [::PDE::initPde1D]
 ::PDE::grid1D $g 20 -10 10

 ::PDE::statevars $g u
 ::PDE::initial   $g u 2.1
 ::PDE::boundary  $g u left dirichlet 0
 ::PDE::boundary  $g u right dirichlet 0

 ::PDE::timestep  $g 0.1

 ::PDE::printPde $g

 set i 0
 while { $i < 10 } {
     ::PDE::zeroDerivs $g
     ::PDE::addTerm $g u 1.0
     ::PDE::addSecondDeriv $g u 1.0 u
     ::PDE::nextStep $g
     incr i
 }

 ::PDE::printPde $g

----
[[
[Category Mathematics]
]]