09/13/2009, 10:55 AM
Hi.
On the thread
http://math.eretrandre.org/tetrationforu...273&page=3
the following integral+continuum product iterative formula was mentioned for tetration to base e:
\( f(x)=f'(0)\int_{-1}^x \prod _x f(x+1)dx \)
So I put together a Pari/GP code to road-test it numerically, to see if maybe it would converge on something. I use the formula
\( \prod_x f(x) = \exp\left(\sum_x \log(f(x))\right) \)
to convert the product to a sum (which also has the nice effect of changing the x+1 to x), and then the Faulhaber's formula (expressed using Bernoulli polynomials/coefficients) is used to do the continuous sum on a set of Taylor coefficients. The Bell polynomials can be used to compute the exponential exp of the series.
When I run the code, I notice something strange. For a small number of terms in the test array (I think around 15), it seems to converge and give like 4-5 digits of accuracy. For a larger number (say 28 or more), however, it starts converging but then it blows up.... Not sure if this is a problem with the code, or what... That's what I want to find out. As ultimately I'd like to be able to explore the general formula for other bases, including the really intriguing case of \( 0 < b < e^{-e} \)...
To use the code, run tetiter several times over a vector containing coefficients for an initial guess (all zeroes will work), and then use tsum to evaluate it at a given x-value. f0 should be set equal to a numerical value for the derivative at 0 (see the original thread for the value I used).
On the thread
http://math.eretrandre.org/tetrationforu...273&page=3
the following integral+continuum product iterative formula was mentioned for tetration to base e:
\( f(x)=f'(0)\int_{-1}^x \prod _x f(x+1)dx \)
So I put together a Pari/GP code to road-test it numerically, to see if maybe it would converge on something. I use the formula
\( \prod_x f(x) = \exp\left(\sum_x \log(f(x))\right) \)
to convert the product to a sum (which also has the nice effect of changing the x+1 to x), and then the Faulhaber's formula (expressed using Bernoulli polynomials/coefficients) is used to do the continuous sum on a set of Taylor coefficients. The Bell polynomials can be used to compute the exponential exp of the series.
When I run the code, I notice something strange. For a small number of terms in the test array (I think around 15), it seems to converge and give like 4-5 digits of accuracy. For a larger number (say 28 or more), however, it starts converging but then it blows up.... Not sure if this is a problem with the code, or what... That's what I want to find out. As ultimately I'd like to be able to explore the general formula for other bases, including the really intriguing case of \( 0 < b < e^{-e} \)...
Code:
/* Continuous sum transformation. */
bernoulli(n) = if(n==0, 1, (-1)^(n+1) * n * zeta(1 - n));
contsumcoef(series, k) = if(k==0, 0, k!*sum(n=1,matsize(series)[2], series[n]/(n!)*binomial(n, k)*bernoulli(n-k)));
sumoperator(series) = { local(v=series); for(i=1,matsize(series)[2],v[i] = contsumcoef(series, i-1)); return(v); }
tsum(series, x) = sum(k=0,matsize(series)[2]-1,series[k+1]/k! * x^k);
/* Compute exponential of a Taylor series. */
bellcomplete(coefs) = {
local(n=matsize(coefs)[2]);
local(M=matrix(n,n));
for(i=1,n,\
for(j=1,n,\
M[i, j] = 0;
if(j-i+1 > 0, M[i, j] = binomial(n-i, j-i)*coefs[j-i+1]);
if(j-i+1 == 0, M[i, j] = -1);
);
);
return(matdet(M));
}
expoperator(series) = {
local(v = series);
for(i=1,matsize(series)[2]-1,
v[i+1] = bellcomplete(vector(i,n,series[n+1])));
v[1] = 1;
return(v*exp(series[1]));
}
/* Compute integral of a Taylor series at lower bound a. */
intoperator(series, a) = {
local(v = series);
for(i=1,matsize(series)[2]-1,
v[i+1] = series[i]);
v[1] = -sum(i=2,matsize(series)[2],v[i]*(a^(i-1))/((i-1)!));
return(v);
}
/* Tetration iteration */
tetiter(f0, series) = f0*intoperator(expoperator(sumoperator(series)), -1);
To use the code, run tetiter several times over a vector containing coefficients for an initial guess (all zeroes will work), and then use tsum to evaluate it at a given x-value. f0 should be set equal to a numerical value for the derivative at 0 (see the original thread for the value I used).