Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have an equation for function F[x,y]==0 which first argument x is real and another, y, is complex. All constants in equation are real.

I define y as y=a+I*b.

I need to plot Im[y] as function of x. (in other words, I need a plot of b as function of x. The equality for F for a and b and x can be solved only using FindRoot function.

Here's how I tried to plot what I need:

ContourPlot[{{Im[F] == 0, Re[F] == 0}}, {b, 0.1, 
  10^8}, {x, 0.1, 10^16}]

But nothing is plotted. Any ideas? Thank you for help!


Here is the function and constants:

e2 = 5
e3 = 1
omega1 = 1.1242*10^(16)
omega2 = 4.5896*10^(15)
omega3 = 5.5927*10^(15)
e1 = 1 - omega1^2/x^2
c = 3*10^8
k1 = Sqrt[y^2 - x^2/c^2 e1]
k2 = Sqrt[y^2 - x^2/c^2 e2]
k3 = Sqrt[y^2 - x^2/c^2 e3]
h = 5*10^(-9)
y = a + I b
F = ((k2/e2 + k1/e1)*(k3/e3 + k2/e2))/((k1/e1 - k2/e2)*(k3/e3 - 
       k2/e2)) - Exp[-2 k2 h]

Here is my attemption to plot roots I that I tried with FindRoot. But I am not sure it plots all roots correctly and I don't even know where to search for them:

 findPoint[i_] := 
 Evaluate[Function[
   Assuming[a \[Element] Reals && b \[Element] Reals, 
    FindRoot[{(Re[F] /. x -> (10000*i)) == 
       0, (Im[F] /. x -> (10000*i)) == 0}, {{a, 0.5}, {b, 0.2}}]]]]
list1 := Table[b /. findPoint[i][[1]], {i, 1, 2000}]
list2 := Table[10000*i, {i, 1, 2000}]
PointDone := Inner[List, list1, list2, List]
Show[{Graphics[Point[PointDone], Axes -> False, AspectRatio -> 1/2, 
   Frame -> {{True, False}, {True, False}}, FrameTicksStyle -> 18]}]

Here is a plot I got..

enter image description here

share|improve this question
    
I'm so sorry! I confused letters and there was a mistake. I edited post, and now it seems to be correct. Thank you! –  Gretchen Apr 5 '13 at 16:32
1  
Please show your FindRoot usage –  belisarius Apr 5 '13 at 16:43
    
Now showed, @belisarius :) –  Gretchen Apr 6 '13 at 7:11
    
You pasted a piece of code without checking if it really runs. –  belisarius Apr 6 '13 at 10:10
    
@belisarius, I'm sorry again. I used small f notation, and forgot that in post I used big 'F'. Now everything should work. –  Gretchen Apr 6 '13 at 16:14
show 1 more comment

1 Answer 1

up vote 5 down vote accepted

If you provide a value for a you will get a plot:

Block[{a = 10^7}, ContourPlot[{Re[F] == 0, Im[F] == 0}, {b, 0.1, 10^8}, {x, 0.1, 10^16}]]

enter image description here

You could also do a 3D contour plot including a as one of the axes:

ContourPlot3D[{Re[F] == 0, Im[F] == 0}, {a, 0, 10^8}, {b, 0.1, 10^8}, {x, 0.1, 10^16}, 
 ContourStyle -> {Opacity[0.5, Red], Opacity[0.5, Green]}]

enter image description here

Update

Here's one way to get a plot of $b$ versus $x$. Looking at the first ContourPlot we can see that there's a root near $y = 10^7 + 3 \times 10^7 i$ for $x = 8 \times 10^{15}$. Using that root as a starting point, we can slowly increase (or decrease) $x$, at each step using the previous root as the initial estimate in FindRoot.

First restart the kernel and run the initialization code again, but without setting y = a + I b. The following code then creates a list of {x,y} values where y is the complex root of F[x,y].

root[x1_, rootestimate_] := {x1, y /. FindRoot[F /. x -> x1, {y, rootestimate}]}

start = root[8*10^15, 10^7 + 3*10^7 I];

dataleft = With[{xstep = -10^14}, NestList[root[#[[1]] + xstep, #[[2]]] &, start, 50]];   
dataright = With[{xstep = +10^14}, NestList[root[#[[1]] + xstep, #[[2]]] &, start, 50]];   
alldata = Sort@Join[dataleft, dataright];

The plot of Im[y] against x is then:

ListLinePlot[{#1, Im[#2]} & @@@ alldata]

enter image description here

share|improve this answer
    
Thank you very much for answer! But when I fix a I get not whole b[x] but only one point in 2D case(the intersection of curves). The 3D case is nothing like I need to obtain. –  Gretchen Apr 6 '13 at 6:51
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.