Cheta with base-change: preliminary results

[Gottfried]

Hi Jay,

there is a simple recursion formula, with which you can generate the inverse matrix.
I used the description at wikipedia for the stencils,

Code:

f(x-2h)= f(x) -2h f'(x)/1! + 4h^2 f''(x)/2! - 8h^3*f'''(x)/3! + 16h^4*f''''(x)/4!
f(x-h) = f(x) -1h f'(x)/1! + 1h^2 f''(x)/2! - 1h^3*f'''(x)/3! +  1h^4*f''''(x)/4!
f(x)   = f(x)  0             0                0                  0
f(x+h) = f(x) +1h f'(x)/1! + 1h^2 f''(x)/2! + 1h^3*f'''(x)/3! +  1h^4*f''''(x)/4!
f(x+2h)= f(x) +2h f'(x)/1! + 4h^2 f''(x)/2! + 8h^3*f'''(x)/3! + 16h^4*f''''(x)/4!

To prepare this for matrix-analysis I rewite this as matrix-equation:

Code:

´
|f(x-2h)|= |1 -2  4  -8 16 |  | f(x)/0!    |
|f(x-h) |= |1 -1  1  -1  1 |  | f'(x)/1!   |
|f(x)   |= |1  0  0   0  0 |* | f''(x)/2!  |
|f(x+h) |= |1  1  1   1  1 |  | f'''(x)/3! |
|f(x+2h)|= |1  2  4   8 16 |  | f''''(x)/4!|

and define the following matrices:

Code:

´
Fh(x) = columnvector of [f(x-2h),f(x-h),f(x),f(x+h),f(x+h)]
Fd(x) = columnvector of [f(x)/0!,f'(x)/1!,f''(x)/2!,f'''(x)/3!,f''''(x)/4!]

prefixing them with "d", if taken as diagonal matrices

Surely this definitions are to be adapted for higher dimension with the size-parameter, but it's easier here to stay with simple examples

Code:

´
Fac = diagonalmatrix of factorials diag(0!, 1!,2!,3!,4!)
FacI = Fac^-1 = diag(1/0!,1/1!,1/2!,1/3!,1/4!)

Sc    = matrix of the Stencils coefficients, such that

    Fh(x) = Sc * Fd(x)

What you want, is the inverse of Sc, such that Fd(x) = Sc^-1 * Fh(x) and evaluating the rhs you get the values for the derivatives f'(x)/1!, f''(x)/2! etc.
The computation of the inverse of a square matrix, especially if one is interested in different sizes or even for the limit case of infinite size is tedious or even impossible. But sometimes we can make life easier, if the square matrix can be factored in triangular and diagonal matrix-factors, and then the triangular factors are easier inverted.

Using the lower triangular pascalmatrix P and its inverse PI we find, that we can indeed factor Sc into such factors: Sc = P * Fac * X where we know the entries in the other factors P and Fac and we can find then X by X = FacI * PI * Sc

Now an independent description of entries of X would be nice, but is not obvious. On the other hand, the inversion of Sc is easier now, because the inversion of X, which is now of the simple triangular form, needs significantly less computation

Code:

´
     Sc^-1 = X^-1 * FacI * PI

So this alone would save computational effort.

But the really good news is, that I've found an independent description for the entries in X as well for X^-1 (call it XI, since we can compute it independently), so that the inversion can be omitted at all.

Assume an odd dimension for the problem dim=9, then use d2 = (dim-1)/2=4
Reference a cell as a(r,c), index beginning at zero. If an index r or c exceeds the matrix, let the value of a(r,c):=0. Work only for the upper triangular submatrix.

Code:

´
Init X and XI as identity-matrix
Further cells are (where we need only c>r)

   in X  :   a(r,c) = a(r-1,c-1) - (d2 -  r  )*a(r,c-1)  
   in XI :   a(r,c) = a(r-1,c-1) + (d2 -(c-1))*a(r,c-1)

Then

Code:

Fd(x) =     Sc^-1      * Fh(x)
          = XI * FacI * PI * Fh(x)

without additional inversions...

Data-examples: I add only a picture here - I'm lazy to retranslate the formatting from my winword-table-cells at moment

   

   

   

I can send you the Pari/Gp-code if this is needed.


Gottfried