PDEs and the Tensor package

Arjen Markus (22 october 2008) Recently Neil McKay announced version 4.0a1 of his Tensor package on c.l.t. I thought it might be a good idea to see if this package that deals with vectors, matrices and their higher-dimensional cousins, could be used for solving partial differential equations (PDEs) and what kind of code that would produce.

The program below is just a quick-and-dirty attempt to get some feeling for the package. Some remarks:

• Tensors that you create are commands and therefore persistent. This means that you will have to clean them up explicitly
• The assignments do not create new tensors, they merely set the values contained in the tensors.
• The tensor::expr command does not like to deal with Tcl variables directly, so I used double quotes (") instead of the more usual {} around my expressions.
• The tensor package seems to lack some functionality still, notably a convenient way to get the data

But apart from that, the experiment, as far as I am concerned, was a success. This package does hold some promises in this area ...

(Note: The package is available for Tcl 8.5 via Teapot and the source code is at http://www.eecs.umich.edu/~mckay/computer/tensor4.0a1.tar.gz )

# tensors.tcl --
#     Use the Tensor package to solve a simple diffusion problem
#
package require Tensor

# initialise --
#     Create the tensors
#
# Arguments:
#     value           Value at cell x=0, t=0
#     noCells         Number of cells
#
# Result:
#     None
#
# Side effects:
#     Creates the tensors solution, work and deriv
#
proc initialise {value noCells} {

    ::tensor::create solution -type double -size [expr {$noCells+2}]
    ::tensor::create deriv    -type double -size [expr {$noCells+2}]
    ::tensor::create work     -type double -size [expr {$noCells+2}]

    set zero [expr {$noCells/2}]
    solution = scalar 0.0
    solution section $zero = scalar $value

    work  = scalar 0.0
    deriv = scalar 0.0
}

# nextTime --
#     Compute the solution for the next time
#
# Arguments:
#     diff            Diffusion coefficient (m2/s)
#     dx              Size of a grid cell
#     dt              Time step (s)
#
# Result:
#     None
#
# Side effects:
#     An updated solution
#
proc nextTime {diff dx dt} {

    set maxCell   [expr {[solution dimensions] - 1}]
    set maxCellM1 [expr {$maxCell - 1}]
    set maxCellM2 [expr {$maxCell - 2}]

    ::tensor::expr "work(i=1:$maxCellM1) = -2.0 * solution(i=1:$maxCellM1)
                                                  + solution(i=0:$maxCellM2)
                                                  + solution(i=2:$maxCell)"
    deriv = scalar 0.0
    ::tensor::expr "deriv(i) = $diff * work(i) / $dx / $dx"
    ::tensor::expr "solution(i) = solution(i) + $dt * deriv(i)"
}

# main --
#     Run it
#
initialise 1.0 10

set dx   0.1
set dt   0.01
set diff 0.1

for { set i 0 } { $i < 10 } { incr i } {
   nextTime $diff $dx $dt
}

#solution writechannel stdout text ;# Did not print anything useful ...

solution put into a
parray a 

NDM (07 January 2009) Just a couple of comments. First, the line

    ::tensor::expr "solution(i) = solution(i) + $dt * deriv(i)"

could be written as

    ::tensor::expr "solution(i) += $dt * deriv(i)"

Second, to get a tensor value as a nested Tcl list, just use the tensor's name as a command, for example

    puts [solution]

TV For a ´package´ it would be nice to define what you mean by 'tensor' and what types of Differential Equations you mean to solve with them. The idea of the combination of a space spanned by linear interpolation and matrix or algebraic manipulations is at least a century old, and I don´t think you´re suddenly going to solve and important instance of this kind of work (if you wish to use normal definitions from that time), namely Einsteins general relativity. The number of possible DEs reminds me of a phone-book thick textbook from my second year (EE) and I don´t think you´ll cover that either, at least not without an unchristian amount of work.

As an 'example' I´m sure it´s fun to take up an ancient idea and put it in a decent computer language for some purpose.

AM (3 april 2009) Tensors are a generalisation of vectors and matrices, with appropriate definitions for adding and multiplying and so on. Nothing really mysterious there and no inherent relation to Einstein ;). As for the demarcation of PDEs: given that some of my colleagues and many other people in the world make a living out of developing efficient and effective methods to solve them, I do not pretend to come up with a Package To Replace All Others. It would just be nice to deal with the simpler ones using Tcl.