Cauchy Integral Experiment

[mike3]

Aagh... it looks like I've noticed additional problems now... geez...

For A = 10 and 209 nodes, and 32 normalized iterations followed by 6 unnormalized ones it does work good, giving residual at 0 of magnitude ~4.68 x 10^-7. Yet if I bump the node count to 801 and A to 24 it gets wrong, residual plummets horrificially to magnitude 0.01 at 32 normalized iterations. 6 unnormalized ones just makes it worse, 0.02.

The graphs attached to this show what happens after 0, 8, 16, 24, and 32 iterations (normalized) respectively. It's interesting. It's like "waves" are being emitted from the center and rippling on out toward the ends. No, seriously, flip between the pics. You can see the "waves" actually moving. So now I have a nice wave tank. Which I suppose is useful as a neat gimmick, but not as a thing for calculating Tetration

Plotting code (v = node array):

Code:

psploth(n=1,801,[real(v[floor(n+0.5)]), imag(v[floor(n+0.5)])])

What is wrong?

Code for the newest program:

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 at point zaround upper/lower semicircle centered at +/-iA. */
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));

/* 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=1,N,
                zi = -I*A + (i-1)*dz;
                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);
                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 = UpperSemicircle(L, A, a);
            I3 = -CauchyIntegrateArr(-1-I*A, -1+I*A, logv, a);
            I4 = LowerSemicircle(conj(L), 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));
}