Numerical Analysis

Archive for November, 2009

Sun-specific BLAS??

Does such a beast exist?


_______________________________________________________________________________
Ajay Shah, (213)734-3930, ajays…@usc.edu
                             The more things change, the more they stay insane.
_______________________________________________________________________________

.
posted by admin in Uncategorized and have Comments (6)

Looking for combination algorithm

Dear MathNetters:

I’m looking for an algorithm to generate all possible _combinations_ for
an number of elements (n) taken m at a time, that is n choose m.  There are
lots of permutation algorithms, but I don’t see combinations.
Any suggestions in the literature or code frags greatly appreciated.

Cheers,
David S. McCormick

MIT-EAPS Geology
Cambridge, MA  02139
d…@athena.mit.edu
617-253-5747

posted by admin in Uncategorized and have No Comments

Wanted: integration method for 1/sqrt(f(x)) dx

Hi folks
i would need a proper integration method to solve the following definite
integral: 1/sqrt(a+b*cos(x)+c*sin(x)) dx.

Please send any replies to mich…@hsa.e-technik.tu-muenchen.de
 Thanks a lot
michael

posted by admin in Uncategorized and have No Comments

Request for problems in computing….

I am planning a seminar for computer science grad students in which I want
to focus on programming languages, compiling, and numerical processes.
I want to include lots of examples of pathologies (like cancellation) and
no-no’s (don’t apply algebraic rules when optimizing floating point).

I would appreciate any examples that the readers are aware of which fit
into these categories, or other categories that you think should be emphasized
in such a seminar. The more insidious, the better for examples. Flames accepted
as long as you don’t get personal :-)

steve

===============================================================================
Steve (really "D. E.") Stevenson           st…@hubcap.clemson.edu
Department of Computer Science,            (803)656-5880.mabell
Clemson University, Clemson, SC 29634-1906

posted by admin in Uncategorized and have No Comments

SUNA, Triangle Isoparametrics

Isoparametric transformation is the standard approach where the Finite Element
Method relies on when it has to deal with curved geometries: it’s strong point.
But let’s consider first the simplest non-trivial finite element shape in two
dimensions, the linear triangle, and see where isoparametric transformations
                       do come from. Function behaviour is approximated inside
          3            such a triangle by a _linear_ interpolation between the
        /   \          function values at the vertices, or nodal points.
      /       \        Let  T  be such a function, and  x,y  coordinates,
    /           \      then:
  1 ————- 2               T = A.x + B.y + C

Where the constants  A, B, C  are determined by:

      T1 = A.x1 + B.y1 + C
      T2 = …

But wait. The first of these equations can already be used to eliminate the
constant  C , once and forever:

      T – T1 = A.(x – x1) + B.(y – y1)

Then the constants  A  and  B  are determined by:

      T2 – T1 = A.(x2 – x1) + B.(y2 – y1)
      T3 – T1 = A.(x3 – x1) + B.(y3 – y1)

Two equations with two unknowns. The solution is found by Cramer’s rule:

      A = [ (y3 - y1).(T2 - T1) - (y2 - y1).(T3 - T1) ]/J
      B = [ (x2 - x1).(T3 - T1) - (x3 - x1).(T2 - T1) ]/J

Where:  J = (x2 – x1).(y3 – y1) – (x3 – x1).(y2 – y1) = 2 * triangle area .
So:
         T – T1 = h.(T2 – T1) + k.(T3 – T1)
Where:
              h = [ (y3 - y1).(x - x1) - (x3 - x1).(y - y1) ]/J
              k = [ (x2 - x1).(y - y1) - (y2 - y1).(x - x1) ]/J
Or:
          | h |   | + (y3 – y1)    - (x3 – x1) |   | x – x1 |
          |   | = |                            |/J |        |
          | k |   | – (y2 -y1)     + (x2 – x1) |   | y – y1 |

The inverse of the following problem is recognized herein:

     | x – x1 |   | x2 – x1     x3 – x1 | | h |
     |        | = |                     | |   |
     | y – y1 |   | y2 – y1     y3 – y1 | | k |
Or:
        x – x1 = h.(x2 – x1) + k.(x3 – x1)
        y – y1 = h.(y2 – y1) + k.(y3 – y1)
But:
        T – T1 = h.(T2 – T1) + k.(T3 – T1)

Typical! The _same_ expression holds for the function  T  as well as for the
coordinates x and y . Which is precisely what people mean by an ISOPARAMETRIC
("same parameters") transformation.

Now recall the formulas which express h and k into x and y:

    h = [ (y3 - y1).(x - x1) - (x3 - x1).(y - y1) ]/J
    k = [ (x2 - x1).(y - y1) - (y2 - y1).(x - x1) ]/J

Thus h can be interpreted as: area of the sub-triangle spanned by the vectors
(x – x1 , y – y1) and (x3 – x1 , y3 – y1) divided by the whole triangle area.
And k can be interpreted as: area of the sub-triangle spanned by the vectors
(x – x1 , y – y1) and (x2 – x1 , y2 – y1) divided by the whole triangle area.
This is the reason why  h  and  k  are sometimes called _area-coordinates_.

There are even _three_ of these coordinates in litterature (: Zienkiewicz).
But one of them can be allways be safely discarded. I think it’s a bad habit
to retain a "coordinate" like (1 – h – k), being dependent on the other two.

Instead of area-coordinates, we prefer to talk about _local coordinates_
h and k  of an element, in contrast to the _global coordinates_  x and y .

 3              It is possible that local coordinates coincide with the
 |\             global coordinates. A triangle for which this is the case
 |  \           is called a _parent element_. The portrait of the parent
 |    \         triangle is displayed here: rectangular, two sides equal.
 |______\
 1       2

Let’s reconsider for a moment the expression:

     T – T1 = h.(T2 – T1) + k.(T3 – T1)

Partial differentiation to  h and k  gives:

        dT/dh = T2 – T1
        dT/dk = T3 – T1
So:
        T = T1 + h.dT/dh + k.dT/dk

This is part of a Taylor series expansion around node (1). Such Taylor series
expansions are very common in Finite Difference analysis.
Now rewrite as follows:

     T = (1 – h – k).T1 + h.T2 + k.T3

Here the functions  (1 – h – k) , h , k  are called the SHAPE FUNCTIONS of the
Finite Element. Shape functions  Nk  have the property that they are unity in
one of the nodes (k), and zero in all other nodes. In our case:

     N1 = 1 – h – k   ;   N2 = h   ;   N3 = k

So we have two representations, which are allmost trivially equivalent:

     T = T1 + h.(T2 – T1) + k.(T3 – T1)      : Finite Difference like
     T = (1 – h – k).T1 + h.T2 + k.T3        : Finite Element like

The two sets of functions (1,h,k) and (1-h-k,h,k) are mutually related as
follows:
             F.D. <- F.E.                             F.E <- F.D.
    | 1 |   | 1   1   1 | | 1-h-k |      | 1-h-k |   | 1  -1  -1 | | 1 |
    | h | = | 0   1   0 | | h     |      | h     | = | 0   1   0 | | h |
    | k |   | 0   0   1 | | k     |      | k     |   | 0   0   1 | | k |

The two sets of functions can be conceived as two coordinate systems in an
abstract three dimensional space. Then the above matrices are _transition_
matrices, belonging to F.D. <-> F.E. base transitions. Yes, the above is an
instance of the Unification Principle (Re: SUNA, Manifesto).

Not too difficult, eh? Well, not yet. To be continued …

* 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.deBru…@RC.TUDelft.NL –| Fax: +31 15 78 37 87 —-* Blood

posted by admin in Uncategorized and have No Comments

There must be a simple answer to this!

   Hi!

   I’m secretary of our bowling league and at the beginning of the season
I have to make a seasonal sheet which lists which teams play who each week.
   The way our league works, for each season we have 4 cycles.  A cycle is
defined by the number of weeks it takes a team to play every other team.
So we have 4 cycles and the season lasts for 28 weeks(we have 8 teams).
   Now making up this seasonal sheet is easy.  There are pre-made sheets
which look something like this:

     Week Num     Lane 1-2     Lane 3-4     Lane 5-6     Lane 7-8
     ————————————————————
        1          5-6          7-8          3-4           1-2
        2          2-4          1-3          5-7           6-8
        3          7-3          2-6          1-8           5-4
        ..
        28

We then have team names which correspond to the numbers 1-8.

There are 2 rules:

   1> A team can bowl another team only once in a cycle.
   2> A team can bowl on the same lane only twice in one cycle, and the
       second week should be as many weeks away from the first as possible.
       (There will be one team which will only bowl once on the same lane.)

   I’d like to know what algorithm will generate the table.  Is it a
permutation of some sort?

                                   Any help would be appreciated!

                                                                [Tim]

posted by admin in Uncategorized and have No Comments

Introductory text to stochastic differential equations

Can someone recommend a good introductory text on stochastic differential
equations, the Feynman-Kac formula and its connections to the theory of
operator semigroups? It is no drawback if it is application oriented.


"Kleinphi macht auch Mist"
Reiner Lenz | Dept. EE.                 |
            | Linkoeping University     | email:        rei…@isy.liu.se
            | S-58183 Linkoeping/Sweden |

posted by admin in Uncategorized and have No Comments

SUNA, Convex Quadrilaterals

Let’s consider next the simplest-but-one element shape in 2-D: a quadrilateral.
                        Function behaviour may be approximated inside such a
          3             quadrilateral by a _bilinear_ interpolation between
        /   \           the function values at the vertices, or nodal points.
      /       \         Let  T  be such a function, and  x,y  coordinates,
    /           \       then:
  1      0,0      4               T = A + B.x + C.y + D.x.y        <–there–\
    \           /                                                            |
      \       /         As an example, consider the quadrilateral depicted –/
        \   /           with nodes given by the 2′nd and 3′th row of the matrix
          2             that originates when specifying  T  for the nodes:

| T1 |   | 1   -1    0    0 | | A |      The last column of our matrix is zero.
| T2 | = | 1    0   -1    0 | | B |      Thus A,B,C and D cannot be found here!
| T3 |   | 1    0   +1    0 | | C |      This means that a method as has been
| T4 |   | 1   +1    0    0 | | D |      used for the linear triangle can not
           1    x    y   x.y             be employed for higher order elements.

There is, however, another way out. Let’s assume that the same expression holds
for the function  T  as well as for the coordinates x and y : an ISOPARAMETRIC
transformation again. The next step is to find a good _parent_ element. Well,
the above one surely did’nt work out. Let’s try the following:

   3———4       T = A + B.h + C.k + D.h.k     gives:
   |         |
   |         |       | T1 |   | 1  0  0  0 | | A |
   |         |       | T2 | = | 1  1  0  0 | | B |       F.E. <- F.D.
   |         |       | T3 |   | 1  0  1  0 | | C |       transition
   1———2       | T4 |   | 1  1  1  1 | | D |
                                   h  k
Solving for
A,B,C,D          | A |   | +1   0   0   0 | | T1 |
isn’t too        | B | = | -1  +1   0   0 | | T2 |       F.D. <- F.E.
difficult:       | C |   | -1   0  +1   0 | | T3 |       transition
                 | D |   | +1  -1  -1  +1 | | T4 |

So we have two equivalent representations:

T = T1 + (T2 – T1).h + (T3 – T1).k + (T1 – T2 – T3 + T4).h.k    : F.D.-like
T = (1 – h).(1 – k).T1 + h.(1 – k).T2 + (1 – h).k.T3 + h.k.T4   : F.E.-like

Partial derivatives are:

        dT/dh = T2 – T1
        dT/dk = T3 – T1
and:
        d^2(T)/dh.dk = T1 – T2 – T3 + T4

Shape functions:   (1 – h).(1 – k)  ;  h.(1 – k)  ;  (1 – h).k  ;  h.k
 at nodal point:          1             2                   3       4

Global coordinates are related to local coordinates inside a general curved
quadrilateral by the isoparametric transformation:

   x = (1 – h).(1 – k).x1 + h.(1 – k).x2 + (1 – h).k.x3 + h.k.x4
   y = (1 – h).(1 – k).y1 + h.(1 – k).y2 + (1 – h).k.y3 + h.k.y4

The differential form of an isoparametric transformation is often even more
important (here d*/d* are partial derivatives):

   dx = dx/dh.dh + dx/dk.dk            | dx |   | dx/dh  dx/dk | | dh |
   dy = dy/dh.dh + dy/dk.dk     or:    |    | = |              | |    |
                                       | dy |   | dy/dh  dy/dk | | dh |
Inverse differential transform:

 | dh |   |  dy/dk   -dx/dk |     | dx |
 |    | = |                 | / J |    |
 | dk |   | -dy/dh    dx/dh |     | dx |

The inverse transform exists if and only if the Jacobian determinant is never
zero. Because  J  will be a function of (h,k) this also means that it should
not switch its sign. We shall demand that it is positive:

         J = dx/dh.dy/dk – dx/dk.dy/dh > 0
Here:
        dx/dh = (1 – k).(x2 – x1) + k.(x4 – x3)
        dy/dh = (1 – k).(y2 – y1) + k.(y4 – y3)
        dx/dk = (1 – h).(x3 – x1) + h.(x4 – x2)
        dy/dk = (1 – h).(y3 – y1) + h.(y4 – y2)

Substitution gives:

  J = (1 – h).(1 – k).J1 + h.(1 – k).J2 + (1 – h).k.J3 + h.k.J4

Where:    J1 = (x2 – x1).(y3 – y1) – (x3 – x1).(y2 – y1)
          J2 = (x2 – x1).(y4 – y2) – (x4 – x2).(y2 – y1)
          J3 = (x4 – x3).(y3 – y1) – (x3 – x1).(y4 – y3)
          J4 = (x4 – x3).(y4 – y2) – (x4 – x2).(y4 – y3)

Therefore the Jacobian is expressed by the quadrilateral shape functions into
its values  Jk  at the nodal points, as it should be.
Well, not really. Look at the Finite Difference representation:

   J = J1 + (J2 – J1).h + (J3 – J1).k + (J1 – J2 – J3 + J4).h.k

The last term is zero. This is easiest seen by considering the fact that the
area of the quadrilateral is equal to the sum of the areas of the triangles
(1) and (4), but also to the sum of the areas of the triangles (2) and (3).
Therefore  J1 + J4 = J2 + J3  which proves the above claim. We may conclude
| that  J  is not a bi-linear but only a _linear_ expression in the  Jk ‘s.

Now it was demanded that the Jacobian must be positive over the whole element.
This can only be the case if all its values at the vertices are postive.
Let’s take vertex (1) as an example. Here:   J1 = |_12_|.|_13_|.sin(phi) .
_1x_ is the vector which points from 1 to (x), (phi) is the angle between the
vectors _12_ and _13_. This means that  sin(phi) > 0  or  0 < phi < Pi .
Thus obtuse angles are forbidden. Quadrilateral elements must be convex.

Discliamer: the above is still well known. To be continued …

* 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.deBru…@RC.TUDelft.NL –| Fax: +31 15 78 37 87 —-* Blood

posted by admin in Uncategorized and have No Comments

SUNA, Five Point Star

Consider the well known five point star. Function values (f1,f2,f3,f4,f5) are
                 defined at its nodal points. These values can be interpolated
       5         by a (Finite Element) polynomial, which is defined by:
       |
       |             f = a1 + a2.x + a3.x.x + a4.y + a5.y.y
 2—–1—–3
       |         Specify for the nodes: (0,0), (-1,0), (+1,0), (0,-1), (0,+1).
       |         This results in 5 equations for the unknowns (a1,a2,a3,a4,a5):
       4
                     f1 = a1  ;  f2 = a1 – a2 + a3  ;  f3 = a1 + a2 + a3
                                 f4 = a1 – a4 + a5  ;  f5 = a1 + a4 + a5

From which it follows that:  a1 = f1  ;  a2 = (f3 – f2)/2  ;  a4 = (f5 – f4)/2
                            a3 = (f3 – 2.f1 + f2)/2  ;  a5 = (f5 – 2.f1 + f4)/2

These are well known finite difference schemes for the zero’th, first and second
partial derivatives on the five node star.
By substitution of the a’s, the function  f  is expressed as follows:

              f = f1 + (f3 – f2)/2.x + (f3 – 2.f1 + f2)/2.x.x +
                     + (f5 – f4)/2.y + (f5 – 2.f1 + f4)/2.y.y
Rewrite this:
              f = (1 – x.x – y.y).f1 + 1/2.(- x + x.x).f2 + 1/2.(+ x + x.x).f3
                                     + 1/2.(- y + y.y).f4 + 1/2.(+ y + y.y).f5

Sic! Here we find the (Finite Element) Shape Functions of the five point star.

They are:     N1(x,y) = 1 – x.x – y.y
              N2(x,y) = 1/2.(- x + x.x) ; N3(x,y) = 1/2.(+ x + x.x)
              N4(x,y) = 1/2.(- y + y.y) ; N5(x,y) = 1/2.(+ y + y.y)

It is necessary to introduce now local coordinates (h,k) instead of (x,y).
Then, a (Finite Element) Isoparametric Mapping can be defined as follows:

   x = N1(h,k).x1 + N2(h,k).x2 + N3(h,k).x3 + N4(h,k).x4 + N5(h,k).x5
   y = N1(h,k).y1 + N2(h,k).y2 + N3(h,k).y3 + N4(h,k).y4 + N5(h,k).y5

The partial derivatives  dx/dh, dx/dk, dy/dh, dy/dk  are calculated:

    dx/dh = (1/2 – h).(x1 – x2) + (1/2 + h).(x3 – x1)
    dy/dh = (1/2 – h).(y1 – y2) + (1/2 + h).(y3 – y1)
    dx/dk = (1/2 – k).(x1 – x4) + (1/2 + k).(x5 – x1)
    dy/dk = (1/2 – k).(y1 – y4) + (1/2 + k).(y5 – y1)

Jacobian determinant:  J = (1/2 – h)(1/2 – k).J1 + (1/2 + h)(1/2 – k).J2
                         + (1/2 – h)(1/2 + k).J3 + (1/2 + h)(1/2 + k).J4

Where:   J1 = (x2 – x1).(y4 – y1) – (x4- x1).(y2 – y1)
         J2 = (x4 – x1).(y3 – y1) – (x3- x1).(y4 – y1)
         J3 = (x5 – x1).(y2 – y1) – (x2- x1).(y5 – y1)
         J4 = (x3 – x1).(y5 – y1) – (x5 -x1).(y3 – y1)

The shape functions of the five point star’s Jacobian are therefore identical
to those of a (Finite Element) quadrilateral: see for example O.C. Zienkiewicz’
"The Finite Element Method" or the previous poster on "Convex quadrilaterals".
This quadrilateral defines an area inside the star as indicated in the figure:

              5              These area’s "happen" to be ADJACENT when the F.D.
              |              stars are combined into a (curvilinear) grid.
        J3____|____J4
         |    |    |         The determinants  J1,J2,J3,J4  denote the area’s of
    2____|____1____|____3    triangles which are spanned by the star’s "arms".
         |    |    |         For the *inverse* of an isoparametric mapping to
         |____|____|         exist, it is required that the Jacobian J has a
        J1    |    J2        uniform (positive) sign. Since J is in general a
              |              bilinear interpolation of the Jk’s (k = 1,2,3,4),
              4              it is necessary and sufficient that the latter are
                             positive:  J1 > 0 , J2 > 0 , J3 > 0 , J4 > 0 .
The geometrical meaning of this is that a five point star must not be distorted
in such a rigourous way that its arms "overcross" each other.

Performing the isoparametrics for an arbitrary function f (partial derivatives):

   df/dx = dh/dx.df/dh + dk/dx.df/dk = (dy/dk.df/dh – dx/dk.df/dk) / J
   df/dy = dh/dy.df/dh + dk/dy.df/dk = (- dy/dh.df/dh + dy/dk.df/dk) / J
   The derivatives  d(x,y)/d(h,k)  and  J  were calculated above.

If these formulas are specified for the special case that (h,k) = (+1/2,+1/2)
or any other of the "inscribed" quadrilateral nodes, then they reduce to well
known (Finite Element) schemes for taking derivatives at a triangular element:
Re: SUNA, Triangle isoparametrics.

This means that a devising a special Finite Element Method for five point stars
in a curvilinear grid is rather pointless: such a method would NOT behave very
differently from a conventional F.E. method on ordinary triangles! I feel that
the above is already a nice demonstration of the Unification Principle at work
(Re: SUNA, The Manifesto).

Disclaimer: the above stuff was published earlier in this group. I nevertheless
decided to post it, in order to preserve continuity with past & future.

To be continued …

* 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.deBru…@RC.TUDelft.NL –| Fax: +31 15 78 37 87 —-* Blood

posted by admin in Uncategorized and have No Comments

Symbolic algebraic manipulation algorithm

If this has been discussed before please email, otherwise please post.
I am interested in trying to write a simple (symbolic) algebraic manipulator.
Also some symbolic differentiation and integration is of interest. I want to
write this in C++. Can anyone point me to books, papers, or other references
on this subject. Many thanks.

Bob Kaires

posted by admin in Uncategorized and have No Comments