09/24/2009, 10:07 PM

Hi.

Here's the newest code I put together for running the Cauchy Integral method of Kouznetsov. It's written for the Pari/GP arbitrary precision calculation package, and set up for tetration of base e.

I use Simpson's rule to do the numerical integration. To cure the "drift" which was causing trouble in earlier attempts (not posted), I went and used a technique of "normalizing" the array of imaginary-axis values so that the center node (representing the value of the function at z = 0) is 1 after every iteration.

With 209 nodes, A = 10, 32 normalized iterations followed by 6 unnormalized ones (causes a little "drift"/"sag" of the middle peak but seems to improve accuracy better once this is corrected for via looking for the necessary offset/bias by which to shift the function reconstructed from these imag-axis values through the Cauchy formula so it has value 1 at z = 0) gives a real-axis residual of ~0.00002.

My questions on this are: Can this also be used for the imaginary-periodic bases ? I heard here something about incorporating the values along the real axis and its periodic offset (or presumably any 2 axes separated by 1 period), but what contour does one use to update those arrays?

And if one can do that, could one also use something similar with storing values along the real axis for , such as b = 0.04? As I suspect that such a base should decay to the repelling fixed point (for b = 0.04, ~0.33747) as x -> +ioo, but the behavior is unknown at -ioo. I infer this because the quasiperiod T for b = 0.04 is ~1.9986 + 0.0526i, the line of which (i.e. that parameterized by multiplying T by a real variable t) sloping gently from the lower left to upper right, and toward the left half-plane we should expect the function to (except for the singularities, at z = -2, -3, -4, -5, etc.) decay to a fixed point, and this QP line seems to be a "border" (not a hard and fast one of course as that would be a discontinuity!) of sorts (at least looking at the graphs of bases b = e, b = 2) separating different types of behaviors. So presumably, I'd think it should also decay to the fixed point as z -> +ioo, as that is also in the same QP-line-half of the plane as going left puts you in.

But without a way to update the node values on the real axis, I can't make a program to try it out.

Another question is that if I try using the half-circle contour for the endpoints (as opposed to not updating them at all, which is what I do in the above), the thing diverges wildly. Presumably, the values of the integral around the endpoints should approximate L/2 and conj(L)/2 as the Cauchy integral of a constant function around a half-circle whose center is on the parameter/evaluation point ("a" in the integrand ) is K/2, where K is the value of the constant function, and tet is assumed approximately constant at large z = ix, x real. Why is that?

Here's the newest code I put together for running the Cauchy Integral method of Kouznetsov. It's written for the Pari/GP arbitrary precision calculation package, and set up for tetration of base e.

I use Simpson's rule to do the numerical integration. To cure the "drift" which was causing trouble in earlier attempts (not posted), I went and used a technique of "normalizing" the array of imaginary-axis values so that the center node (representing the value of the function at z = 0) is 1 after every iteration.

With 209 nodes, A = 10, 32 normalized iterations followed by 6 unnormalized ones (causes a little "drift"/"sag" of the middle peak but seems to improve accuracy better once this is corrected for via looking for the necessary offset/bias by which to shift the function reconstructed from these imag-axis values through the Cauchy formula so it has value 1 at z = 0) gives a real-axis residual of ~0.00002.

Code:

`\p 16;`

L = 0.3181315052047641353126542516 + 1.337235701430689408901162143*I; /* fix of log */

/* Compute delta spacing for N equidistant nodes from z1 to z2, so that z1 + (N-1)delta = z2 */

GetDelta(N, z1, z2) = (z2 - z1)/(N-1);

/* Fill up array with initial approximation of tet from -iA to +iA */

tetapprox(ix) = {

local(t=0);

if(ix < -3.0, return(conj(L)));

if(ix > 3.0, return(L));

return(subst(polinterpolate([-3.0, 0, 3.0], [conj(L), 1, L]), x, ix));

}

MakeApprox(N, A) = {

local(v=vector(N));

local(dx=GetDelta(N,-A,A));

local(ix=-A);

for(i=1,N,

ix = -A + ((i-1)*dx);\

v[i] = tetapprox(ix);\

);

return(v);

}

/* Simpson's rule numerical integration of function given at N equidistant nodes from z1 to z2 in vector v */

/* Odd number of nodes only! */

Twiddle(n, N) = if(n == 1 || n == N, 1, if(n%2 == 1, 2, 4));

SimpIntegrateArr(z1, z2, v) = {

local(N=matsize(v)[2]);

local(dz=GetDelta(N, z1, z2));

return((1/3)*dz*sum(k=1,N,Twiddle(k, N)*v[k]));

}

/* Do Cauchy integral using values of function along z1...z2 in v for point a */

CauchyIntegrateArr(z1, z2, v, a) = {

local(N=matsize(v)[2]);

local(dz=GetDelta(N, z1, z2));

local(vp=v);

local(zi=z1);

for(i=1,N,\

zi = z1+((i-1)*dz);\

vp[i] = vp[i]/(zi - a););

return(1/(2*Pi*I)*SimpIntegrateArr(z1, z2, vp));

}

/* Do "Cauchy integral" of a constant function K from z0...z1 with parameter a. */

CauchyIntegralConst(K, z0, z1, a) = (K/(2*Pi*I))*(log(z1 - a) - log(z0 - a));

/* Do an iteration over a vector v representing values of tetration from -iA to +iA */

OneIteration(v, A) = {

local(N=matsize(v)[2]);

local(newv = v);

local(expv = exp(v));

local(logv = log(v));

local(zi=-I*A);

local(dz=GetDelta(N,-I*A,I*A));

local(I1=0);

local(I2=0);

local(I3=0);

local(I4=0);

for(i=2,N-1,

zi = -I*A + (i-1)*dz;

I1 = CauchyIntegrateArr(1-I*A, 1+I*A, expv, zi);

I2 = CauchyIntegralConst(L, 1+I*A, -1+I*A, zi);

I3 = -CauchyIntegrateArr(-1-I*A, -1+I*A, logv, zi);

I4 = CauchyIntegralConst(conj(L), -1-I*A, 1-I*A, zi);

newv[i] = I1+I2+I3+I4);

return(newv);

}

/* Do a normalized iteration (requires odd number of nodes) */

OneIterationNrml(v, A) = {

local(N=matsize(v)[2]);

local(newv = OneIteration(v, A));

return(newv/newv[(N-1)/2+1]);

}

/* Recover value at a point from Cauchy integral and values on imag axis */

PointEvaluate(v, A, a) = {

local(N=matsize(v)[2]);

local(expv = exp(v));

local(logv = log(v));

local(I1=0);

local(I2=0);

local(I3=0);

local(I4=0);

I1 = CauchyIntegrateArr(1-I*A, 1+I*A, expv, a);

I2 = CauchyIntegralConst(L, 1+I*A, -1+I*A, a);

I3 = -CauchyIntegrateArr(-1-I*A, -1+I*A, logv, a);

I4 = CauchyIntegralConst(conj(L), -1-I*A, 1-I*A, a);

return(I1+I2+I3+I4);

}

/* Find where function approximated by v has value 1 */

FindBias(v, A) = {

local(r=0);

local(diff=0);

for(i=1,32,

diff = (PointEvaluate(v, A, r + 1e-12) - PointEvaluate(v, A, r))/(1e-12);\

r=r-((PointEvaluate(v, A, r) - 1)/diff));

return(r);

}

/* Compute residual of an approximation given by imag axis values & Cauchy integral */

Residual(v, A, bias) = log(PointEvaluate(v, A, 0.5 + bias)) - PointEvaluate(v, A, -0.5 + bias);

/* Do general tetration */

Tetrate(v, A, bias, z) = {

if(real(z) < -0.5, return(log(Tetrate(v, A, bias, z+1))));

if(real(z) > 0.5, return(exp(Tetrate(v, A, bias, z-1))));

return(PointEvaluate(v, A, z + bias));

}

My questions on this are: Can this also be used for the imaginary-periodic bases ? I heard here something about incorporating the values along the real axis and its periodic offset (or presumably any 2 axes separated by 1 period), but what contour does one use to update those arrays?

And if one can do that, could one also use something similar with storing values along the real axis for , such as b = 0.04? As I suspect that such a base should decay to the repelling fixed point (for b = 0.04, ~0.33747) as x -> +ioo, but the behavior is unknown at -ioo. I infer this because the quasiperiod T for b = 0.04 is ~1.9986 + 0.0526i, the line of which (i.e. that parameterized by multiplying T by a real variable t) sloping gently from the lower left to upper right, and toward the left half-plane we should expect the function to (except for the singularities, at z = -2, -3, -4, -5, etc.) decay to a fixed point, and this QP line seems to be a "border" (not a hard and fast one of course as that would be a discontinuity!) of sorts (at least looking at the graphs of bases b = e, b = 2) separating different types of behaviors. So presumably, I'd think it should also decay to the fixed point as z -> +ioo, as that is also in the same QP-line-half of the plane as going left puts you in.

But without a way to update the node values on the real axis, I can't make a program to try it out.

Another question is that if I try using the half-circle contour for the endpoints (as opposed to not updating them at all, which is what I do in the above), the thing diverges wildly. Presumably, the values of the integral around the endpoints should approximate L/2 and conj(L)/2 as the Cauchy integral of a constant function around a half-circle whose center is on the parameter/evaluation point ("a" in the integrand ) is K/2, where K is the value of the constant function, and tet is assumed approximately constant at large z = ix, x real. Why is that?