Cauchy Integral Experiment

[mike3]

Hello mike3. I looked at your code.
I would mention two the defects, that may be important:

1. The code does not provide any graphical output.
It is difficult to localize the step responsible for the instability you report.

No, but Pari/GP's "plot" and "ploth" can be used to draw 1-dimensional graphs.

2. The routine CauchyIntegralConst looks suspicious.
I doubt if it treats well the cutlines of the logarithmic function.
I suggest that you plot the result of the evaluation with CauchyIntegralConst
as function of position of the point of evaluation in the complex plane,
the cut lines are supposed to be seen at such a plot;
and let us compare it to the figure 3 from my paper in Mathematics of Computation.

CauchyIntegralConst just gives the integral of a constant function along a segment. In this case it's used along a horizontal segment that runs just above/below (gets closer and closer to as the number of nodes is increased, though) the singularity of the Cauchy integrand. Such a segment does not cross the cutline of the log (at least if you use a branch whose cutline is along the real axis, in this case I use the principal branch), it is always parallel to it.

Can this also be used for the imaginary-periodic bases \( 1 < b < e^{1/e} \)?

I doubt if you can generalize the algorithm to other bases before you make it work well for base b=e.
Also, there should be additional problems with the closure of the contour of integration. They are discussed at
http://math.eretrandre.org/tetrationforu...hp?tid=248

Yeah, I suppose I should get it working for base e first

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?

That is because you try to debugg your code in "one piece", without to arrange the routines for the internal tests of your intermediate steps and service functions.
Consider to upload the plot of your approximation of tetration along the imaginary axis, just before and just after the "wild divergence" you report;
then, I hope, it will be possible to reveal the cause of such a behavior.

Here's the initial approximation, generated by

Code:

v = MakeApprox(209, 10.0);

giving 209 nodes at A = 10.0. I then graph with

Code:

ploth(n=1,209,[real(v[floor(n+0.5)]),imag(v[floor(n+0.5)])]);

(or psploth to dump to postscript file.)

This yields:
   

After one iteration, done with

Code:

v = OneIteration(v, 10.0);

it looks like:
   

After another iteration:
   

This is done with the half-circle end contour, which means you must change the OneIteration routine to:

Code:

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=1,N,
                zi = -I*A + (i-1)*dz;
                I1 = CauchyIntegrateArr(1-I*A, 1+I*A, expv, zi);
                I2 = L/2; /* see this? */
                I3 = -CauchyIntegrateArr(-1-I*A, -1+I*A, logv, zi);
                I4 = conj(L)/2; /* and this? */
                newv[i] = I1+I2+I3+I4);
            return(newv);
}

It grows up faster the more we do.

If I use the normalized one (OneIterNrml) it doesn't seem to "wildly diverge" now that I look at it (I'm not sure what I was looking at, maybe I had a tweak in there that I don't know), but it still does not give the right function -- here's a graph after I iterate enough normalized iterations for it to settle (number of nodes is 209, A=10.0):
   

You can see the wrongness of it.

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.

Could you plot this residual as function of the imaginary part of the argument?
How does this residual depend on the order of update of the values approximation tetration in your array?

Well I ran (using the usual iteration routine not the one with the half-circle business as that doesn't work yet, see above)

Code:

v = MakeApprox(209, 10.0);
for(i=1,32,v = OneIterationNrml(v, 10.0));
for(i=1,6,v = OneIteration(v, 10.0));
bias = FindBias(v, 10.0); /* find offset so function evaluated at 0 yields 1 exactly */
psploth(X=-9.0,9.0,abs(Residual(v, 10.0, bias + I*X)));

Residual at zero is ~ 0.00002033 in magnitude. Value of tetration at z = 0.5 is approximately 1.646456 (looks good to almost 4 places past the point ("right" value is 1.646354).).

The graph generated by the above shows the magnitude of the residual from z = -9.0i to +9.0i:
   

I did not run the Cauchi with the Simpson approximation of integrals.
I used the approcimation with "rectangles" (second order of approximation)
and, after to get the residual of order of 10^(-5), I switched to the Gauss-Legendre (infinite order of approximation).

I'm using the Simpson because it has more accuracy than the rectangular (so I don't need so many nodes, and that makes it faster).

In my code, I used to update first the only even nodes: 0,2,4,6,..
(The odd points were updated authomatically due to the symmetry);
after to see that it provides the residual at the level of the rounding errors,
I stoped to play with the order of updates;
I did not test many other ways to choose the order of uprates.
However, I expect, that the final result does not depend on the way we approximate the integrals. If the procedure converges, it converges to the holomorphic function, satisfying the equation tor the superexponential.

I'd think so. The normalized iteration doesn't _seem_ to need a special order of evaluation, however.

How does the residual depend on the step of integration?
How does it depend on the parameter A ?

Changing "A" from 10 to 15 (while holding node number at 209 and 32 normalized+6 unnormalized iterations) causes residual at 0 to leap to 0.1415. Ow.

And what's "step"? Number of nodes?