3

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!

Oleksandr R.
  • 23,023
  • 4
  • 87
  • 125
Leo Fang
  • 754
  • 7
  • 17
  • 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.'s missing motivation Apr 30 '13 at 16:48
  • @Szabolcs, yeah I saw the message warning me that it had hit the requested precision/accuracy goal but couldn't find a sufficient decrease in the merit function. Is this a warning that the solution is definitely incorrect?

    And thanks, I know I could define a:=..., I just forgot to refine this part of code. :)

     – Leo Fang Apr 30 '13 at 16:54
  • @J.M., thanks a lot for the link! Let me look at it and see if I can implement it! – Leo Fang Apr 30 '13 at 16:56
  • 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
  • @Szabolcs, wow that's a very helpful observation. Thank you so much! I did know I was using Newton's method here but I didn't think in this way...But if I do a similar analysis to the second case, it seems not working because in this case I scan the roots along the real axis and far away from those shown roots. Do I miss something? – Leo Fang Apr 30 '13 at 17:02
  • by the way, found another related thread seems also useful on this problem!

    http://mathematica.stackexchange.com/questions/5663/about-multi-root-search-in-mathematica-for-transcendental-equations

     – Leo Fang Apr 30 '13 at 22:10
  • @Szabolcs, I plotted both real part as well as imaginary part of the second case, and now I totally understand what you were saying! Thanks so much! – Leo Fang May 01 '13 at 00:43
  • @J.M., thanks for letting me know the FindAllCrossings2D[] function! It is amazing! – Leo Fang May 02 '13 at 15:21
  • You're welcome, Leo. Thank Stan Wagon for coming up with it. – J. M.'s missing motivation May 02 '13 at 15:22
  • Sorry guys, I found myself puzzled by Wagon's solution...I don't have enough reputation to leave a comment here, so forgive me to ask in this post.

    I don't understand how the line seeds = Flatten[(...), 1]; works. It seems that it generates a list of seeds by finding crossings of {f==0,g==0} and feed into FindRoots to get more accurate numbers. But can somebody elaborate how exactly it works? To me, for example, the line signs*RotateRight[signs] is rather confusing...

     – Leo Fang May 04 '13 at 19:31

0 Answers0