hi @cswiercz sorry about my late response, I should be checking gitter more often
No worries. Haven't looked at csympy since my last post. Need the right environment setup. I've been working on "normal sympy" things.
I definitely want to get csympy working on all platforms. C++ is actually pretty portable, as long as one doesn't use any system dependent things, and we don't. We have added tests to CMake that check that your compiler is modern enough
It should work both with gcc and clang
Essentially, it either fails in cmake with a nice error message, or it doesn't and then
make must work. Otherwise that's a bug.
From what I last saw I think it might be something with cmake finding the required libraries. I'll poke around tomorrow and ping you with anything I find. I doubt it's something on csympy' send.
Btw, @cswiercz I independently discovered your work with Riemann surfaces. I had a few questions about that:
1) Is there a way in Sage or SymPy to plot a Riemann surface of a complex function, like f(z) = log(z), f(z) = sqrt(z^2-1), etc.?
2) What is the algorithm to visualize it? I read online that there are some kind of differential equations that one solves numerically, but I couldn't find any detailed descriptions of this.
3) How exactly is the Riemann surface defined? I guess once I understand 2), it will become clear. I understand that you move in the complex plane as x-y, and then you plot something on the z- axis, but clearly it's not just the Re(f) or |f|.
1) Is there a way in Sage or SymPy to plot a Riemann surface of a complex function, like f(z) = log(z), f(z) = sqrt(z^2-1), etc.?
2) What is the algorithm to visualize it? I read online that there are some kind of differential equations that one solves numerically, but I couldn't find any detailed descriptions of this.
3) How exactly is the Riemann surface defined? I guess once I understand 2), it will become clear. I understand that you move in the complex plane as x-y, and then you plot something on the z- axis, but clearly it's not just the Re(f) or |f|.
Wow! I didn't know there was interest! :)
Well, it comes from this long topic on sage-devel: https://groups.google.com/d/topic/sage-devel/6j-LcC6tpkE/discussion
1) To my knowledge, no.
2) The naive way is to do a complex plot of slices. However, the interpretation I take is to use local parameterizations of Riemann surfaces. I'll explain in 3:
3) I define a Riemann surface via a desingularization and compactification of a complex plane curve. In particular, I think of the surface as a y-covering of the x-sphere where C:f(x,y)=0 defines the curve...
2) The naive way is to do a complex plot of slices. However, the interpretation I take is to use local parameterizations of Riemann surfaces. I'll explain in 3:
3) I define a Riemann surface via a desingularization and compactification of a complex plane curve. In particular, I think of the surface as a y-covering of the x-sphere where C:f(x,y)=0 defines the curve...
3 cont.) However, the map from the Riemann surface to the curve is not 1-1 at the discriminant and singular points of the curve. At those points I use a Puiseux series expansion of the curve to separate the corresponding places on the Riemann surface.
Essentially one either treats complex functions as multivalued, but then one needs good understanding of Riemann surfaces (which I don't have), but then things like this hold:
log(a*b) = log(a)+log(b), or as single valued on a principal branch (I understand this very well, I put my notes here: http://www.theoretical-physics.net/dev/math/complex.html), but then the correct formula is log(a*b) = log(a) + log(b) + 2*pi*i*floor((pi-arg(a)-arg(b))/(2*pi)), i.e. there is a correction.
(A Puiseux series is a type of parameterization of the form (x(t),y(t)))
So I only work with finite genus Riemann surfaces. The nice thing about them is that they can all be written as the desingularization and compactification of a plane curve.
Though, there are other ways to parameterize them.
3 cont.) Anyway, you can plot these Puiseux series, giving a local picture of the surface at a place. That's the short of it.
Can your algorithm plot surfaces like these: http://commons.wikimedia.org/wiki/File:Riemann_surface_cube_root.jpg
So is my understanding correct, that Riemann surface is not very well defined and everybody defines it a bit differently?
Probably. At least close to the center. A Puiseux series is a kind of Taylor series but with some branching information in it. They unfortunately have pretty slow convergence rates.
I can write something up for you tomorrow in a gist or something.
I was hoping of having a routine, to which I pass a complex function and it (quickly) gives me a 3D representation of the Riemann surface, so that I can plot it.
I would create the Riemann surface obtained via the curve C:y3−x=0, compute the Puiseux series representation at the center, ...., and then profit. I can see how it would be done and it would probably be easier to demonstrate via code.
so if I have f(z) = log(z), how do I obtain the Riemann surface from your curve?
Ooooh. That's good to know. I'm very interested in visualization methods for Riemann surfaces. I have some preliminary ideas that I'd like to try out. Been working on other things in my package, though.
You can't. log(z) is infinite genus.
ok, let's do f(z) = z^(1/3)
- Define the algebraic curve C:f(x,y)=y3−x=0
(But I would like to be able to plot f(z) = log(z) and other similar functions too eventually.)
2) Compute the Puiseux series expansions centered at (x,y)=(0,0).
(Btw, somehow I have a feeling that this involves quite a bit of calculations, to obtain a 3D visualization, so it could be a great application of csympy.)
3) The series are given in terms of x, which I can vary on the complex disc. Plotting each of these series will give the "intersecting" surface.
Is the step 1. the same for other functions like f(z) = sqrt(z^2-1)
The Sympy / CSympy part is involved in computing the Puiseux series expansions of the curve / surface. I've jumped though a lot of hoops to try to make it as fast as possible. I'm really really looking forward to CSympy.
Once you try CSymPy out, let me know a list of things that you need implemented.
So that I can prioritize it.
For f(z)=√z2−1 I would define the curve C:y2−(x−1)=0 and Puiseux expand at (x,y)=(1,0).
As I mentioned I can try to write something up tomorrow. Send me some finite genus (algebraic) examples you'd like to see. It'll be a good way to test my code on some simple test cases.
Ok, I'll do so. I know wolfram alpha can do it, but it looks weird: http://www.wolframalpha.com/input/?i=riemann+surface+z%5E%281%2F2%29
it gives two surfaces, and they do not correspond to what I would imagine for sqrt(x).
I read the first 10 or so links in google about riemann surfaces and I still don't understand how it is defined.
w2−z is a two-sheeted surface. It looks like Mathematica is plotting each sheet individually...which is kinda weird because you need to make an arbitrary choice of where you branch cut lies.
Or maybe it's plotting the real and imaginary parts of w separately.
Probably, because x^(1/3) should have 3 sheets, shouldn't it? http://www.wolframalpha.com/input/?i=riemann+surface+z%5E%281%2F3%29
Yep!
And it only plots 3 images.
I wonder if this plot: http://commons.wikimedia.org/wiki/User:Jan_Homann/Mathematics#mediaviewer/File:Riemann_surface_sqrt.jpg is simply the absolute value of the real and imaginary parts that wolfram alpha returns, and then colored by a complex phase.
Actually, looking at the example you just sent I think it's the real and imaginary parts of w in w3=z. You can kinda see three separate values of w per z in the pictures.
The vertical direction is the real part of w. The coloring is the argument. The continuity of the color through the "supposed intersection" is supposed to demonstrate the the surface isn't really intersecting itself (it is a complex manifold, after all).
So could it be that Riemann surface is all possible complex values for the given function, e.g. for
f(z) = log(z), we get f(z) = log|z| + i*arg(z) + 2*pi*i*n for all integer "n", and that we just plot this in some way, for example we plot |f(z)| for all "n" and color by arg(f(z))?
Sort of. In a way the Riemann Surface Plotting Problem is a "how do I plot a four-dimensional thing" problem. The link between these two comes from the interpretation of an RS as trying to write w as a function of z from the equation f(z,w)=0. (For example, w3−z=0.) z is complex. However, w is complex as well. re(z)vs.im(z)vs.re(w) is a three-dimensional plot. But you're missing the imaginary part of w.
I haven't yet read any of the literature attempting to plot w in a way other than plotting the real and imaginary parts separately. I wouldn't be surprised if there are some clever attempts out there.
So, to answer you original question, my approach is to plot the real and imaginary parts of w separately.
Well but that's the same problem as how you plot any complex function, even if just the principal, single valued, branch.
True. But with Riemann surfaces you have the added complication of "unraveling" the apparent intersections...which, again, don't exist.
Right.
I didn't realize that Riemann surface is just all the complex values of the multivalued function
so for each "z" you get many (even infinitely many sometimes) values of f(z) and both "z" and f(z) is complex, so there is plenty of ways to visualize this in 3D or 2D.
I suppose you could do some kind of polar plot where, instead of ranging from 0 to 2π you would range from 0 to 2πB where B is the branching number of a place on the Riemann surface at which you're plotting. (Which is at most the w-degree of the curve.)
Exactly, I was thinking of introducing a 2D mesh in the complex x-y plane.
I might play with this approach in my plots.
then loop around it once and plot f(z) in 3D
then continue looping twice, and keep plotting f(z)
(i.e. by plotting f(z), I mean plotting |f(z)| and color by arg(f(z)))
then tree times etc.
for log(z), you just need to stop this by hand.
for z^(1/3) you arrive at the initial position in 3 loops.
Sure. Although, what I was suggesting is that a single rotation in the plot represents d rotations in the z-plane where d is the w-degree.
Otherwise, you get the intersections in the plot.
ok
I think you can see the 2D mesh in this plot: http://commons.wikimedia.org/wiki/User:Jan_Homann/Mathematics#mediaviewer/File:Riemann_surface_sqrt.jpg
So it seems to me you just need to evaluate at the vertices of the quads
That I can definitely do. Perhaps even without Pusieux series... (In fact, I did it once before but I, of course, ended up deleting that iPython notebook.)
and then I can use a tool like Paraview to plot the 3D surface.
Thanks for your time. I'll think about this.
We can work on this together.
Sure! (Though, I'm not as cool as you are. I still use matplotlib.)
Matplotlib is fine too.
It'd be nice to have a function that does this.
Like I said, I can ping you tomorrow with some plots and code. Possibly sans Puiseux series.
Excellent, that would be awesome!
Though, I would eventually like to have some sort of plotting routines for my
RiemannSurface object in abelfunctions.
Thanks for the chat!
Definitely.
Thanks for the chat as well.
@cswiercz here I found Matlab code that does this I think -- I just need to translate to ipython notebook. ;) http://people.math.sc.edu/sharpley/math552/Homework.html
and here is good introductory text for plotting of complex functions: https://www.tu-chemnitz.de/mathematik/preprint/2009/PREPRINT_12.pdf
See also plenty of examples of the phase plots here: http://commons.wikimedia.org/wiki/User:Jan_Homann/Mathematics
@certik Here's a preliminary demonstration. The key idea is to pass into polar coordinates in the complex z-plane. https://gist.github.com/cswiercz/1fde0a82f8e9e1b0660a
nice, this is it!
I tried to also plot abs(w) for the sqrt(z), but that doesn't look like the Riemann surface at all
See the cell [9] in http://nbviewer.ipython.org/gist/certik/f4b9990452ba8221211e
@cswiercz do you think that the Riemann surface in http://commons.wikimedia.org/wiki/File:Riemann_surface_sqrt.jpg is just the re(w)? It looks like it is, as the caption says "riemann surface of Sqrt[z], projection from 4dim C x C to 3dim C x Re(C), color is argument"
And the WA images (http://www.wolframalpha.com/input/?i=riemann+surface+z^%281%2F2%29) seem to be simply re(w) and im(w).
Color is currently not the argument, but the color is continuous across the "intersection" so it clearly shows which branch is which. I'll see if I can get color by argument working.
Let's try to create just a 2D plot by the argument. I wonder how to do this one: http://commons.wikimedia.org/wiki/File:Z4_over_z4minus1.jpg
As you can see, the color function is given there, and it is not just arg(z), it somehow also depends on "r"
I.e. here is the color function: http://commons.wikimedia.org/wiki/File:Z4_over_z4minus1.jpg#mediaviewer/File:Complex_coloring.jpg
it's not just a function of arg(z).
I can't figure out how that works.
Maybe you pick the color for the value "f(z)" from the color diagram simply by picking the x-y coordinates as re(f(z)) and im(f(z))
i.e. you pick the angle using arg(f(z)), but then depending on abs(f(z)) you pick how far you are from the center.
That makes a lot of sense to me.
I'm not particularly convinced by colored complex plots
they make for pretty pictures, but not much else
What's a better way @asmeurer?
I don't know
if you're going to do it, watch https://www.youtube.com/watch?v=rkDgBvT-giw
the choice of color map is very important
if you have a question you are trying to answer, you can maybe pick a good visualization, but I don't know of a good way to "overall" visualize four dimensions
Here is an interesting complex colorplot technique I came across a while ago. Also might be an option: http://nbviewer.ipython.org/github/empet/Math/blob/master/DomainColoring.ipynb
@cswiercz that notebook is very good. It shows that no single coloring is the best.
@asmeurer thanks for the video.
It seems to me that the right way to do this is to implement a few common options (e.g. all the methods in the notebook that @cswiercz just sent) and then allow the user to provide his own colormap.
Oh! The video is of Kristen! (I know her from UW.)
As to plotting of the Riemann surface, I can see now that sometimes you want to plot re(f(z)), sometimes im(f(z)) in the z-dimension.
@cswiercz is this your understanding as well?
And then you want to color it by any of the methods that we talked about for 2D plots.
there is also this talk that happened right before that one that was less informative but more entertaining https://www.youtube.com/watch?v=Alnc9E1RnD8
@certik Yes. Plotting Re(z) vs. Im(z) vs Re(w) or Re(z) vs. Im(z) vs Im(w) is what I've always known to be a complex 3d plot / Riemann surface plot.
I'm updating the gist I posted with a custom colormap option. Kinda messy at the moment but the idea is there.
@cswiercz can there be any other Riemann surface plot than those two? I am sorry for all these questions, I read http://en.wikipedia.org/wiki/Riemann_surface, but there is not a single word that you are just plotting Re(z), Im(z), Im(f(z)) or Re(z), Im(z), Im(f(z))
There's a discontinuity at the end of two full rotations because the jet colormap isn't periodic.
I know that people have attempted all sorts of clever ways to plot Riemann surfaces but I don't know any of those techniques. It's something that I've been meaning to investigate for a while but I keep getting caught up in other projects.
I see, as we talked about, it's essentially a plot of Re(z), Im(z), Re(f(z)), Im(f(z)), and then you are just trying to project it to 3D or 2D.
I checked your update, it looks very good
This paper discusses a related approach: http://www.mi.fu-berlin.de/en/math/groups/ag-geom/publications/db/2010-AutomaticGenerationOfRiemannSurfaceMeshes.pdf It's not too far beyond what we're already doing but it adds some more modulus information to the plots.
Let me throw a fly in the ointment for the "naive approach" I wrote up: this probably only works for Riemann surfaces in which there is only one place lying above the curve at (w.l.o.g) $z=0$.
Er, I guess I meant when there is more than one w-value lying above your center z=0. It may still work but it will be a messy picture.
(In the w2=z example there is exactly one w lying above z=0.)
I see. Yes, I read that paper yesterday too.
With Puiseux series you can separate the sheets by how you can get from one sheet to another. For example, you can have a Riemann surface where over the branch point z=0 you can have four w-sheets but you can only swap between sheets 1 and 2 and between 3 and 4 by rotating around the branch point. (i.e. there is no way to get from sheet 1 to sheet 3 by just going around that one branch point. Since the RS is a manifold you can get from sheet 1 to sheet 3 but only by going around some other branch point.)
Interesting.
Here you can plot a log(z) with your script: w = numpy.log(R) + 1j*Theta
it works like a charm.
and then you do branching_number = 3 to set how many revolutions you want to do
Neat!
So in most literature on Riemann surfaces (at least on finite genus ones) they always draw the cartoon of a genus g torus saying that it's the Riemann surface. I always wondered if it was possible to plot something like that. It would be especially challenging since the points / places at infinity are on that torus cartoon...
...anyway. Just a thought.
Can you post a picture?
I don't think I've seen a torus as a Riemann surface. Of which function?
@asmeurer anyway, @cswiercz code is really fast, so I think we should add this functionality to sympy
I would imagine just lambdifying the expression using numpy, and then use @cswiercz code.
For now, see slides 12-14 of my general exam talk: http://www.cswiercz.info/2014/03/31/general-exam.html I'll look for some better pictures.
There is a cute process by which you "glue together" the branch cuts on each sheet of your surface to demonstrate that it's homeomorphic to a genus g torus.
(Again, this is only for the finite genus case. In the infinite genus case like ew=z there are infinitely many branch cuts.)
@cswiercz you have a typo on your webpage "My reserach is in applications of algebraic geometry to nonlinear partial differential equations with my advisor, Bernard Deconinck."
I'm not sure what "reserach" is. It's probably what I've been doing for these past few years instead of "research".
Clearly nobody noticed
Because nobody is reading it. :)
Thanks for the typo catch, @asmeurer !
Okay. I have a whole mess of end-of-the-quarter homework to grade. I'll check back in later today.
@certik : The first few sections of Gamelin's "Complex Analysis" textbook illustrates the branch cutting and pasting operation fairly well.
@cswiercz thanks, I need to check it out. Here is the nb with the log(z) plot, cell [11]: http://nbviewer.ipython.org/gist/certik/f90c2bef899bcd01b7e7
Thanks @cswiercz for your help, this was very useful. I can now see how to code up and reproduce pictures from the literature.
One last thing. (Grading is boring.) Here's a document with similar sketches from Gamelin's text: http://math.berkeley.edu/~teleman/math/Riemann.pdf. Check out pages 1-5.
Also, your log plot looks really nice!
I've seen this pdf as well yesterday. ;)
As far as integrating into sympy, the only hard part is to figure out how to rewrite log(z) -> w = numpy.log(R) + 1j*Theta, so that it can be plotted. This should do it in this case:
In [4]: log(z).as_real_imag()
Out[4]: (log(│z│), arg(z))But we need something general that works for any expression.
In general, the strategy might be to do a change of variables z=reiθ.
With exponents: w=zp/q=rp/qeiθp/q.
With log: w=log(z)=log(r)+iθ
Compute
r and theta domains, transform to z and then plot?
This might be a start:
w = sqrt(z)
w = w.subs(z,r*exp(I*theta), positive=True)
w = sympy.expand_power_base(w, force=True)
w = sympy.powdenest(w)
print wOutput:
sqrt(r)*exp(I*theta/2)
Nice!
Then I blindly tried this and received an interesting error!
w = sqrt(z)
w = w.subs(z,r*exp(I*theta), positive=True)
w = sympy.expand_power_base(w, force=True)
w = sympy.powdenest(w)
R = sqrt(x**2+y**2)
Theta = sympy.atan2(y,x)
W = w.subs({r:R,theta:Theta})
Wreal = sympy.re(W)
Wimag = sympy.im(W)
sympy.plotting.plot3d(Wimag, (x,-1,1), (y,-1,1))/Users/cswiercz/anaconda/lib/python2.7/site-packages/sympy/plotting/experimental_lambdify.py:164: UserWarning: The evaluation of the expression is problematic. We are trying a failback method that may still work. Please report this as a bug.
warnings.warn('The evaluation of the expression is'...and then it produced the "naive" image with a jump...so that's not the right way to do it.
I like the phrase "The evaluation of the expression is problematic." :)
Wow.
@asmeurer ^^^