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 got stuck on FindRoot and I didn't see any similar problem posted, so let me explain what I am trying to do and what problem I meet here.

I try to find roots of a particular function, which in the lowest order is just a second-order polynomial and hence can be solved by hand easily. However I want to generalize my code to compute situations of higher orders, in which will be transcendental equations, so it is still necessary to see things work correctly from the very beginning.

Here is the code and the output in the lowest order:

a = FindRoot[(W -100 + I/2 (1 + 1/10))^2 + 1/4 == 0, {W, i I}, AccuracyGoal -> Infinity, PrecisionGoal -> MachinePrecision, WorkingPrecision -> MachinePrecision] // Hold
c = Table[{Re[(a // ReleaseHold)[[1, 2]]-100], Im[(a // ReleaseHold)[[1, 2]]]}, {i, -2, 0, 0.001}];
c = DeleteDuplicates[c, (Abs[#1[[2]] - #2[[2]]] < 0.1 &)];
ListPlot[c, PlotRange -> {{-4.2, 4.2}, {-1.5, 0.5}}, PlotMarkers -> Automatic]

poles of case 1

so basically I create a table of different initial values, each one differs by 0.001, and use them to "scan" the position of roots on the complex plane by putting them into FindRoot one by one. Finally I get rid of duplicates by certain criteria and do a ListPlot. The result I expect is simply two complex numbers: 100-i/20 and 100-21i/20.

However, no matter how I rewrite my code I always get three solutions (the middle one in the result of ListPlot, in which I shift the origin by 100). Only when I enlarge the grid size from 0.001 to 0.1 can I kill the extra solution which is obvious wrong. But this is really bothering me because when I go to higher orders, I expect to have more points in the region I plot. By that time I would certainly need a smaller grid size in order to do the scanning.

For example, the next order I want to consider is the following:

a = FindRoot[(W -100 + I/2 (1 + 1/10))^2 + Exp[I W Pi/100]/4 == 0, {W, i}, AccuracyGoal -> Infinity, PrecisionGoal -> MachinePrecision, WorkingPrecision -> MachinePrecision] //Hold
c = Table[{Re[(a // ReleaseHold)[[1, 2]]-100], Im[(a // ReleaseHold)[[1, 2]]]}, {i, 100-4, 100+4, 10^-3}];
c = DeleteDuplicates[c, (Abs[#1[[2]] - #2[[2]]] < 0.1 &) || Abs[#1[[1]] - #2[[1]]] < 10^-9 &]
ListPlot[c, PlotRange -> {{-4.2, 4.2}, {-1.5, 0.5}}, PlotMarkers -> Automatic]

poles of case 2

This time there will be infinite number of roots since now it's transcendental, but in the region I plot I still expect only two roots for physical reasons. However once again you can see that there is a third one right on the imaginary axis which shouldn't be there. (One can easily show that no root can locate on the imaginary axis by setting W=100+yi and separating the real part and imaginary part, but let me omit the analysis here.)

Does anyone ever meet a similar situation? I really have no idea what's wrong here...I appreciate any suggestion/hint/answer!

ps. this is my first time posting a question here, please let me know how I can improve my question to be more well-organized or well-defined. Thanks!

share|improve this question
2  
FindRoot is not guaranteed to find a valid solution, but it will always return a number. The numerical method may stop before it has found a solution and will will warn you that the solution is not correct. –  Szabolcs Apr 30 '13 at 16:41
2  
Unrelated: you can improve your code by using a:=... (not a=...) and getting rid of the Hold/ReleaseHold. –  Szabolcs Apr 30 '13 at 16:43
1  
You might be interested in FindAllCrossings2D[] to assist you in searching for roots in the complex plane. –  J. M. Apr 30 '13 at 16:48
1  
@Leo You need to understand the numerical method that's being used. I believe it uses Newton's method, which is a very simple method (look it up). If you give a starting point with -0.55 I as the imaginary part, the imaginary part will be stable under the Newton-type iteration: if you substitute -0.55 I into your equation, you'll see that the LHS will have 0 as its imaginary part. While this is an unstable fixed point, the numerical errors seem not to be big enough to let the method leave it, so it's stuck there forever. –  Szabolcs Apr 30 '13 at 16:58
3  
You can check the result by backsubstituting it into the equation, but beware of numerical errors! The warning you get is not a guarantee, just a warning that something may be wrong and you should investigate further. Generally, with numerical methods like this often there are no guarantees about the correctness or completeness of the solution. But there are ways to diagnose the solution and get hints about what's going on. –  Szabolcs Apr 30 '13 at 17:00

Your Answer

 
discard

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

Browse other questions tagged or ask your own question.