Tag Archives: Vepstas

Vepstas Regions 4

Further to the topic of Vepstas regions, defined for positive R as V(R) = the set of complex z for which f(z)=z^2/(z-1) has magnitude less than R. Linas Vepstas, in his software package ANANT for arbitrary precision number theoretical calculations, implements a method of calculating the polylogarithm of complex order s of complex argument z. He defines a complex data type cpx which is simply a pair of pointers to two arbitrary precision real values, ie the real and imaginary parts of a complex number. The routine to calculate polylog(s,z) is cpx_polylog in the mp_polylog.c module of ANANT.

Looking at the code for that module, one sees that, after some preliminary qualification of the argument, the evaluation for z in the unit circle (or on its boundary) comes down to a very neat algorithm. Here, to start things off, is a picture of the boundary of the V(1.225) region (1.225 is approximately square root of 1.5).

If z is in V(1.225), then polylog(s,z) is evaluated directly using an algorithm due to Peter Borwein, which I’ll call the PB algorithm. Otherwise, if z is in the unit circle U but not in V(1.225), a “duplication formula” for polylog is used to transform the calculation of polylog(s,z) to that of a linear combination of polylog(s,-z) and polylog(s,z^2). Suppose that z is complex, in U-V(1.225). Then z has positive real part at least 0.6, and so -z has negative real part, so lies in V(1.225) and polylog(s,-z) is readily calculated with the PB algorithm. Moreover, z^2 is at twice the phase angle of z, so swings up (or down) from the real axis. If z^2 now lies in V(1.225), then the PB algorithm can be used to calculate polylog(s,z^2), and the duplication formula can be used to obtain polylog(s,z) as a linear combination.

Otherwise, suppose that z^2 lies in U but not in V(1.225). Then use the duplication formula again, to express polylog(s,z^2) as a linear combination of polylogs of -z^2 and z^4. The phase angle of z^4 will be twice that of z^2. And so on. Eventually, one gets z^(2^n) which lies in V(1.225), and can complete the calculation. The whole process is a very nice illustration of how recursive coding can simplify a programming task conceptually.

Perhaps in another post I’ll discuss how polylog(s,z) is calculated for z outside the unit circle.

I guess I’d better say something here about the situation of real z between about 0.6 and 1.0. For z less than 1.0, successive applications of the duplication formula will produce z^2, z^4, etc tending towards the origin, so eventually get some z^(2^n) within V(1.225). For z equal to 1, the function polylog(s,1) is identical to the Riemann zeta function, so one can use some other code to calculate Riemann zeta.

Some references so you can follow up on this yourself:

Vepstas paper: Available at the URl http://arxiv.org/abs/math/0702243

ANANT code: Can be downloaded from https://launchpad.net/anant

Best wishes,
Ken Roberts
21-Apr-2014

Vepstas Regions 3

Leave a reply

Carrying on with the exploration of Vepstas regions. The function f(z) = z^2/(z-1). R is positive. V(R) is the set of complex numbers z for which abs(f(z)) < R. The boundary B(R) is the set of z for which abs(f(z)) = R. Prior post gave the solutions x (reals) for which z = x is on the boundary. If we imagine the B(R) curve drawn in the complex z-plane, those x values are the points where B(R) intersects the real axis.

Can we generalize to get a better understanding of the shape of the boundary curve B(R)? I can think of three ways right off:

(1) One might ask where some other line in the z-plane intersects the boundary curve B(R). For example, pass a line thru the origin, say the imaginary axis, or perhaps a line thru the origin at some other slope. Or maybe, because z=1 is a special point, being a pole of the function f(z) and hence guaranteed to be outside V(R), pass lines thru the point z=1. Or, to be a bit more general, pass non-straight curves thru the origin or thru the point z=1, and ask where they intersect the boundary B(R). All those involve considering some function g(z)=0 which defines the straight line or curve, and solving the simultaneous equations abs(f(z))=R and g(z)=0. This seems quite feasible, but perhaps tedious. It is like attempting to determine the shape of some hidden object, say one wrapped in a huge ball of yarn, by passing a knitting needle into the ball. If you think that sounds a bit like computer axial tomography (CAT) scans, or maybe particle scattering experiments, you are correct. This is in general known as an inverse problem — to attempt to characterize a shape by looking at how things bounce off it or the absorption (loss of energy) they experience when passing thru it. There are many interesting aspects of such problems. Maybe another time. Personally, I'm going to set method 1 aside, with regard to the Vepstas Region boundary determination, and look for a simpler approach.

(2) Another method might be to prepare a 3-D graph of the function h(z) = abs(f(z)). That is a real valued function of a complex variable, so is quite easily graphed if one has suitable software. Especially software which allows one to rotate the 3-D graph to look at it from various points of view. Each R value will produce a contour line, which represents the boundary region B(R). Thus one might also want to use a contouring software if available, and see all the R-contours nested one within another. I am going to set that method aside as well, for the present, and assume such 3-D graphing or coutouring software might not be readily available to you. We can perhaps look at the 3-D graph and R-contours of the function h(z) in a later post. Right now let's try to get something further done by alternative methods.

(3) The method I'm going to pursue for the present, is an extension of the method used to solve abs(f(x))=R for x on the real line, in the prior two posts. This third method is actually not the only alternative — there are more than three ways to determine the B(R) curves — but let's stick with this method 3 for now. That involved writing f(x)=x^2/(x-1) where x is real, and asking what would make the absolute value of f(x) equal to R. There are two possibilities: f(x) is either +R or -R, which was written for convenience as f(x)=pR where p denotes either +1 or -1. Then the equation x^2/(x-1)=pR was multiplied out, giving a quadratic equation in x, which could be solved to get a formula for x in terms of R and p.

Suppose that x were not restricted to be a real number, but could be an arbitrary complex number. Then the equation abs(f(x))=R is equivalent to f(x)=pR, where now we are letting p represent any complex number of magnitude 1. That is, p is any point on the unit circle U about the origin. We have x^2/(x-1)=pR, and just as before, we get the quadratic equation x^2-pRx+pR=0, which can be solved using the quadratic formula to obtain x=(1/2)*(pR+sqrt(p^2*R^2-4pR)). One change here is that I’ve written ‘sqrt” to denote the complex square root function. The complex square root is a multi-valued function, two values in fact, and hence the formula for the solution actually represents two values, one for each possible value of the sqrt function. The two-valued interpretation of the sqrt function replaces the use of the variable q which, in the real-x case, represented +1 or -1.

We can no longer remove p^2 from that formula, because p^2 for p on the unit circle is another point on that circle, not just 1. We can talk about the phase angle t of the value p (the point p), by writing p=exp(it). Then p^2 is the point which is twice as far rotated around the circle, ie has phase angle 2t, so p^2=exp(2it). Within the square root, the term p^2*R^2-4pR represents the result of taking the real number R^2 (bigger than 4R if R is greater than 4, and lesser than 4R if R lies between 0 and 4), rotating it by phase angle 2t, and then subtracting 4pR, which is rotated only by phase angle p. The resulting complex number D=p^2R^2-4pR is something with a phase angle intermediate between the two terms p^2R^2 and 4pR. The magnitude of D is also intermediate between those two terms. Calculating sqrt(D) gives two results. One has a phase angle half of that of D, so a phase angle less than t. And the other lies opposite the first, ie with phase angle t+pi. (Oops — pi here represents 3.14… not symbol p times symbol i. Please make allowances for unfortunate choice of letter p in the foregoing, combined with use of “pi” rather than greek letter pi!). And the magnitude of sqrt(D) is between 1 and abs(D). It is all perfectly definite when done geometrically. It’s only the words that produce trouble when describing how to find D and sqrt(D).

So then what’s x=(pR+sqrt(D))/2 look like? It is the midpoint of the line between the points pR and sqrt(D). Remembering that there are two points sqrt(D) because sqrt is two-valued, so we are actually considering two midpoints x.

Now things are getting interesting. We are almost able to make a geometric construction for the solution and to build a device — a set of linkages and sliders — that will draw the boundary set B(R). Let’s start by replacing the expression pR in the x-formula with a variable S, where S=pR. The variable S represents a point on the circle of radius R about the origin. We can imagine S moving around that circle, being moved by hand or by machine, as the device (still to be designed in detail) traces out the B(R) curve. The formula for x, with the substitution S=pR, becomes x = (S+sqrt(S^2-4S))/2. We already understand how to do the (S+sqrt(D))/2 part; it is the midpoints of the lines between S and the two points represented by sqrt(D). How to calculate S^2? It is at twice the phase angle t of S, and has magnitude equal to R^2, so is the point on the circle of radius R^2 which lies at phase angle 2t. From that we subtract 4S. The point 4S will be the point on the circle of radius 4R which lies at the phase angle t. We subtract vector B from vector A by adding the negative of vector B to vector A. So we want to have some linkage which traces out -4S; that is easy, just extend the vector S 4 times its length in the other direction from the origin. Now add the vectors S^2 and -4S. That vector addition is achieved by using the parallelogram law; we want a linkage for that. I’m going to assume such can be built. I don’t know how to do that at the moment, but believe I know where to look. I’ll hold the details of mechanical construction for such a vector-adding device to a later post.

Here is an outline of how to construct the B(R) curve.
(a) Pick a point S on the R-circle, ie a vector S in the complex plane. Let t be the phase angle of S, that is, the angle between the S vector and the positive real axis, measured counterclockwise.
(b) Construct the point A = S^2 on the circle of radius R^2 which has phase angle 2t.
(c) Construct the point C = -4S on the circle of radius 4R which has phase angle t+pi, that is, in the direction opposite to the vector S.
(d) Use the parallelogram law to find the vector D=A+C=S^2-4S.
(e) Take the square root Q of abs(D). These vectors have magnitude equal to the square root of the distance of D from the origin; finding square roots is a standard ruler and compass construction. Also bisect the angle between the positive real axis and D, to determine the phase angle u which is half the phase angle of D.
(f) Locate the points Q1 and Q2 which are on the circle of radius Q about the origin, and have phase angles equal to half the phase angle of D (that’s the point Q1) and to pi plus half the phase angle of D (that’s the point Q2). The point Q2 is opposite the point Q1, the same distance Q from the origin.
(g) Construct the midpoints of the lines between S and Q1, and between S and Q2. Those two midpoints lie on the boundary B(R) of the Vepstas region V(R). Let’s call those points B1 and B2.
(h) Repeat steps (a) thru (g) as for as many of points S on the R-circle as is desirable to draw B(R) to suit.

It’s about time for a diagram, don’t you think? Here we are, first showing the construction for one particular S point. The R value is 3.9, and the S point is at phase angle 80 degrees. The two points B1 and B2 are labelled.

And here is a diagram showing the B(3.9) boundary. This has been determined by using the construction technique for 100 points S evenly spaced around the R-circle, at 3.6 degree increments of the phase angle.

I didn’t know, when starting this post, that it would lead to such an enjoyable geometric construction. Somehow I thought I was going to end up with some sort of algebraic transformation of the f(z) formula. So there is perhaps more to be said about Vepstas regions — of course there is! — but I’ll leave that for later posts.

Best wishes,
Ken Roberts
15-Apr-2014

Vepstas Regions 2

Leave a reply

Continuing with the topic of Vepstas regions … Did you spot my mistake in the prior post? To recap, f(z) = z^2/(z-1), where z is complex. Given positive R, the set of z for which abs(f(z) < R is called the Vepstas region V(R), and the boundary set B(R) is the set of z for which abs(f(z)) = R. I'm interested in the shape of Vepstas regions. As a preliminary step, I want to calculate the real values x = x(R) for which B(R) intersects the real axis in the z-plane. That means solving the equation f(x) = x^2/(x-1) = +R or -R, since we are restricting to the case of x being real. But I missed one of the solutions. Too many plus/minus alternatives on the go, I guess. So let's try again at solving abs(f(x)) = abs(x^2/(x-1)) = R, restricting x to be real.

To be systematic, let's use p to represent either of +1 or -1, so that we are solving, to find real x for which x^2/(x-1) = pR. That takes in both alternatives for +R and -R. We can think of p as being one of the two square roots of 1. Now write that equation as a quadratic in the form x^2 = pRx – pR, or x^2 – pRx + pR = 0. Using the quadratic formula, the solutions can be written as (1/2) times (pR plus or minus the square root of a quantity D). The quantity D is called the discriminant), and in this case the discriminant is D = p^2R^2 – 4pR. We can use a symbol q to represent either of +1 or -1, so that taking the square root of D will produce a result that can be written q*sqrt(D), and we do not have to say "plus or minus" when reading the equation. The solution (for real x) is x = (1/2)*(pR + q*sqrt(D)). Moreover, p^2 in the discriminant D is equal to 1 (because p is a square root of 1), and so D equals R^2 – 4pR. So the resulting solution simplifies to x = (1/2)*(pR + q*sqrt(R^2 – 4pR)). In this solution, each of p and q represents either of +1 or -1. There are four alternatives; p might be +1 whereas q is -1, or vice versa, or both p and q could be +1, or both could be -1. On the other hand, the two appearances of p in the solution are the same value. It seems a bit cumbersome to write symbols p and q, but one does not have to consider separate cases quite so much; some of the work, manipulating the quadratic equations, is combined into a common portion of the reasoning. Something like computer programming.

Here are the graphs of the x-intercepts of the boundary curves B(R) with the x-axis, for various R values.

If you compare with the prior post, you will see that there is an additional curve at the top. For R greater than 4, the boundary B(R) of a Vepstas region will intersect the x-axis four times. For example, if R is 5, the boundary curve B(5) will intersect the x-axis at x equal to any of four values, approximately x = -6, +3.7, 0.8 and 1.4 — just approximate, reading off the graph. You can calculate precise values from the formula for x, taking the various alternatives for p and q.

Now, about mistakes … how did I know I made a mistake in my prior post? It just didn’t feel right, but why? It seemed strange to have three solutions (for general R larger than 4), not four solutions. Walking along the x-axis, from very large negative x to very large positive x, one expects to be outside V(R), then inside by the time you reach the origin. Then outside as you approach x=1, where f(x) is undefined (plus or minus infinity). Then (for R bigger than 4) inside again as you reach x = 2, then outside once again for very large positive x. Outside, inside, …(some intermediate wobbles perhaps) … outside — that should involve an even number of transitions through the boundary B(R). So it was the kinethestics of the solution that bothered me. I don’t know how you feel about complex variables and complex functions, but for me they are genuinely FELT — movements of points and surfaces. When one considers the map z –> w = z^2 between two complex planes, it is a wrapping (a physical wrapping) of a surface (the z-plane) around the origin of the w-plane, and it must wrap around twice. You could do it with bedsheets, in principle, except for certain complexities as the sheet has to connect back to itself after two trips around the origin. Anyway, all I can recommend is to cultivate a geometric or physical “feeling” for complex variables and complex functions. For me, it is kinesthetic, and probably involves movement of my muscles (at least the hands and arms) when I am visualizing something.

More later about Vepstas regions.

Best wishes,
Ken Roberts
12-Apr-2014

Vepstas Regions

Leave a reply

The term Vepstas Region is a phrase I’m using for a region of the complex plane which appears in a paper by Linus Vepstas. It’s my term, not his! I just need something to use when speaking about regions of this type. A Vepstas region V(R) is the set of complex points z such that abs(z^2/(z-1)) is less than some constant R, let’s say R = 4 typically. The points of V(4) are the ones for which a computational algorithm devised by Vepstas will converge. The points of V(R) for R meaningfully less than 4, let’s say V(3) or V(2), are the ones for which his algorithm will converge raapidly. In this cluster of posts, I just want to play around a bit with the regions themselves, to look at their shapes. To have some fun with the underlying geometric objects, not worry about convergence of an algorithm.

Let f(z) = z^2/(z-1). The boundary of V(R) is the set B(R) for which abs(f(z)) = R. To get an initial handle on the regions, note that as z tends to 1 (in complex plane), f(z) behaves like 1/(z-1) so has a simple pole at z = 1. The point z = 1 is not going to be in any V(R) region. On the other hand, f(0) = 0 so the point z = 0 is going to be in V(R) for any positive R. The regions V(R) are not empty, and they have boundary points (in the set B(R)) along any path from z = 0 to z = 1. I guess I should mention that f(z) is analytic (except at z = 1), if you haven’t remembered that already, so it is continuous, has derivatives of all orders, is a conformal map, has a power series representation about any point in its domain (that is, except at z = 1 which is not in the domain of f(z)), and so on. The usual properties of analytic functions.

Given fixed positive R, there is a point (or points) z = x along the real axis, between 0 and 1, for which the absolute value of f(x) is R. Let’s find those points. Since we’re using only real values of z = x for the moment, the value of f(x) will also be real, so there are only two possibilities: f(x) is either +R or -R.

If f(x) = +R, we have to solve x^2/(x-1) = R, or x^2 = Rx – R, or x^2 – Rx + R = 0, which is quadratic, and has the solution x(R) = (1/2) * (R plus or minus square root of (R^2 – 4R)). In order to get a real result from that formula, we must have R^2 – 4R non-negative; that is there is a solution f(x) = +R only for R at least 4. Here’s a little graph of those solutions, showing x(R) as a function of R, for R at least 4:

Well, it appears that none of those points x(R) which solve f(x) = +R lies on the segment (0,1) of the real axis. They are all further out than 1. For large R, they get close to 1, never reaching 1. For R a little bit more than 4, they get close to 2, and do reach 2. So the boundary B(R) intersects the real axis somewhere between 1 and 2, provided R is at least 4. Let’s just check: f(2) is 2^2/(2-1) = 4, so that’s right; the point x = 2 is in V(R) for R equal to or greater than 4.

Now, what about the other case, f(x) = -R. Does that have solutions for x between 0 and 1? If f(x) = -R, we have to solve x^2/(x-1) = -R, or x^2 = -Rx + R, or x^2 + Rx – R = 0, which is quadratic, and has the solution x(R) = (1/2) * (-R plus or minus square root of (R^2 + 4R)). That formula will always give a real result, two of them. The square root of (R^2 + 4R) will be greater than R, so one of the results will be negative (outside the interval (0,1), and the other will be positive, and equal to (1/2)*(sqrt(R^2+4R) – R). Here’s a graph of those x(R) values. This time, I’m going to consider any positive R, not just R at least 4.

Actually, I’ve shown three x(R) curves on that graph. The top curve is the previous one, squished. The horizontal line is just to show the line with ordinate 1, so you can see the asymptotic behaviour. The middle curve is the position of the boundary point of B(R) on the x-axis between 0 and 1, and the lower curve is the position of the boundary point of B(R) for negative x.

That’s enough for now. I encourage you to play a bit with determining the shape of the V(R) regions for various R values. More later…

Best wishes,
Ken Roberts
10-Apr-2014

ps. Are there mistakes in the above? Perhaps. I make lots of mistakes; it is how I learn. See post some time ago about Mistakes being an essential aspect of the learning process. You will learn the most if you catch me in a mistake before I get around to correcting it (if ever).