¿La solución formal de la serie de potencias de $f(f(x))= \sin( x)$ ¿converger?

Question

He pasado algún tiempo utilizando gp-pari. Hay, por supuesto, un serie de potencia formal solución a $f(f(x)) = \sin x.$ Se muestra a continuación, identificado con el símbolo $g$ porque no estoy del todo seguro de que esté en función de algo.

Por otro lado, si los coeficientes siguen disminuyendo (en general), esto sugiere un radio de convergencia no nulo. Si el radio de convergencia es distinto de cero, entonces dentro de eso, no sólo es una función definida y, ya sabes, analítica, sino que la ecuación funcional se satisface. De hecho, todo lo que se necesita es un radio de convergencia estrictamente mayor que $\frac{\pi}{2}$ debido a ciertas simetrías. Por ejemplo, dado mi polinomio $g,$ parece que tenemos $g=1$ en torno a $x \approx 1.14.$ Entonces parece que tenemos un máximo local en $x =\frac{\pi}{2},$ y aparentemente hay $g \approx 1.14,$ estrictamente mayor que 1, lo cual es un punto importante. Así que todo caería en su lugar con un radio de convergencia suficientemente grande y no nulo.

$$\begin{array}{lll} g & = & x - \frac{x^3 }{ 12} - \frac{x^5 }{ 160} - \frac{53 x^7 }{ 40320} - \frac{23 x^9 }{71680} - \frac{92713 x^{11}}{1277337600} - \\\ & & \\\ & & \frac{742031 x^{13} }{79705866240} + \frac{594673187 x^{15} }{167382319104000} + \frac{329366540401 x^{17} }{91055981592576000} + \\\ & & \\\ & & \frac{104491760828591 x^{19} }{62282291409321984000} + \frac{1508486324285153 x^{21} }{4024394214140805120000} + \cdots \end{array}$$

Obsérvese que el polinomio $g$ es menor que $x$ pero más grande que $\sin x,$ para, digamos, $0 < x \leq \frac{\pi}{2}.$

Entonces, esa es la pregunta, ¿la serie de potencia formal que comienza con $g$ convergen en cualquier lugar que no sea $x = 0$ ?

EDIT: nota que los términos después de la inicial $x$ mismo han resultado ser todos $$\frac{a_{2 k + 3} x^{2 k + 3} }{2^k ( 2 k + 4)!}$$ donde cada $a_{2 k + 3}$ es un número entero. Esto parece demostrable, aunque todavía no lo he intentado.

EDIT, viernes 12 de noviembre de 2010. Ahora parece realmente improbable que este problema en particular tenga una respuesta analítica. Sospecho que la respuesta es $C^\infty$ y analítica a trozos, con fallo de analiticidad sólo en los puntos "parabólicos" donde la derivado tiene un valor absoluto tan grande como 1, siendo esos puntos $0,\pi, 2 \pi, \ldots.$ Sin embargo, necesitamos el punto de anclaje en el punto fijo 0, si no, ¿cómo empezar? Y creo que la serie de potencias servirá como expansión asintótica alrededor de 0.

Dado el problema del tamaño de la derivada, ahora espero grandes cosas, y una solución obviamente periódica y analítica, a la variante más fácil $f(f(x)) = g(x) = (1/2) \sin x.$ Me gustaría tanto una buena serie de potencias como una buena respuesta por métodos sumando iterados $g^{[k]}(x),$ que por el momento es un método totalmente misterioso para mí, pero atractivo para las funciones objetivo periódicas ya que la periodicidad sería automática.

Answer 1

Esta no es una nueva respuesta sino información adicional para la respuesta de Will Jagy sobre el cálculo de la función de Abel con el método de J. Écalle.

• [actualización] Con algunos ejemplos, se observa que el cálculo "estándar" de la iteración fraccionaria a través de la serie de potencias formal para el logaritmo iterativo y su truncamiento a los términos principales proporciona los mismos valores que el método de Écalle utilizando la función de Abel como se describe a continuación. Con el mismo truncamiento a 64 términos de la serie de potencias y el mismo desplazamiento de z_0 a z_h hacia cero el diferencia entre los dos métodos es menor que 1e-40 que es la precisión que también se consigue con cada método por separado. [fin actualización]

He calculado la serie formal de Laurent para la función de Abel de tipo Écalle utilizando Pari/GP para 509 coeficientes en números racionales exactos, lo que significa que los coeficientes de $z^{-2}$ a $z^{506}$ .

El último coeficiente tiene numerador con 1423 dígitos y el denominador con 1247 dígitos que suman un valor absoluto de aproximadamente 175 dígitos, aproximadamente -2.66945040282 E175 , por lo que la serie tiene del mismo modo que la serie comparable para la iteración fraccionaria de $\exp(x)-1$ radio de convergencia cero y cuando trazamos la curva que muestra el número de dígitos de los coeficientes no nulos por $\log_{10}(|a_k|)$ obtenemos la forma típica del hockeystick vertical.

Aquí están los principales 11 términos no nulos de la serie de Laurent para la función de Abel (que aquí llamo "función incompleta de Abel" hasta ahora porque el "completo" La función Abel necesita también el término del logaritmo y el término de la altura de iteración h (este es el índice n en la respuesta de Will):

  Laurent-series in z:                             3  *z^-2
                                           + 79/1050  *z^2
                                           + 29/2625  *z^4
                                    + 91543/36382500  *z^6
                              + 18222899/28378350000  *z^8
                             + 88627739/573024375000  *z^10
                        + 3899439883/142468185234375  *z^12
              - 32544553328689/116721334798818750000  *z^14
                 - 4104258052789/1554729734250000000  *z^16
 - 119345896838809094501/141343700374629565312500000  *z^18
      + 745223773943295679/3505548124370772949218750  *z^20
             + O(z^22)

$\qquad \qquad$ (Observación: véase un resumen que caracteriza la tasa de crecimiento al final (§2))

Esto da la "función incompleta de Abel" en términos de su coeffs hasta un cierto truncamiento n (todas las fórmulas en notación Pari/GP):

 abel_inc(z,n=64) = sum(k=1,n, coeff[k]*z^(k-3) )     

La función completa de Abel es entonces :

 {abel(z,h=32,n=64) = local(z_h,a); \\ give some sufficient default values 
                                    \\ in h and n for the required numerical 
                                    \\ precision of the approximate results
        z_h = sin_iter(z,h);        \\ sin_iter prev. defined as iterable sin()
      a = abel_inc(z_h,n) + 6/5*log(z_h) - h ;
      return(a); }                      

La función abel inversa debe ser implementada por algún solucionador de raíces. En Pari/GP utilicé lo siguiente, donde la función inversa de Abel se incluye en el cuerpo de la fraccionable completa sin_h() función:

    {sin_h (h = 0,z_0=1) = local(a_0,z_h,a_h);  \\ restriction abs(h)<1 
               a_0 = abel(z_0, 32, 64);         \\ get the Abel-value for z_0
                                                \\ with meaningful precision
               a_h = a_0 + h ;                  \\ comp Abel-value for z_h
                                 \\ the following is the implementation of
                                 \\ the inverse Abel-function:
         z_h = solve(z = sin(z_0),z_0, abel(z,32,64) - a_h);
         return(z_h); }

Se hace lo siguiente para aplicar lo anterior a algún ejemplo, reproduciendo la aditividad de las alturas de iteración 0.5 y 0.5 a la altura integral 1 con más de 40 precisión de los dígitos:

                               \\  Pari-output
   z_0  = 1                    \\  %529 = 1
   z_05 = sin_h(0.5,z_0 )      \\  %530 = 0.908708429743
   z_1  = sin_h(0.5,z_05)      \\  %531 = 0.841470984808
   z_1 - sin(z_0)              \\  %532 = -6.38920219348 E-42

A continuación muestro la lista de cálculos recalculados en la respuesta de Will con 40 dígitos correctos:

   step   z0=Pi/2 - step   abel(z0)   z05=sin_h(0.5,z0) z1=sin_h(0.5,z05)  z1 - sin(z0)
   0.00   1.57079632679  2.08962271973   1.14017947617  1.00000000000   -2.89445031739E-41
   0.05   1.52079632679  2.09536408453   1.13806963935  0.998750260395  -2.86591796888E-41
   0.10   1.47079632679  2.11273622895   1.13178674818  0.995004165278  -2.78164697945E-41
   0.15   1.42079632679  2.14218948912   1.12146458427  0.988771077936  -2.64553725829E-41
   0.20   1.37079632679  2.18449553252   1.10730765183  0.980066577841  -2.46383393292E-41
   0.25   1.32079632679  2.24078077607   1.08956885996  0.968912421711  -2.24476553049E-41
   0.30   1.27079632679  2.31257688904   1.06852649593  0.955336489126  -1.99807394218E-41
   0.35   1.22079632679  2.40189260763   1.04446448663  0.939372712847  -1.73446474837E-41
   0.40   1.17079632679  2.51131312355   1.01765794736  0.921060994003  -1.46500647333E-41
   0.45   1.12079632679  2.64413616528  0.988364216777  0.900447102353  -1.20050550750E-41
   0.50   1.07079632679  2.80455803137  0.956818478819  0.877582561890  -9.50882282773E-42
   0.55   1.02079632679  2.99792899241  0.923232674366  0.852524522060  -7.24576289372E-42
   0.60  0.970796326795  3.23110684637  0.887796468526  0.825335614910  -5.28014362671E-42
   0.65  0.920796326795  3.51295197372  0.850679308887  0.796083798549  -3.65188391373E-42
   0.70  0.870796326795  3.85503037983  0.812032915560  0.764842187284  -2.37402622132E-42
   0.75  0.820796326795  4.27262886030  0.771993802047  0.731688868874  -1.43260703471E-42
   0.80  0.770796326795  4.78624925852  0.730685613103  0.696706709347  -7.89576195851E-43
   0.85  0.720796326795  5.42385666222  0.688221187210  0.659983145885  -3.89074331205E-43
   0.90  0.670796326795  6.22434753781  0.644704322722  0.621609968271  -1.66626510284E-43
   0.95  0.620796326795  7.24305478745  0.600231264287  0.581683089464  -5.96979699941E-44
   1.00  0.570796326795  8.56077779381  0.554891942675  0.540302305868  -1.69831000319E-44

Una imagen de $y=\sin(x)$ el medio iterado $y=\sin^{\circ 0.5}(x)$ , $y=\sin^{\circ 1/3}(x)$ y $y=x$ :

Observación: en x , donde sin(x)=0 el cálculo de la función Abel se ejecuta en singularidades y el valor de la función se supone (interpolado desde su vecindad) cero.

---

En la función de Abel abel(z,h=32,n=64)=... existe el parámetro h que permite controlar la calidad de la aproximación. El sistema formal exacto la solución se da como límite cuando h va al infinito, pero aquí sólo utilizamos aproximaciones finitas. La clave es que h controla la iteración implícita del argumento z hacia el punto fijo cero, por lo que la evaluación numérica del (truncado a n coeficientes) La serie Laurent da una mejor aproximación al valor real - ¡aunque en realidad el radio de convergencia sigue siendo cero! El propósito de esas iteraciones que cambian z_h hacia cero es desplazar la posición, desde donde la serie de Laurent con el argumento z_h comienza a divergir, a índices más altos y, por tanto, a obtener más precisión. Una combinación de h=32 y n=64 para los argumentos $|z| \le 1$ es aparentemente suficiente para 40 dígitos correctos. (véase la observación (§1))

Por último, para mostrar el efecto de la h=32 iteración en el trabajo proporciono a continuación las sumas parciales de las series de Laurent para z=1 en comparación con h=4 .

En el primer ejemplo utilizo h=4 y en el segundo ejemplo utilizo h=32 .
En la tabla k es el índice del coeficiente hasta donde se calculan las sumas parciales. ps_k indica la suma parcial utilizando z_h que es el h 'th iterate from z_0=1 . Pero por conveniencia el término para el logaritmo y el h -se incluyen siempre para que podamos comparar la suma hasta este término con el valor exacto a_1 para la función Abel en z_1 :

   k    ps_k              error:  a_1 - ps_k    iteration height h=4
   0  3.05810608515        -0.0315166345810
   2  3.05810608515        -0.0315166345810
   4  3.08773833843       -0.00188438129901
   6  3.08945198975      -0.000170729978211
   8  3.08960570369     -0.0000170160371392
  10  3.08962115403    -0.00000156570332243
  12  3.08962261968   -0.000000100050871450
  14  3.08962272183  0.00000000210083986271
  16  3.08962272142  0.00000000169099804938
  18  3.08962271989       1.62746538183E-10
  20  3.08962271970      -2.97721970306E-11
  ... 
  50  3.08962271973      -3.98604755990E-18
  52  3.08962271973       7.74229820435E-19
  54  3.08962271973       1.21098784690E-18
  56  3.08962271973      -6.22150631919E-20
  58  3.08962271973      -3.98357488277E-19
  60  3.08962271973      -3.38541477910E-20
  62  3.08962271973       1.42850133024E-19

Vemos que con la altura de iteración h=4 llegamos a un error absoluto menor que 1e-18 en el 64. término. Y a continuación, la altura de iteración h=32 proporciona una precisión con un error absoluto inferior a 1e-40 con eso 64 términos utilizados:

   k    ps_k              error:  a_1 - ps_k    iteration height h=32
   0  3.08337725463         -0.00624546510435
   2  3.08337725463         -0.00624546510435
   4  3.08954701281       -0.0000757069234782
   6  3.08962130264      -0.00000141708899538
   8  3.08962269011    -0.0000000296188288642
  10  3.08962271915  -0.000000000581829933894
  12  3.08962271972        -8.31025344698E-12
  ... ...                  ...
  52  3.08962271973        -3.06907747463E-37
  54  3.08962271973         5.27409063179E-37
  56  3.08962271973         2.10119895640E-38
  58  3.08962271973        -6.82487772781E-39
  60  3.08962271973        -5.39925105785E-40
  62  3.08962271973         9.44571568505E-41

---

(§1): una suma de Noerlund, como la tengo propuesta en algunos tratados para la evaluación de los iterados fraccionarios de la $\exp(x)-1$ podría dar aproximaciones arbitrarias también, pero ahora me parece que tal procedimiento de suma era a lo sumo necesario aquí por razones teóricas para demostrar la sumabilidad de la serie de Laurent para la función de Abel.

(§2): Un breve resumen sobre los primeros 512 coeficientes de la serie abel_inc():

 index      value             index   value            index    value               index    value
    0           3.000000000   47                     0   92       -0.005185699555  496       4.633504372E168 
    1                     0   48  -0.00000003870320993   93                     0  497                     0
    2                     0   49                     0   94         0.01347223160  498      -4.983759375E169
    3                     0   50  0.000000006386371562   95                     0  499                     0
    4         0.07523809524   51                     0   96         0.03559427183  500      -8.187596780E170
    5                     0   52   0.00000006229599636   97                     0  501                     0
    6         0.01104761905   53                     0   98        -0.06747379661  502       8.333103850E171
    7                     0   54   0.00000001451248843   99                     0  503                     0
    8        0.002516127259   55                     0  100         -0.2528544049  504       1.467790435E173
    9                     0   56   -0.0000001074166810  101                     0  505                     0
   10       0.0006421408926   57                     0  102          0.3439480705  506      -1.412786474E174
   11                     0   58  -0.00000007200630916  103                     0  507                     0
   12       0.0001546666126   59                     0  104           1.879638019  508      -2.669450403E175
   13                     0   60    0.0000001982539503  105                     0  ...          ...
   14      0.00002737060121   61                     0  106          -1.706858981
   15                     0   62    0.0000002440284845  107                     0
   16   -0.0000002788226624   63                     0  108          -14.69827943
   17                     0   64   -0.0000003845753696  109                     0
   18    -0.000002639853064   65                     0  110           7.295584305
   19                     0   66   -0.0000007917263057  ...             ...         
   20   -0.0000008443665796   67                     0
   ...                  ...   ...                  ...

Answer 2

Si se crea la matriz de Bell para la función $f(x) = \sin(x)$ , digamos que SI entonces se puede calcular el logaritmo matricial de SI , digamos que SIL \= MLog(SI) . A continuación, un poder formal de SI es SIP(h) \= MExp(h*SIL) y la matriz de Bell para la función dependiente de la altura $\operatorname{sin\_iter}(x,h)$ que tiene polinomios en h para los coeficientes en x . SI comienza con

  1        .     .       .
  0        1     .       .
  0        0     1       .
  0     -1/6     0       1
  0        0  -1/3       0
  0    1/120     0    -1/2
  0        0  2/45       0
  0  -1/5040     0  13/120

donde la columna 1 contiene los coeficientes de la serie de potencias $\sin(x)$ , columna 2 que para $(\sin(x))^2$ , columna 0 que para $(\sin(x))^0 = 1$ y de forma similar para todas las demás columnas k .

La matriz-logaritmo SIL comienza con

  0         .      .      .     .     .  .  . 
  0         0      .      .     .     .  .  . 
  0         0      0      .     .     .  .  . 
  0      -1/6      0      0     .     .  .  . 
  0         0   -1/3      0     0     .  .  . 
  0     -1/30      0   -1/2     0     0  .  . 
  0         0  -1/15      0  -2/3     0  0  . 
  0  -41/3780      0  -1/10     0  -5/6  0  0 

Aquí la columna k es el k 'th múltiplo de la columna 1 desplazado $k-1$ fila hacia abajo.

A continuación, la columna 1 de SIP(h) \= MExp(h*SIL) es

                              0
                              1
                              0
                         -1/6\*h
                              0
                1/24\*h^2-1/30\*h
                              0
  -5/432\*h^3+1/45\*h^2-41/3780\*h

y la función $\operatorname{sin\_iter}$ es

$$\operatorname{sin\_iter}(x,h) = 1 x - h \cdot {x^3 \over 3!} + (5 h^2-4 h) \cdot {x^5 \over 5!} - (...) \cdot {x^7 \over 7!} + O(x^9)$$

Inserción de $h={1 \over 2}$ te da la serie de potencias para la media iteración.

Utilizando 64 parece que el radio de convergencia para $h=\frac 12$ será 1 ya que los valores absolutos de los coeficientes parecen estabilizarse en el intervalo $\pm 1E-7$ pero voy a ver esto más adelante en el día.

[Actualización]

utilizando 256 se produce un claro crecimiento de los coeficientes. Mirando el logaritmo de los valores absolutos de esos coeficientes tenemos una impresión aproximada. Véase aquí:

http://go.helms-net.de/math/images/sincoeff_c.png

Estos son los coeficientes en $x^{123},x^{125},x^{127}$ y $x^{251}, x^{253}, x^{255}$ :

c_123     -2156.72733764089915  // 4 digits
c_125     31313.42875545542423  // 5 digits
c_127     34859.64557727596911  // 5 digits 
...    
c_251       -35365220492708296140377087748804440170254492009.570  // 46 digits    
c_253     -1378449672866233726070664896135098313484573633108.4    // 48 digits    
c_255       987848122496441964413343332623221752473112662017.00   // 47 digits    

Las diferencias de los logaritmos son también cocientes de los coeficientes Por el gráfico de las diferencias obtenemos también una tendencia de aumento logarítmico. (Si las diferencias siguen aumentando entonces el radio de convergencia de la serie de potencias es cero, porque la tasa de crecimiento de los valores absolutos de los coeficientes es hipergeométrica)

http://go.helms-net.de/math/images/sincoeff_d.png

[fin de la actualización]

Pari/GP calcula esto bastante rápido, tardó, digamos 5 segundos en manejar matrices de 64 términos.

[actualización2 , febrero de 2016]
Un método muy sencillo para obtener la serie de potencias formal para el medio iterado de la función seno es combinar la función de expansión de taylor interna de Pari/GP, el serreverse() con el algoritmo Newton para la raíz cuadrada. Para un escalar $t$ como raíz de un determinado $z$ es $t_{k+1}=(z/t_k+t_k)/2$ y aquí interpretamos $t$ y $z$ como serie de potencias, donde $t$ también es invertible.
Este es el protocolo de la sesión Pari/GP:

t=x + O(x^12)   \\ Initialization of the Newton-algorithm with a simple power series
 %76 = x + O(x^12)  \\ the protocol that Pari/GP shows in the dialog

t = (sin(serreverse(t))+t)/2   \\ first iteration
 %77 = x - 1/12*x^3 + 1/240*x^5 - 1/10080*x^7 + 1/725760*x^9 - 1/79833600*x^11 + O(x^12)

t = (sin(serreverse(t))+t)/2   \\ secons iteration
 %78 = x - 1/12*x^3 - 1/160*x^5 - 11/5040*x^7 - 11/17920*x^9 - 2425/12773376*x^11 + O(x^12)

t = (sin(serreverse(t))+t)/2
 %79 = x - 1/12*x^3 - 1/160*x^5 - 53/40320*x^7 - 341/1935360*x^9 + 44311/638668800*x^11 + O(x^12)

t = (sin(serreverse(t))+t)/2
 %80 = x - 1/12*x^3 - 1/160*x^5 - 53/40320*x^7 - 23/71680*x^9 - 138913/1277337600*x^11 + O(x^12)

t = (sin(serreverse(t))+t)/2
 %81 = x - 1/12*x^3 - 1/160*x^5 - 53/40320*x^7 - 23/71680*x^9 - 92713/1277337600*x^11 + O(x^12)

t = (sin(serreverse(t))+t)/2  
 %82 = x - 1/12*x^3 - 1/160*x^5 - 53/40320*x^7 - 23/71680*x^9 - 92713/1277337600*x^11 + O(x^12) // the solution becomes stable for the first coefficients

Los coeficientes pueden ampliarse de forma muy sencilla a un índice arbitrario, sólo hay que ajustar la expansión de la serie de potencias por defecto a la precisión deseada y definir la inicialización de t en consecuencia.

Answer 3

(Esto debería ir como comentario, pero ha sido imposible). @Aaaron: He subido una lista de los primeros 128 coeficientes no nulos, ver:

http://go.helms-net.de/math/tables/sinxcoeffs.htm

También hay una rutina para Pari/GP para calcular el sqrt de una matriz de Bell triangular inferior (la matriz SI en mi respuesta anterior) Con esto se puede calcular la serie de potencias para la media iteración (por la columna 1 del sqrt de SI) en un segundo, incluso si el tamaño de la matriz es 256x256.

\\\\ square-root of a lower triangular Bell-matrix
\\\\ only implemented for operator/"Bell"-matrices for functions
\\\\ where f(x) = ax  + bx^2+ cx^3 + ... with a>0
\\\\
 trisqrt(m) = local(tmp, rs=rows(m), cs=cols(m), c);
  tmp=matrix(rs,cs,r,c,if(r==c,sqrt(m\[r,r\])));
  for(d=1,rs-1,
       for(r=d+1,rs,
          c=r-d;
          tmp\[r,c\]=(m\[r,c\]-sum(k=c+1,r-1, tmp\[r,k\]\*tmp\[k,c\]) )/(tmp\[c,c\]+tmp\[r,r\])
          );
      );
 return(tmp);