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}{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
And if one can do that, could one also use something similar with storing values along the real axis for
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
