1

I have an equation F(x,y;a,b,c,...)=0, where x and y are variables, a,b,c,... - parameters. y is a real number, x - complex. For every given y I need to solve this equation, i.e. to find Im(x) and Re(x) (this equation may have many solutions). After solving this equation I need to plot two graphs: the dependense of y from Im(x) and y from Re(x).

For example, F(x,y)=x^2 + sin(xy).

If y and x were real numbers, it would be possible to use ContourPlot. However, as my variables are complex, I have to solve this equation numerically and after that plot my graphs.

I have never come across such a problem, so could anyone tell me how to solve it?

David G. Stork
  • 41,180
  • 3
  • 34
  • 96
Artem
  • 43
  • 6
  • Solve, FindRoot ? Post the code you have written so far. – Sektor Apr 14 '15 at 19:49
  • Thanks. My code is too long, but I can explain. I tried to solve this problem in other way using ContourPlot: Manipulate[ContourPlot[Re[F(x,y,a,b)]==0,{x,0,5},{y,0,5}],{a,0,1},{b,0,1}] (or Im[F]). However, eventually I have realised that the unswer is wrong and, moreover, the calculations are too long and are not accurate at all. I understand how I should solve this problem (as I described above), but I have never done such things before, so I need advice. – Artem Apr 14 '15 at 20:07
  • Are you sure you know your way around Mathematica enough to tackle this problem ? – Sektor Apr 14 '15 at 20:08
  • I am not sure, that's why I have decided to ask. – Artem Apr 14 '15 at 20:15
  • I suggest you open up the documentation and read how to define functions, call them, pass different arguments, etc. – Sektor Apr 14 '15 at 20:44
  • 1
    Make the substitution x = x1 + I * x2 and solve for real x1 and x2. Of course, very soon you will have too many dimensions to plot on a flat computer screen :-) If you can obtain an analytical solution, try using Plot3D with y as a function of x1 and x2. – LLlAMnYP Apr 14 '15 at 20:45

1 Answers1

3

New way

The OP mentioned ContourPlot but its behavior is V10 makes my original solution practically unusable except for a very rough plot. Another approach is to solve the equation for all the roots in a given region. From the ContourPlot, one can see there are two types, ones that cross y == -5 and ones that cross y == 5. We can use NDSolve to solve the equation by first differentiating it (ode) and then integrating it (sols[y0]) using the intersections with the two planes as initial conditions (ics[y0]). This approach gives us one of the solutions twice, but we can delete by comparing values of x[y] at y == -5. This approach is much quicker and produces nice, differentiable, and accurate solutions.

eqn = x^2 + Sin[x y] == 0
ode = D[eqn /. x -> x[y], y]
(*
  x^2 + Sin[x y] == 0
  2 x[y] Derivative[1][x][y] + Cos[y x[y]] (x[y] + y Derivative[1][x][y]) == 0
*)

xmin = -5 - 5 I; xmax = 5 + 5 I; ymin = -5; ymax = 5;

ics[y0_] := NSolve[eqn && Re[xmin] <= Re[x] <= Re[xmax] && Im[xmin] <= Im[x] <= Im[xmax] /. y -> y0, x];

Clear[sols]; sols[y0_] := sols[y0] = Quiet[ First@NDSolve[{ode, x[y0] == (x /. #), WhenEvent[Re@x[y] < Re[xmin], "StopIntegration"], WhenEvent[Re@x[y] > Re[xmax], "StopIntegration"], WhenEvent[Im@x[y] < Im[xmin], "StopIntegration"], WhenEvent[Im@x[y] > Im[xmax], "StopIntegration"]}, x, {y, ymin, ymax}] & /@ ics[y0], NDSolve::ndsz]

tolerance = 8; mysols = DeleteDuplicates[Join[sols[ymin], sols[ymax]], Quiet[SetPrecision[x[ymin] /. #1, tolerance] == SetPrecision[x[ymin] /. #2, tolerance]] &];

Show[ ParametricPlot3D[{y, Re@x[y], Im@x[y]} /. #, Evaluate@Flatten[{y, x["Domain"] /. #}]] & /@ mysols, PlotRange -> All, AxesLabel -> {y, Re[x], Im[x]} ]

Mathematica graphics

Note that tolerance is used to determine how equal (up to what precision) two initial conditions have to be for them to be considered the same. In this case, they're identical, but I wanted to show a more robust solution.

Note also that as far as robustness goes, we made some assumptions about where the initial points on the solutions could be found.

Original answer (ContourPlot)

Something like this, perhaps? It's too slow for Manipulate, though.

plot = ContourPlot3D[
  {0 == Re[(x1 + I x2)^2 + Sin[(x1 + I x2) y]], 
   0 == Im[(x1 + I x2)^2 + Sin[(x1 + I x2) y]]},
 {y, -5, 5}, {x1, -5, 5}, {x2, -5, 5},
 ContourStyle -> None, Mesh -> None, MaxRecursion -> 2, 
 BoundaryStyle -> {1 -> None, 2 -> None, {1, 2} -> {{Lighter@Blue, Thick}}}, 
 AxesLabel -> Automatic]

Mathematica graphics

The BoundaryStyle trick is from this answer by Daniel Lichtblau.

Also, it seems faster in V9 than V10 (< 80 sec. vs. > 1629 sec. with MaxRecursion -> 1).


The parts of the coordinates corresponding to {Re[x], y} are {2, 1}. One can project the coordinates and ancillary options of Graphics3D to Graphics with the following.

With[{parts = {2, 1}},
 plot /.
  {GraphicsComplex[pts_, stuff__] :> GraphicsComplex[pts[[All, parts]], stuff],
   Graphics3D -> Graphics,
   HoldPattern[PlotRange -> pr_] :> (PlotRange -> pr[[parts]]),
   HoldPattern[AxesLabel -> label_] :> (AxesLabel -> label[[parts]])}
 ]

Mathematica graphics

To get {Im[x], y}, use parts = {3, 1}.

Michael E2
  • 235,386
  • 17
  • 334
  • 747
  • Thank you @Michael for the link for the BoundaryStyle trick - I have been trying to find it for several weeks now. (+1) – kglr Apr 14 '15 at 21:03
  • Thank you @Michael for this interesting idea. However, I'm not sure that this solution is appropriate for me, because you plot two graphs: Re(F[x,y]) = 0 and Im(F[x,y])=0, but I need the solution of the equation F[x,y]=0, i.e. Re(F[x,y])=0 and Im(F[x,y])=0 simultaneously. I mean, for example, I take y_1=0 and find x_1 = Im[x_1]+iRe[x_1] (maybe more solutions x_1', x_1'',...), then I take y_2=0.1 and find x_2 = Im[x_2]+iRe[x_2] (maybe also x_2',...) and so on. After that (using these numbers) I plot the dependense of y_i on Im[x_i] (also on Im[x_i'],...) and on Re[x_i] (again on Re[x_i'],...). – Artem Apr 14 '15 at 21:30
  • @Artem The plot is the intersection where the real and imaginary simultaneously vanish, i.e., where the complex value of F[x, y] is zero. Remove the option ContourStyle -> None and you will see the two surfaces, one for Re[F[x, y]] == 0 and one for Im[F[x, y]] == 0. You can get the actual points from the plot, if you need them. – Michael E2 Apr 14 '15 at 23:36
  • Thank you @Michael for your help. Can I further extract two 2D graphs from this 3D graph? I mean you plot F[Re[x],Im[x],y]=0, but I would like to get two projections of this graph Re[x],y and Im[x],y. – Artem Apr 15 '15 at 05:45
  • Dear Michael, thank you so much for your kind help! – Artem Apr 15 '15 at 13:23
  • @Artem You're welcome. :) – Michael E2 Apr 15 '15 at 13:27
  • Dear @MichaelE2, may I ask you another question? I have tried to use your code for a very simple function and I have received a strange result. (* x-[Sqrt[y^2+I*y]]==0,{y,0,1}], denote Re[Sqrt[y^2+I*y]]=A and Im[Sqrt[y^2+I*y]]=B, so Re[x]=A and Im[x]=B) Plot[{Re[Sqrt[y^2 + Iy]], Im[Sqrt[y^2 + I*y]]}, {y, 0, 1}]. The result which I receive using this straightforward method is essentially differet from the result I get using your code. (sorry, I don't know how to add pictures or .nb files here to demonstrate what I mean, but I can send it to you if you write your email). – Artem Apr 15 '15 at 15:40
  • Your Plot output and combined projections of the ContourPlot, for {y, -5, 5} and parts = {1, 2} and parts = {1, 3}. They seem to agree. – Michael E2 Apr 16 '15 at 02:40
  • Dear @Michael, I am sorry, this is my mistake. – Artem Apr 16 '15 at 17:11
  • Dear @MichaelE2, could you tell me what should I do in order to get "good" output? Mathematica accuratelu plots 3D graph, but backgroung of projections ir red (here is a link for this file). link – Artem Apr 16 '15 at 17:41
  • And could you also suggest how can I improve the quality of the graph without changing MaxRecursion -> 1, because my computer aborts calculations when I use MaxRecursion -> 2 – Artem Apr 16 '15 at 17:45
  • @Artem If you're willing to give up ContourPlot3D and manage a list of individual solutions, see my update. – Michael E2 Apr 16 '15 at 20:04
  • Dear @Michael, but could you tell me what version of Mathematica do you use? The thing is that I can't use your code with my version 7 Mathematica. – Artem Apr 17 '15 at 13:36
  • @Artem I'm on V10.1. (You might find it helpful mention you're using V7 in future questions, since a lot has changed since then.) Instead of WhenEvent, you have to use the "EventLocator" method in V7. Look it up in the docs for NDSolve or google it. I forget the exact syntax. You should be able to translate it. I think the rest of it should work. – Michael E2 Apr 17 '15 at 13:44
  • @Michael, ok, thank you! – Artem Apr 17 '15 at 15:06