These are chat archives for symengine/symengine

4th
Dec 2014
Ondřej Čertík
@certik
Dec 04 2014 05:11
hi @cswiercz sorry about my late response, I should be checking gitter more often
Chris Swierczewski
@cswiercz
Dec 04 2014 05:13
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
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
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
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
Wow! I didn't know there was interest! :)
Ondřej Čertík
@certik
Dec 04 2014 05:24
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
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
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
(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
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
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
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
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
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
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
ok, let's do f(z) = z^(1/3)
Chris Swierczewski
@cswiercz
Dec 04 2014 05:37
  1. Define the algebraic curve C:f(x,y)=y3x=0
Ondřej Čertík
@certik
Dec 04 2014 05:37
(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
2) Compute the Puiseux series expansions centered at (x,y)=(0,0).
Ondřej Čertík
@certik
Dec 04 2014 05:38
(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
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
Is the step 1. the same for other functions like f(z) = sqrt(z^2-1)
Chris Swierczewski
@cswiercz
Dec 04 2014 05:39
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
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
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
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
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
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
Yep!
Ondřej Čertík
@certik
Dec 04 2014 05:45
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
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
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
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
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
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
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
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
Exactly, I was thinking of introducing a 2D mesh in the complex x-y plane.
Chris Swierczewski
@cswiercz
Dec 04 2014 05:56
I might play with this approach in my plots.
Ondřej Čertík
@certik
Dec 04 2014 05:57
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
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
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
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
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
Sure! (Though, I'm not as cool as you are. I still use matplotlib.)
Ondřej Čertík
@certik
Dec 04 2014 06:02
Matplotlib is fine too.
It'd be nice to have a function that does this.
Chris Swierczewski
@cswiercz
Dec 04 2014 06:02
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
Excellent, that would be awesome!
Chris Swierczewski
@cswiercz
Dec 04 2014 06:03
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
Definitely.
Thanks for the chat as well.
Ondřej Čertík
@certik
Dec 04 2014 06:16
@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
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
@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
nice, this is it!
Ondřej Čertík
@certik
Dec 04 2014 16:39
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
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
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
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
What's a better way @asmeurer?
Aaron Meurer
@asmeurer
Dec 04 2014 16:56
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
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
@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
Oh! The video is of Kristen! (I know her from UW.)
Ondřej Čertík
@certik
Dec 04 2014 17:07
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
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
@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
@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
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
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
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
I see. Yes, I read that paper yesterday too.
Chris Swierczewski
@cswiercz
Dec 04 2014 17:25
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
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
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
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
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
@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
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
Clearly nobody noticed
Chris Swierczewski
@cswiercz
Dec 04 2014 17:37
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
@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.
Chris Swierczewski
@cswiercz
Dec 04 2014 18:39
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
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.

Chris Swierczewski
@cswiercz
Dec 04 2014 18:49
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

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
Nice!
Chris Swierczewski
@cswiercz
Dec 04 2014 19:21

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
Wow.
@asmeurer ^^^