[Update] Comparision of 5 methods of interpolation to continuous tetration

[Gottfried]

Hi Sheldon,

first to your second question, number 1): yes, the function is real on the real axis because we use just the eigendecomposition of a real matrix - no complex coefficients anywhere.
About number 2) I've no idea yet, perhaps its too late for me tonight...
I'll come back to this tomorrow.

We do not have a true Schröder- and inverse Schröder-function because the diagonal after the diagonalization (the set of eigenvalues) does not contain the consecutive powers of the log of the fixpoint but are more or less arbitrary - depending on the size of the matrix. So I don't see at the moment how to model the sexp() and inverse sexp()-functionality.


Gottfried

Here is code to reproduce the behaviour - no optimization is made, just the basic functionality. If you want to exceed the range of convergence of the "power series" you should add explicite iteration over the integer-part of the heights.

Code:

n=24   \\ set size for matrices, note this needs ~ 200 digits real precision,
       \\ n=64 needs about 4000 digits precision and 2 hours to be computed!
default(precision,400) \\ when n~ 24

\\ procedures ===========================================================
\\    several global variables/matrices will be created/changed
{init(base,matsize,prec) = local(a);
   n=matsize;  \\ set size for matrices, note size=24x24 needs ~ 200 digits real precision,
               \\ n=64 needs about 4000 digits precision and 2 hours to be computed!
   default(precision,prec,1); \\ when n~ 24 then prec ~ 200



      b = base;bl=log(b); \\ sets the base for tetration x-> 4^x
                    \\ bl must have the required precision depending on matrix-size n
                    \\ thus must be recomputed if size and precision is increased!

      B = matrix(n,n,r,c,((c-1)*bl)^(r-1)/(r-1)!) ;
                 \\ the Carlemanmatrix for x-> b^x
                 \\ the resulting coefficients for the power-series
                 \\ are in the second column

     \\ Diagonalization such that B = M * D * W (where W=M^-1)
      M = mateigen(B);
      W = M^-1 ;
      D = diag(W * B * M);

      POLY= M * matdiagonal(W[,2])  \\ this matrix is a "kernel" which needs only the parameter
                    \\ for the height to give the associated power series

      print("matrices are created");
   return(1);}


f(x) = sum(k=0,n-1,x^k * B[1+k,2])  \\ just to show how this works for the
                                    \\ function f(x) itself



        \\ then we can the h'th fractional power  of Bb approximate by M * D^h * W (= B^h)
         make_powerseries(h) = return( POLY * vectorv(n,r, D[r]^h) );

{ tet_polynomial(x,h) = local(pc);
       pc = POLY*vectorv(n, r , D[r]^h ) ;  \\ creates coefficients for power series
       result = sum(k=0,n-1, x^k * pc[1+k]);  \\ evaluates the powerseries
       return(result);}



\\ note that M and W do not really implement the Schröder-functions because
\\ the entries in D are not consecutive powers of one argument


\\ then you start:=================================================================

init(4,24,200) \\ for matrixsize 24x24 you'll need at least 200 digits precision

x0 = 0.1*I
x05= tet_polynomial( x0 , 0.5 )   \\ halfiterate from x0=0.1*I
x1 = tet_polynomial( x05, 0.5 )   \\ halfiterate from x05 should equal one whole iteration:

x1 - b^x0    \\ check the difference  

\\ bigger matrixsize
init(4,32,600) \\ for matrixsize 32x32 you'll need at least 600 digits precision

x0 = 0.1*I
x05= tet_polynomial( x0 , 0.5 )   \\ halfiterate from x0=0.1*I
x1 = tet_polynomial( x05, 0.5 )   \\ halfiterate from x05 should equal one whole iteration:

x1 - b^x0    \\ check the difference  


\\ bigger matrixsize
init(4,64,4000) \\ for matrixsize 64x64 you'll need at least 4000 digits precision
                 \\ this needs 2 hours for computation!
x0 = 0.1*I
x05= tet_polynomial( x0 , 0.5 )   \\ halfiterate from x0=0.1*I
x1 = tet_polynomial( x05, 0.5 )   \\ halfiterate from x05 should equal one whole iteration:

x1 - b^x0    \\ check the difference