These are chat archives for symengine/symengine

4th
Dec 2014 Ondřej Čertík
@certik
Dec 04 2014 05:11 UTC
hi @cswiercz sorry about my late response, I should be checking gitter more often Chris Swierczewski
@cswiercz
Dec 04 2014 05:13 UTC
No worries. Haven't looked at csympy since my last post. Need the right environment setup. I've been working on "normal sympy" things. Ondřej Čertík
@certik
Dec 04 2014 05:15 UTC
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. Chris Swierczewski
@cswiercz
Dec 04 2014 05:19 UTC
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. Ondřej Čertík
@certik
Dec 04 2014 05:19 UTC
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|. Chris Swierczewski
@cswiercz
Dec 04 2014 05:23 UTC
Wow! I didn't know there was interest! :) Ondřej Čertík
@certik
Dec 04 2014 05:24 UTC
Well, it comes from this long topic on sage-devel: https://groups.google.com/d/topic/sage-devel/6j-LcC6tpkE/discussion Chris Swierczewski
@cswiercz
Dec 04 2014 05:25 UTC
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...
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. Ondřej Čertík
@certik
Dec 04 2014 05:27 UTC
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. Chris Swierczewski
@cswiercz
Dec 04 2014 05:27 UTC
(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. Ondřej Čertík
@certik
Dec 04 2014 05:31 UTC
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? Chris Swierczewski
@cswiercz
Dec 04 2014 05:33 UTC
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. Ondřej Čertík
@certik
Dec 04 2014 05:35 UTC
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. Chris Swierczewski
@cswiercz
Dec 04 2014 05:35 UTC
I would create the Riemann surface obtained via the curve C:y3x=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. Ondřej Čertík
@certik
Dec 04 2014 05:36 UTC
so if I have f(z) = log(z), how do I obtain the Riemann surface from your curve? Chris Swierczewski
@cswiercz
Dec 04 2014 05:36 UTC
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. Ondřej Čertík
@certik
Dec 04 2014 05:37 UTC
ok, let's do f(z) = z^(1/3) Chris Swierczewski
@cswiercz
Dec 04 2014 05:37 UTC
1. Define the algebraic curve C:f(x,y)=y3x=0 Ondřej Čertík
@certik
Dec 04 2014 05:37 UTC
(But I would like to be able to plot f(z) = log(z) and other similar functions too eventually.) Chris Swierczewski
@cswiercz
Dec 04 2014 05:37 UTC
2) Compute the Puiseux series expansions centered at (x,y)=(0,0). Ondřej Čertík
@certik
Dec 04 2014 05:38 UTC
(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.) Chris Swierczewski
@cswiercz
Dec 04 2014 05:38 UTC
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. Ondřej Čertík
@certik
Dec 04 2014 05:39 UTC
Is the step 1. the same for other functions like f(z) = sqrt(z^2-1) Chris Swierczewski
@cswiercz
Dec 04 2014 05:39 UTC
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. Ondřej Čertík
@certik
Dec 04 2014 05:40 UTC
Once you try CSymPy out, let me know a list of things that you need implemented.
So that I can prioritize it. Chris Swierczewski
@cswiercz
Dec 04 2014 05:40 UTC
For f(z)=z21 I would define the curve C:y2(x1)=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. Ondřej Čertík
@certik
Dec 04 2014 05:42 UTC
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. Chris Swierczewski
@cswiercz
Dec 04 2014 05:44 UTC
w2z 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. Ondřej Čertík
@certik
Dec 04 2014 05:45 UTC
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 Chris Swierczewski
@cswiercz
Dec 04 2014 05:45 UTC
Yep! Ondřej Čertík
@certik
Dec 04 2014 05:45 UTC
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. Chris Swierczewski
@cswiercz
Dec 04 2014 05:46 UTC
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). Ondřej Čertík
@certik
Dec 04 2014 05:48 UTC
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))? Chris Swierczewski
@cswiercz
Dec 04 2014 05:52 UTC
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, w3z=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. Ondřej Čertík
@certik
Dec 04 2014 05:53 UTC
Well but that's the same problem as how you plot any complex function, even if just the principal, single valued, branch. Chris Swierczewski
@cswiercz
Dec 04 2014 05:54 UTC
True. But with Riemann surfaces you have the added complication of "unraveling" the apparent intersections...which, again, don't exist. Ondřej Čertík
@certik
Dec 04 2014 05:55 UTC
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. Chris Swierczewski
@cswiercz
Dec 04 2014 05:56 UTC
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.) Ondřej Čertík
@certik
Dec 04 2014 05:56 UTC
Exactly, I was thinking of introducing a 2D mesh in the complex x-y plane. Chris Swierczewski
@cswiercz
Dec 04 2014 05:56 UTC
I might play with this approach in my plots. Ondřej Čertík
@certik
Dec 04 2014 05:57 UTC
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. Chris Swierczewski
@cswiercz
Dec 04 2014 05:58 UTC
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. Ondřej Čertík
@certik
Dec 04 2014 05:59 UTC
ok
So it seems to me you just need to evaluate at the vertices of the quads Chris Swierczewski
@cswiercz
Dec 04 2014 06:00 UTC
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.) Ondřej Čertík
@certik
Dec 04 2014 06:00 UTC
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. Chris Swierczewski
@cswiercz
Dec 04 2014 06:02 UTC
Sure! (Though, I'm not as cool as you are. I still use matplotlib.) Ondřej Čertík
@certik
Dec 04 2014 06:02 UTC
Matplotlib is fine too.
It'd be nice to have a function that does this. Chris Swierczewski
@cswiercz
Dec 04 2014 06:02 UTC
Like I said, I can ping you tomorrow with some plots and code. Possibly sans Puiseux series. Ondřej Čertík
@certik
Dec 04 2014 06:02 UTC
Excellent, that would be awesome! Chris Swierczewski
@cswiercz
Dec 04 2014 06:03 UTC
Though, I would eventually like to have some sort of plotting routines for my RiemannSurface object in abelfunctions.
Thanks for the chat! Ondřej Čertík
@certik
Dec 04 2014 06:03 UTC
Definitely.
Thanks for the chat as well. Ondřej Čertík
@certik
Dec 04 2014 06:16 UTC
@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 Ondřej Čertík
@certik
Dec 04 2014 06:46 UTC
See also plenty of examples of the phase plots here: http://commons.wikimedia.org/wiki/User:Jan_Homann/Mathematics Chris Swierczewski
@cswiercz
Dec 04 2014 15:44 UTC
@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 Ondřej Čertík
@certik
Dec 04 2014 16:34 UTC
nice, this is it! Ondřej Čertík
@certik
Dec 04 2014 16:39 UTC
I tried to also plot abs(w) for the sqrt(z), but that doesn't look like the Riemann surface at all
@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). Chris Swierczewski
@cswiercz
Dec 04 2014 16:45 UTC
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. Ondřej Čertík
@certik
Dec 04 2014 16:50 UTC
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"
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. Aaron Meurer
@asmeurer
Dec 04 2014 16:54 UTC
I'm not particularly convinced by colored complex plots
they make for pretty pictures, but not much else Ondřej Čertík
@certik
Dec 04 2014 16:55 UTC
What's a better way @asmeurer? Aaron Meurer
@asmeurer
Dec 04 2014 16:56 UTC
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 Chris Swierczewski
@cswiercz
Dec 04 2014 16:59 UTC
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 Ondřej Čertík
@certik
Dec 04 2014 17:04 UTC
@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. Chris Swierczewski
@cswiercz
Dec 04 2014 17:07 UTC
Oh! The video is of Kristen! (I know her from UW.) Ondřej Čertík
@certik
Dec 04 2014 17:07 UTC
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. Aaron Meurer
@asmeurer
Dec 04 2014 17:08 UTC
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 Chris Swierczewski
@cswiercz
Dec 04 2014 17:10 UTC
@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. Ondřej Čertík
@certik
Dec 04 2014 17:12 UTC
@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)) Chris Swierczewski
@cswiercz
Dec 04 2014 17:13 UTC
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. Ondřej Čertík
@certik
Dec 04 2014 17:15 UTC
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 Chris Swierczewski
@cswiercz
Dec 04 2014 17:16 UTC
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.) Ondřej Čertík
@certik
Dec 04 2014 17:21 UTC
I see. Yes, I read that paper yesterday too. Chris Swierczewski
@cswiercz
Dec 04 2014 17:25 UTC
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.) Ondřej Čertík
@certik
Dec 04 2014 17:26 UTC
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 Chris Swierczewski
@cswiercz
Dec 04 2014 17:27 UTC
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. Ondřej Čertík
@certik
Dec 04 2014 17:30 UTC
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. Chris Swierczewski
@cswiercz
Dec 04 2014 17:33 UTC
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.) Aaron Meurer
@asmeurer
Dec 04 2014 17:35 UTC
@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." Chris Swierczewski
@cswiercz
Dec 04 2014 17:36 UTC
I'm not sure what "reserach" is. It's probably what I've been doing for these past few years instead of "research". Ondřej Čertík
@certik
Dec 04 2014 17:37 UTC
Clearly nobody noticed Chris Swierczewski
@cswiercz
Dec 04 2014 17:37 UTC
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. Ondřej Čertík
@certik
Dec 04 2014 17:43 UTC
@cswiercz thanks, I need to check it out. Here is the nb with the log(z) plot, cell : 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. Chris Swierczewski
@cswiercz
Dec 04 2014 18:39 UTC
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! Ondřej Čertík
@certik
Dec 04 2014 18:41 UTC
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 : log(z).as_real_imag()
Out: (log(│z│), arg(z))

But we need something general that works for any expression. Chris Swierczewski
@cswiercz
Dec 04 2014 18:49 UTC
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? Chris Swierczewski
@cswiercz
Dec 04 2014 19:05 UTC

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 w

Output:

sqrt(r)*exp(I*theta/2) Ondřej Čertík
@certik
Dec 04 2014 19:12 UTC
Nice! Chris Swierczewski
@cswiercz
Dec 04 2014 19:21 UTC

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." :) Ondřej Čertík
@certik
Dec 04 2014 19:21 UTC
Wow.
@asmeurer ^^^ Ondřej Čertík
@certik
Dec 04 2014 23:56 UTC