Cauchy Integral Experiment

[mike3]

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.

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?

[Kouznetsov]

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.

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.

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

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.

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?

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).
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.

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

[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?

[Kouznetsov]

Ok, Mik. Thank you for the pics.
They confirm my suspections about the ends of the interval: a lot of variation there.

Does your iteration keep the function which is just constant L everywhere?
Could you check this?

Does your iteration keep the function which is precise tetrational at the first?
Could you check this?

As for the residual, the single value at zero does not brung much information.
You should plot the residual along the imaginary axis.

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

I mean the distance between nodes. As I understand, your nodes are equidistant.

[mike3]

Ok, Mik. Thank you for the pics.
They confirm my suspections about the ends of the interval: a lot of variation there.

Yeah. It's the endpoints that are apparently causing all the trouble Which I find extremely peculiar, esp. the circle contour, since mathematically that _should_ work.

Does your iteration keep the function which is just constant L everywhere?
Could you check this?

OK, now that's weird. A massive disturbance develops in the initially flat line (line pair, actually, due to there being both real and imag parts) when I do that:

   

(one iteration of OneIteration(v, 10.0) on 209 nodes all set to L)

As you can see there's that whacko curve and that shouldn't be there. The pass should hardly change the array if change it at all. There's something very very bad happening near the left endpoint... I'm betting this is the source of the problem... (note that I don't update the two endpoints themselves as they lie right on the contour, remember the semicircle thing was failing... but maybe it wasn't, maybe there's something wrong with my integration routine, I'll have to check it some more)

Does your iteration keep the function which is precise tetrational at the first?
Could you check this?

I suppose I could, but given the problem I just mentioned above, I think it'd be best to address that first.

As for the residual, the single value at zero does not brung much information.
You should plot the residual along the imaginary axis.

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

I mean the distance between nodes. As I understand, your nodes are equidistant.
[/quote]

So if we hold A constant then that would translate to increasing the node number?

[mike3]

OK, that big failure was just because I forgot to change conj(L) to L in the OneIteration routine when running it over the constant-L array... it was still assuming it was tetration But when that's done it keeps the array mostly constant (more or less up to the error of the Simpson rule), so I'm still perplexed.

[mike3]

Oh dear. I finally see why the circle contour wasn't working! Geez, why didn't I think of this earlier...

The expression K/2 for the half-circle Cauchy integral of a constant function _only_ makes sense when the parameter (the "a" in the integrand) is at the center of the circle, i.e. when a = +/-iA. But that's only when you're evaluating the endpoints... when you're not, then a is something else, namely whatever point you're evaluating at. Otherwise, the integral is more complicated.

We have

\( K = \frac{1}{2 \pi i} \oint \frac{K}{z - a} dz \)

Now the cricle is parameterized by \( C(t) = e^{it} + iA \) for the top and \( C(t) = e^{it} - iA \) for the bottom. This does NOT cancel a, so we get for the semicircles

\( I_2 = \frac{L}{2 \pi i} \int_{0}^{\pi} \frac{ie^{it}}{e^{it} + iA - a} dt \)
(upper)

and

\( I_4 = \frac{\bar{L}}{2 \pi i} \int_{\pi}^{2\pi} \frac{ie^{it}}{e^{it} - iA - a} dt \)
(lower)

Now I see why you were warning about branchcuts and log and all that as this involves log now!!! Let's see how it goes. The indefinite integral of \( I_2 \) is (substitution \( u = e^{it} \) and \( du = ie^{it} dt \))

\( \int \frac{e^{it}}{e^{it} + iA - a} dt = \int \frac{du}{u + iA - a} = \log(u + iA - a) + const = \log(e^{it} + iA - a) + const. \)

So that for \( I_4 \) is

\( \int \frac{ie^{it}}{e^{it} - iA - a} = \log(e^{it} - iA - a) + const. \).

Now of course as you mention, the devil is in the branch cut. For the upper integral, there is no problem, as at the "lowest" value for a, i.e. -iA, we get \( log(e^{it}) \) whcih equals it for t in (-pi, pi] and t goes from 0 to pi. For the lower integral, however, there is a problem, as we go from pi to 2pi or alternatively -pi to 0 but -pi is not in the principal branch of logarithm. Thus when a is +iA, it breaks. The alternative then is to use a branch of log that is "continuous from below", not from above, at the negative real axis, which I'll denote as \( blog \) for "below-log". This equals the principal log(z) at all z except negative reals, where it instead equals log(z) - 2pi i. Thus we get

\( \int_{0}^{\pi} \frac{ie^{it}}{e^{it} + iA - a} dt = \log(-1 + iA - a) - \log(1 + iA - a) \)

\( \int_{\pi}^{2\pi} \frac{ie^{it}}{e^{it} - iA - a} dt = \int_{-\pi}^{0} \frac{ie^{it}}{e^{it} - iA - a} dt = \mathrm{blog}(1 - iA - a) - \mathrm{blog}(-1 - iA - a) \)

and so the sum of the integrals at the upper/lower part of the contour is

\( K(z) = \frac{L}{2\pi i} \left(\log(-1 + iA - z) - \log(1 + iA - z)\right) + \frac{\bar{L}}{2\pi i} \left(\mathrm{blog}(1 - iA - z) - \mathrm{blog}(-1 - iA - z)\right) \).

So we get the Pari/GP code

Code:

blog(z) = if(imag(z) == 0 && real(z) < 0, log(z) - 2*Pi*I, log(z));
UpperSemicircle(K, A, z) = (K/(2*Pi*I))*(log(-1 + I*A - z) - log(1 + I*A - z));
LowerSemicircle(K, A, z) = (K/(2*Pi*I))*(blog(1 - I*A - z) - blog(-1 - I*A - z));

and we change the I1...I4 in OneIteration and PointEvaluate to

Code:

(OneIteration)
I1 = CauchyIntegrateArr(1-I*A, 1+I*A, expv, zi);
I2 = UpperSemicircle(L, A, zi);
I3 = -CauchyIntegrateArr(-1-I*A, -1+I*A, logv, zi);
I4 = LowerSemicircle(conj(L), A, zi);

(PointEvaluate)
I1 = CauchyIntegrateArr(1-I*A, 1+I*A, expv, a);
I2 = UpperSemicircle(L, A, a);
I3 = -CauchyIntegrateArr(-1-I*A, -1+I*A, logv, a);
I4 = LowerSemicircle(conj(L), A, a);

With this, the problem appears solved and I now can obtain exquisite convergence. The graph on the imag. axis is finally correct:

   

Residual is extremely good -- at z = 0 (i.e. of corrected function biased so it equals 1 at z = 0), 209 nodes, A = 10 its magnitude is just a hair over 0.0000004 after 32 normalized iterations and 6 unnormalized ones. Value of tetration at z = 0.5 is 1.6463536, which when rounded to 6 places past the point yields \( ^{1/2}e \approx 1.646354 \). Even at z = 9i and z = -9i the residual still only has magnitude a little over 0.0000006.

[Kouznetsov]

OK, that big failure was just because I forgot to change conj(L) to L in the OneIteration routine when running it over the constant-L array... it was still assuming it was tetration But when that's done it keeps the array mostly constant (more or less up to the error of the Simpson rule), so I'm still perplexed.

If it keeps the constant L, good!
Could you plot please, the deviation the function at the first iteration from the constant?
Could you plot it along the imaginary axis from -iA-i to iA+i ?
We'll see, where the approximation is worse, at the center or at the end.

If it keeps at least 4 digits, you may run the second teat I suggested:
for the precise initial, plot the initial and the first iteration in the range -iA-i to iA+i .

[mike3]

If you read the latest post above, I also got the circle contour working, and got it to converge properly for the tetration now (as the endpoint problem was resolved).

[Kouznetsov]

If you read the latest post above, I also got the circle contour working, and got it to converge properly for the tetration now (as the endpoint problem was resolved).

I did not read, I was too slow typing. Sorry. Now I see.
As I understand, now you confirm that the algorithm runs well and converges, do you?
Congratulations!
Could you print your version of table 1 from
http://www.ils.uec.ac.jp/~dima/PAPERS/20...pRepri.pdf
?
Let us compare the values.