Thread Rating:
  • 0 Vote(s) - 0 Average
  • 1
  • 2
  • 3
  • 4
  • 5
Cauchy Integral Experiment
#11
(09/28/2009, 07:20 AM)mike3 Wrote: With this, the problem appears solved and I now can obtain exquisite convergence.

Wah, congratulations! I did fail in several tries, probably because all those little nifty oversights.
Reply
#12
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. Smile Which I suppose is useful as a neat gimmick, but not as a thing for calculating Tetration Smile

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));
}


Attached Files Image(s)
                   
Reply
#13
(09/28/2009, 10:33 PM)mike3 Wrote: ...It's like "waves" are being emitted from the center and rippling on out toward the ends...
What is wrong?
Mike, there is nothing wrong; you seem to do well.
you can make a beautiful animation with these waves.
I also had the similar waves in first versions of my code.

In order to get convergence and evaluate the tetrational,
update first the even nodes, going from the center to the ends,
and then the odd nodes, going from the ends to the center.
(In my case, the number of nodes is even, so, I can do both in the single loop,
but it is analogy of that I suggest above.)
This dumps the waves you mention.
Please, plot the first few iterations, and post them.
Reply
#14
Well I improved the code with this technique:

Code:
/* Do an iteration over a vector v representing values of tetration from -iA to +iA */
/* Odd number of nodes required */
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);
            local(node1=0);
            local(node2=0);
            local(ctrn=0);
            local(nidx1=0);
            local(nidx2=0);

            /* Update every other node from center to ends */
            ctrn = (N-1)/2+1;
            nidx1 = ctrn;
            nidx2 = ctrn;
            while(nidx1 >= 1,
                  zi = -I*A + (nidx1-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);
                  node1 = I1+I2+I3+I4;

                  zi = -I*A + (nidx2-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);
                  node2 = I1+I2+I3+I4;

                  newv[nidx1] = node1; /* simultaneously update both */
                  newv[nidx2] = node2;

                  expv[nidx1] = exp(node1);
                  expv[nidx2] = exp(node2);

                  logv[nidx1] = log(node1);
                  logv[nidx2] = log(node2);

                  nidx1 = nidx1 - 2;
                  nidx2 = nidx2 + 2;
               );

            /* Update remaining nodes from ends to center */
            if(ctrn%2 == 0,
               nidx1 = 1;
               nidx2 = N;,
               nidx1 = 2;
               nidx2 = N-1;
              );
            while(nidx1 <= nidx2,
                  zi = -I*A + (nidx1-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);
                  node1 = I1+I2+I3+I4;

                  zi = -I*A + (nidx2-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);
                  node2 = I1+I2+I3+I4;

                  newv[nidx1] = node1; /* simultaneously update both */
                  newv[nidx2] = node2;

                  expv[nidx1] = exp(node1);
                  expv[nidx2] = exp(node2);

                  logv[nidx1] = log(node1);
                  logv[nidx2] = log(node2);

                  nidx1 = nidx1 + 2;
                  nidx2 = nidx2 - 2;
               );

            return(newv);
}

I update the nodes starting from the center one then every other one out, the I go from the end nodes inward (or the next closest nodes to those end nodes if they were already done in the first pass). It seems to work good for node numbers like 211, 215, 219, etc. but not 209, 213, 217, etc. where it now diverges badly. Why is that?
Reply
#15
Mik, I looked at your code, and I do not understand how does it work.
1. If you have calculated the new value at some point z,
you have no need to calculate the same for z^*.
I have no compiler for the language you use, but I suspect,
it has operation of complex conjugation, why you do not use it?

2. Will you post the picture you got with this code?
It is easier to look than to read the code.

3.
Quote:Well I improved the code with this technique
You do not know, do you improve or make worse,
while you do not plot the residual.
Reply
#16
Well, I say "improved" because when it does converge it will converge for very high number of nodes. Perhaps though that isn't accurate, because it now doesn't work for numbers of nodes like 409, 413, 417, etc. I just realized that 209 is not enough, the node amount needs to be higher than that.

For node amounts like 411, 415, 419, etc. even very high ones, the code converges to the smooth graph, so I do not see the need to bother with graphing that case. At 811 nodes with A = 24 it converges exquisitely with the residual after 32 normalized followed by 6 unnormalized iterations having a magnitude smaller than 10^-9.

If you want, the graphs for case 409 nodes and A = 14 are shown below, after 6 to 10 iterations. After the first 5 it seems to be settling, then it does, well, you can see.

Finally, as for not using the conjugate, I do this because I'm thinking about attempting to use it on certain complex-number bases eventually, which will not have the conjugate symmetry. So I'd like to see if it is possible to get it to work without that.


Attached Files Image(s)
                   
Reply
#17
(10/01/2009, 07:57 PM)mike3 Wrote: Well, I say "improved" because when it does converge it will converge for very high number of nodes.
After some tens of iterations, do your curves become smooth?
You may have a problem related with the different weight of nodes in the Simpson,
the even nodes are twice "heavier" than the odd ones. The updating of each third node could boost the convergence, but the code becomes cumbersom. With the even amount of nodes I did not have such a problem, the weight of even nodes is similar to that of hte odd ones (except the tips, but there the weight of each node is small, so, it does no matter) and I could do it in a single loop.

Quote:..At 811 nodes with A = 24 it converges exquisitely with the residual after 32 normalized followed by 6 unnormalized iterations having a magnitude smaller than 10^-9.
Congratulations! I see, improving the precision, you go by the "extensive" way: you increase the number of nodes instead of to change for some high-order approximation;
each 3 addifional significant digits cost you an order of magnitude in the CPU time.
If you plan to get a precision better than that with complex<double> (15 digits),
try to finish your calculus in this century.
Another note: do you still evaluate the residual in few single points?
Could you plot the residual F(z)-exp(F(z-1)) in the complex z-plane?
Reply
#18
Yes, it turns smooth if the number of nodes meets the requirement I mentioned, otherwise it fails as the posted graphs show.

It's interesting you mention about the weight, though I do not update nodes during the Simpson procedure, and why does it seems that odd node numbers that converge alternate with those that don't?

As for using a higher order quadrature, yeah it would be better, but as you can see there are still a few issues to work out, then I'll switch to a more sophisticated method (cubic, maybe the Gauss-Legendre even). That's what I'm after here, trying to get out the bugs so they don't bite me later Smile
Reply
#19
(10/02/2009, 02:59 AM)mike3 Wrote: Yes, it turns smooth if the number of nodes meets the requirement I mentioned, otherwise it fails as the posted graphs show.
1. Can you explicitly formulate the requirement?
What value of A and how many nodes should one use in order to get 3 decimal digits, to get 9 decimal digits, to get 12 decimal digits, etc?

2. Will you print the table of values of your approximation of tetration along the imaginary axis?

Quote:It's interesting you mention about the weight, though I do not update nodes during the Simpson procedure, and why does it seems that odd node numbers that converge alternate with those that don't?
This may be due to the jumping distribution of weights in the nodes of the Simpson formula.
The Simpsons in the cartoon are better than in the precise evaluation of the tetrational.
Quote:..are still a few issues to work out, then I'll switch to a more sophisticated method (cubic, maybe the Gauss-Legendre even).
You already use the cubic one.
You may take the algorithm for the Gauss-Legendre at
http://en.citizendium.org/wiki/GauLegExample/code
The routine for the evaluation takes only 11 lines at very beginning of the code;
it is easy to translate to any language.
Quote:That's what I'm after here, trying to get out the bugs so they don't bite me later
I do not understand how do you treat the conjugation, some bug may be there.
if you do not use f(z^*)=f(z)^* ,
then the simultaneous update of F(z) and f(z^*) does not have any sense.

3. In order to catch bugs, calculate the residual at some sufficiently dense mesh and plot it. The mesh is really good tool to catch bugs.
Reply
#20
The explicit requirement, to be "guaranteed" of convergence, appears to be that the node number be an odd number like 211, 215, etc. i.e. be odd and also for (N-1)/2+1 to be even (or equivalently (N-1)/2 to be odd). One may be able to get away with other numbers if they are sufficiently small (I think 209 works, but 409 does not). Not sure where it starts to fail at precisely.

I'm not sure how calculating the values on the imag axis for where it succeeds in converging is useful for debugging this, though.

But if you really want it, here's what I get using 811 nodes (as this meets the given requirement) and 32 normalized followed by 6 unnormalized iterations with A = 24:
Code:
residual mag at 0: 0.0000000007856242620156700
residual mag at 2.2*I: 0.0000000009481759175558504
bias: 0.000000009638303087908689 - 1.135168310191088 E-16*I
tet(-2.2*I): 0.4620597753654270 - 1.295163164718386*I
tet(-2.1*I): 0.4798492754148191 - 1.283606567671135*I
tet(-2.0*I): 0.4993847992607353 - 1.269767235743538*I
tet(-1.9*I): 0.5207318688313386 - 1.253292769813972*I
tet(-1.8*I): 0.5439302090049257 - 1.233797070129863*I
tet(-1.7*I): 0.5689857178432989 - 1.210863129463098*I
tet(-1.6*I): 0.5958616741755516 - 1.184048130413531*I
tet(-1.5*I): 0.6244695115404595 - 1.152891336259040*I
tet(-1.4*I): 0.6546596527447695 - 1.116925217149086*I
tet(-1.3*I): 0.6862130782596795 - 1.075690135485824*I
tet(-1.2*I): 0.7188344751404906 - 1.028752709130253*I
tet(-1.1*I): 0.7521479533413454 - 0.9757276681381615*I
tet(-1.0*I): 0.7856963882146419 - 0.9163026217258810*I
tet(-0.9*I): 0.8189454132907991 - 0.8502646768632284*I
tet(-0.8*I): 0.8512929111838792 - 0.7775273409233825*I
tet(-0.7*I): 0.8820845103121792 - 0.6981556650513622*I
tet(-0.6*I): 0.9106350914594263 - 0.6123872292055688*I
tet(-0.5*I): 0.9362556726991349 - 0.5206464296113079*I
tet(-0.4*I): 0.9582843404237904 - 0.4235496900871865*I
tet(-0.3*I): 0.9761192253768883 - 0.3218997327286443*I
tet(-0.2*I): 0.9892510000207400 - 0.2166679077342638*I
tet(-0.1*I): 0.9972921071940174 - 0.1089647240256238*I
tet(0.0*I): 1.000000000000000 - 4.743384504624082 E-20*I
tet(0.1*I): 0.9972921071940174 + 0.1089647240256238*I
tet(0.2*I): 0.9892510000207400 + 0.2166679077342638*I
tet(0.3*I): 0.9761192253768883 + 0.3218997327286443*I
tet(0.4*I): 0.9582843404237904 + 0.4235496900871865*I
tet(0.5*I): 0.9362556726991349 + 0.5206464296113080*I
tet(0.6*I): 0.9106350914594263 + 0.6123872292055688*I
tet(0.7*I): 0.8820845103121792 + 0.6981556650513622*I
tet(0.8*I): 0.8512929111838792 + 0.7775273409233825*I
tet(0.9*I): 0.8189454132907991 + 0.8502646768632284*I
tet(1.0*I): 0.7856963882146419 + 0.9163026217258810*I
tet(1.1*I): 0.7521479533413454 + 0.9757276681381616*I
tet(1.2*I): 0.7188344751404906 + 1.028752709130253*I
tet(1.3*I): 0.6862130782596795 + 1.075690135485824*I
tet(1.4*I): 0.6546596527447695 + 1.116925217149086*I
tet(1.5*I): 0.6244695115404595 + 1.152891336259040*I
tet(1.6*I): 0.5958616741755516 + 1.184048130413531*I
tet(1.7*I): 0.5689857178432989 + 1.210863129463098*I
tet(1.8*I): 0.5439302090049257 + 1.233797070129863*I
tet(1.9*I): 0.5207318688313386 + 1.253292769813972*I
tet(2.0*I): 0.4993847992607353 + 1.269767235743538*I
tet(2.1*I): 0.4798492754148191 + 1.283606567671135*I
tet(2.2*I): 0.4620597753654270 + 1.295163164718386*I

And why does one need to bother with the residual, when one can plainly see the failure as the graph veers way off what it's supposed to be?

Also, isn't Simpson's rule based on a quadratic, not a cubic? At least that's what my calc. text said.
Reply




Users browsing this thread: 1 Guest(s)