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 \( 1 < b < e^{1/e} \)? 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 \( 0 < b < e^{-e} \), 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 \( \frac{f(z)}{z - a} \)) 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 \( 1 < b < e^{1/e} \)? 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 \( 0 < b < e^{-e} \), 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 \( \frac{f(z)}{z - a} \)) 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?