sci.math.numanalysis
SUNA59, Numerical Differentiation
=================================
In article 56 of the Series on Unified Numerical Approximations I wrote:
> A piece of postprocessing is done also, in order to complete the AVS UCD file
> (which is meant for a course I'm preparing) with "cell data". Three equations
> are set up, according to the theory of "Stream Function & Flow Potential"
> (14), when extended to three dimensions. All three are of the form:
>
> P2  P1 = (x2  x1).u + (y2  y1).v + (z2  z1).w
>
> In order to have a reasonable condition of the 3 equations, the midpoints of
> the brick's faces are used for evaluating the potential function values.
The twodimensional analogue of the above is depicted below:
3 4
o  *  o [1/2(x2+x4)1/2(x1+x3)].u + [1/2(y2+y4)1/2(y1+y3)].v =
   1/2(P2+P4)1/2(P1+P3)
*  u,v  *
   [1/2(x3+x4)1/2(x1+x2)].u + [1/2(y3+y4)+1/2(y1+y2)].v =
o  *  o 1/2(P3+P4)1/2(P1+P2)
1 2
> \^/ MAPLE V
> ._\ /_. Copyright (c) 19811990 by the University of Waterloo.
> \ MAPLE / All rights reserved. MAPLE is a registered trademark of
> <____ ____> Waterloo Maple Software.
>
>> with(linalg);
>> A := array([[ (x2+x4x1x3)/2 , (y2+y4y1y3)/2 ],
>> [ (x3+x4x1x2)/2 , (y3+y4y1y2)/2 ]]);
>> det(A);
>
> 1/2 x2 y4  1/2 x2 y1 + 1/2 x4 y3  1/2 x4 y2  1/2 x1 y3 + 1/2 x1 y2
>
>  1/2 x3 y4 + 1/2 x3 y1
The determinant is sum of 8 products. It is used when employing Cramer's rule:
substitute P for (x,y) , giving the nominator for (u,v) , respectively.
DET = 1/2 (x2.y4  x2.y1 + x4.y3  x4.y2  x1.y3 + x1.y2  x3.y4 + x3.y1)
u = 1/2 (P2.y4  P2.y1 + P4.y3  P4.y2  P1.y3 + P1.y2  P3.y4 + P3.y1) / DET
v = 1/2 (x2.P4  x2.P1 + x4.P3  x4.P2  x1.P3 + x1.P2  x3.P4 + x3.P1) / DET
We repeated this exercise with MAPLE for the threedimensional case, which was
significant for our theory in article (56).
>> A := array([[ (x5+x6+x7+x8x1x2x3x4)/4 ,
>> (y5+y6+y7+y8y1y2y3y4)/4 ,
>> (z5+z6+z7+z8z1z2z3z4)/4 ],
>> [ (x3+x4+x7+x8x1x2x5x6)/4 ,
>> (y3+y4+y7+y8y1y2y5y6)/4 ,
>> (z3+z4+z7+z8z1z2z5z6)/4 ],
>> [ (x2+x4+x6+x8x1x3x5x7)/4 ,
>> (y2+y4+y6+y8y1y3y5y7)/4 ,
>> (z2+z4+z6+z8z1z3z5z7)/4 ]]);
>> det(A);
Alas, this results in something that isn't much fun anymore (:
DET = 1/16 * sum of +/ 192 (!) products
The straightforward numerical procedure, employed in SUNA56, is certainly more
efficient than brainlessly using such a CAS result.
So far so good. I recently (re)invented that Numerical differentiation can also
be carried out with triangles instead of quadrilaterals, and tetrahedra instead
of bricks. Again, the potential function takes it's values at the vertices, and
velocitycomponents are defined at the barycenter. At first sight, it seems not
to work: three equations with two unknows in the 2D case. But there is still a
possibility that one of the equations is dependent upon the other two. Indeed:
[ x1(x2+x3)/2 ].u + [ y1(y2+y3)/2 ].v = P1(P2+P3)/2
[ x2(x3+x1)/2 ].u + [ y2(y3+y1)/2 ].v = P2(P3+P1)/2
 +
minus { [ x3(x1+x2)/2 ].u + [ y3(y1+y2)/2 ].v = P3(P1+P2)/2 }
Hence, two of these equations are necessary and sufficient. Using MAPLE:
> > A := array([[ x1(x2+x3)/2 , y1(y2+y3)/2 ],
> > [ x2(x3+x1)/2 , y2(y3+y1)/2 ]]);
> > det(A);
> 3/4 x1 y2  3/4 x1 y3 + 3/4 x2 y3  3/4 x2 y1  3/4 x3 y2 + 3/4 x3 y1
As expected, there exists a beautiful symmetry in the (x,y) coordinates and the
numbering of the vertices: it shouldn't matter which two of the three equations
are used. Employing Cramer's rule, we find the solution:
DET = 3/4.(x1.y2  x1.y3 + x2.y3  x2.y1  x3.y2 + x3.y1)
u = 3/4.(P1.y2  P1.y3 + P2.y3  P2.y1  P3.y2 + P3.y1) / DET
v = 3/4.(x1.P2  x1.P3 + x2.P3  x2.P1  x3.P2 + x3.P1) / DET
This exercise has to be repeated for three dimensions: a tetrahedron.
> > A := array([[ x1(x2+x3+x4)/3 , y1(y2+y3+y4)/3 , z1(z2+z3+z4)/3 ],
> > [ x2(x3+x4+x1)/3 , y2(y3+y4+y1)/3 , z2(z3+z4+z1)/3 ],
> > [ x3(x4+x1+x2)/3 , y3(y4+y1+y2)/3 , z3(z4+z1+z2)/3 ]]);
> > det(A);
>
> 16 16 16 16 16
>  x2 y4 z3 +  x2 y3 z1   x2 y3 z4   x4 y2 z3   x2 y1 z3
> 27 27 27 27 27
[ ... stuff deleted. Why doesn't MAPLE bring the 16/27 outside parentheses? ]
The infamous "vi" editor was used for cleaning up the CAS result, resulting in:
DET = 16/27 *
( x2*y4*z3 + x2*y3*z1  x2*y3*z4  x4*y2*z3  x2*y1*z3 + x1*y2*z3
 x1*y2*z4 + x1*y3*z4  x1*y3*z2  x1*y4*z3 + x1*y4*z2 + x3*y2*z4
 x3*y2*z1  x3*y4*z2 + x3*y1*z2 + x4*y2*z1  x4*y3*z1 + x4*y3*z2
+ x4*y1*z3  x4*y1*z2  x2*y4*z1 + x2*y1*z4 + x3*y4*z1  x3*y1*z4 )
Cramer's rule is realized by making the following substitutions, with "vi" :
"s/x/P/g" for u , "s/y/P/g" for v , "s/z/P/g" for w . OK, let's doit:
u = 16/27 *
( P2*y4*z3 + P2*y3*z1  P2*y3*z4  P4*y2*z3  P2*y1*z3 + P1*y2*z3
 P1*y2*z4 + P1*y3*z4  P1*y3*z2  P1*y4*z3 + P1*y4*z2 + P3*y2*z4
 P3*y2*z1  P3*y4*z2 + P3*y1*z2 + P4*y2*z1  P4*y3*z1 + P4*y3*z2
+ P4*y1*z3  P4*y1*z2  P2*y4*z1 + P2*y1*z4 + P3*y4*z1  P3*y1*z4 ) / DET
v = 16/27 *
( x2*P4*z3 + x2*P3*z1  x2*P3*z4  x4*P2*z3  x2*P1*z3 + x1*P2*z3
 x1*P2*z4 + x1*P3*z4  x1*P3*z2  x1*P4*z3 + x1*P4*z2 + x3*P2*z4
 x3*P2*z1  x3*P4*z2 + x3*P1*z2 + x4*P2*z1  x4*P3*z1 + x4*P3*z2
+ x4*P1*z3  x4*P1*z2  x2*P4*z1 + x2*P1*z4 + x3*P4*z1  x3*P1*z4 ) / DET
w = 16/27 *
( x2*y4*P3 + x2*y3*P1  x2*y3*P4  x4*y2*P3  x2*y1*P3 + x1*y2*P3
 x1*y2*P4 + x1*y3*P4  x1*y3*P2  x1*y4*P3 + x1*y4*P2 + x3*y2*P4
 x3*y2*P1  x3*y4*P2 + x3*y1*P2 + x4*y2*P1  x4*y3*P1 + x4*y3*P2
+ x4*y1*P3  x4*y1*P2  x2*y4*P1 + x2*y1*P4 + x3*y4*P1  x3*y1*P4 ) / DET
Hmm, little hard work involved here. Sounds much like composing House music ;)

* Han de Bruijn; Applications&Graphics  "A little bit of Physics * No
* TUD Computing Centre; P.O. Box 354  would be NO idleness in * Oil
* 2600 AJ Delft; The Netherlands.  Mathematics" (HdB). * for
* Email: Han.deBruijn@RC.TUDelft.NL  Fax: +31 15 78 37 87 * Blood