Iteration exercises: f(x)=x^2 - 0.5 ; Fixpoint-irritation...

[Gottfried]

Ah, so I assumed well.    

Ok, here is the solution for the upper fixpoint:

Code:

fz2=x^2+(1+sqrt(3))*x       \\ using upper fixpoint  
F = mkCarleman ( polcoeffs(fz2,32) )  \\ carleman matrix of size 32x32 (power series will have 32 coeffs)
F_eigen = tri_eigen(F) ;
Ser(F_Eigen[4][,2])+O(x^6)   \\ inverse Schröder funktion at upper fixpoint
 -----------------
%82 = x + 0.2113248654*x^2 + 0.0239322565748*x^3 + 0.00174634543176*x^4 + 0.00009103437192*x^5 + O(x^6)

Added: here are some basic Pari/GP procedures used so far:

Code:

{ polcoeffs(lp, flgd=-2, flgNoConst=0) = local(llp, lpd, lv, lv1, maxd, oldps, s);
      maxd=if(flgd==-1|flgd==-2,n,flgd);
            oldps=default(seriesprecision,,1);  \\ save old series-precision and adjust for current query
            s=default(seriesprecision,maxd,1);

      llp=Pol(lp);lpd=poldegree(llp);    \\ convert input into formal POL
            s=default(seriesprecision,oldps,1);  \\ reset old series precision

      if(flgd==-2,maxd=max(0,lpd+1));
      lv=vector(maxd);
           for(k=0,min(lpd,maxd-1),lv[1+k]=polcoeff(llp,k));
      if(flgNoConst,lv[1]=0);      \\ make a true ZERO if constant should vanish

      return(lv); }

\\ polcoeffs(lp, flgd=-2,flgNoConst=0) : uses a scalar entry containing a
\\ polynomial, converts it into a vector of coefficients (calling implicitely "polcoeff()")
\\ flgd=-2: return poldegree+1 coefficients
\\ flgd=-1: return default (=n) coefficients
\\ flgd>=0: return flgd coefficients.


{ mkCarleman(lvec, serprec=n) = local(dim=length(lvec), tmpv, tmpser, tmpM, oldserprec);
              oldserprec=default(seriesprecision,,1);
              if(serprec<>oldserprec,default(seriesprecision,serprec,1));

     tmpv = Ser(lvec);
     tmpser = tmpv + O(x^dim);
     tmpM = matid(dim);
        for(c=2,dim,
             tmpM[,c]=polcoeffs(tmpser,dim)~ ;
             tmpser = tmpser*tmpv  + O(x^dim)
           );
     if(serprec<>oldserprec,default(seriesprecision,oldserprec,1));
     return(tmpM); }

\\ mkCarleman (vector-of-coefficients, seriesprecision = n)
\\ creates from a vector of coefficients (assumed to be of the power series of some function f(x))
\\ a matrix of the Carleman-type (transposed to the wikipedia-convention
\\ serprec : if dimension should other than order of supplied polynomial, give it in "serprec"


{ tri_eigen(U, dim=9999, numfmt=1) = local(u=U[2,2], UEW, UEWi, t=exp(u), uv);
       dim=max(1,min(rows(U),dim));

       uv=vectorv(dim,r,U[r,r]);
       UEW=numfmt*matid(dim);  \\ is also the Carlemanmatrix for Schröder-function
           for(c=1,dim-1,
                for(r=c+1,dim,
                      UEW[r,c]=sum(k=c,r-1,U[r,k]*UEW[k,c])/(uv[c]-uv[r])
              ));

       UEWi=numfmt*matid(dim); \\ UEWi shall become the inverse of UEW, is also the
                                      \\ Carleman-matrix of the inverse Schröder-function
       for(r=2,dim,
            forstep(c=r-1,1,-1,
                 UEWi[r,c]=sum(k=0,r-1-c,U[r-k,c]*UEWi[r,r-k])/(uv[r]-uv[c])
           ));

       return([[u,t,exp(u/t)],UEW,uv,UEWi]);}

\\ tri_eigen(U, dim=9999, numfmt=1) - performs diagonalization of lower triangular matrices
\\ if inputmatrix is of Carlemantype, then the Eigenmatrices are of the Carlemantype as well
\\ returns vector of [coeffs, M, D, M^-1] M being eigenvectors, D being diagonalmatrix of eigenvalues
\\                                        coeffs being [u=log(t), t, u/t = log(b)] where b=t^(1/t)
\\    diagonalization can also be performed if input matrix is Carleman and has parameter u only symbolic
\\
\\ dim : if 9999 use size of input matrix
\\     : if smaller use size of top-left sumbatrix
\\
\\ numfmt : 1 : return integer values if possible ; sometimes faster, sometimes slower than with real values
\\          1.0 : compute in float format (with stadard-precision in Pari/GP). If large factorials
\\                occur in large dimension matrices Ut then using numfmt=1.0 prevents gigantic integers
\\                or rational numbers with gigantic components.

{ fmt(realprec=-1, dispdigits=-1, flg=1) = local(oldp, oldf);
     oldp=precision(0.0);
     oldf=default(format,,1);
     if(realprec>-1,default(realprecision,realprec,1);oldp=realprec);
     if(dispdigits>-1,oldf=Str("g0.",dispdigits));
     oldf=default(format,oldf,1);
     if(flg,print("prec=",oldp," display=",oldf)); }
\\ fmt(200,12) sets Pari/GP real precision to 200 and display digits to 12
\\ if only one parameter should be changed, the other stays untouched. fmt(,8) ; fmt(400) etc


\\ additional remarks:
\\  I use the global variable "n" for default size for matrices and vectors. (This is "n"
\\  simply I'm an old statistician... )
\\  By default, I set in my session-initialization n=32 . With this size all standard
\\  matrix-constants are also precomputed to be permanently available in the session.
\\  If "n" should be changed either then the session-initialization should
\\  be recalled so the matrix-constants are recomputed too (if they are needed at all).
\\
\\  I've implemented the short routine "fmt(internal-precision,display-precision)" and
\\  initialize sessions with "fmt(200,12)" to have by default 200 dec digits precision for
\\  real-arithmetic. This can be changed at any time, but remember that possibly used
\\  real constant (like bases for exponentiation, computed fixpoints etc) should then
\\  recomputed as well.
\\
\\  The full set of routines for my (dynamical) matrix-toolbox-procedures is too large
\\  to be added here, I might send it on personal request.

Gottfried