sci.math.num-analysis 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 two-dimensional 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) 1981-1990 by the University of Waterloo. > \ MAPLE / All rights reserved. MAPLE is a registered trademark of > <____ ____> Waterloo Maple Software. > >> with(linalg); >> A := array([[ (x2+x4-x1-x3)/2 , (y2+y4-y1-y3)/2 ], >> [ (x3+x4-x1-x2)/2 , (y3+y4-y1-y2)/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 three-dimensional case, which was significant for our theory in article (56). >> A := array([[ (x5+x6+x7+x8-x1-x2-x3-x4)/4 , >> (y5+y6+y7+y8-y1-y2-y3-y4)/4 , >> (z5+z6+z7+z8-z1-z2-z3-z4)/4 ], >> [ (x3+x4+x7+x8-x1-x2-x5-x6)/4 , >> (y3+y4+y7+y8-y1-y2-y5-y6)/4 , >> (z3+z4+z7+z8-z1-z2-z5-z6)/4 ], >> [ (x2+x4+x6+x8-x1-x3-x5-x7)/4 , >> (y2+y4+y6+y8-y1-y3-y5-y7)/4 , >> (z2+z4+z6+z8-z1-z3-z5-z7)/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 velocity-components are defined at the barycenter. At first sight, it seems not to work: three equations with two unknows in the 2-D 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 * E-mail: Han.deBruijn@RC.TUDelft.NL --| Fax: +31 15 78 37 87 ----* Blood