2
$\begingroup$

The Left hand side of the equation has only real part. But the right hand side has real and imaginary parts. So i am using a contour plot. The intersection of lefthand side and righthand side of plot gives my solution. But can anyone tell me how to find the intersection points? I tried using Findroot. But that gives me an error message "The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient \ decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. Can anyone help me with this problem? The code is as follows

 d1[x_, r_] := 2^(1/3)/m^(1/3)*AiryAi[-2^((1/3))/m^(1/3)*((a*x*r) - m)];

    e1[x_, r_] := 
      2^(2/3)/m*AiryAiPrime[-(2^((1/3))/m^((1/3)))*((a*x*r) - m)];
    f1[x_, r_] := d1[x, r] + e1[x, r];
    g1[x_, r_] := 2^(1/3)/m^(1/3)*AiryAi[-2^((1/3))/m^(1/3)*((b*x*r) - m)];
    h1[x_, r_] := 
      2^(2/3)/m*AiryAiPrime[-(2^((1/3))/m^((1/3)))*((b*x*r) - m)];
    i1[x_, r_] := g1[x, r] + h1[x, r];
    j1[x_, r_] := -2^((1/3))/m^(1/3)*
       AiryBi[-2^((1/3))/m^(1/3)*((b*x*r) - m)];
    k1[x_, r_] := -2^((2/3))/m*
       AiryBiPrime[-(2^((1/3))/m^((1/3)))*((b*x*r) - m)];
    l1[x_, r_] := j1[x, r] + k1[x, r];
l2[x_, r_] := I*l1[x, r];
p1[x_, r_] := i1[x, r] + l2[x, r];
p2[x_] := p1[x, r] /. r -> 175*10^(-6);
f2[x_] := f1[x, r] /. r -> 175*10^(-6);
q1[x_, r_] := D[f1[x, r], r];
q2[x_] := q1[x, r] /. r -> 175*10^(-6);
q3[x_, r_] := D[p1[x, r], r];
q5[x_] := q3[x, r] /. r -> 175*10^(-6);

    FindRoot[q2[x]/f2[x] == q5[x]/p2[x], {x, 8.9*10^(14)}, 
     AccuracyGoal -> Infinity]
    ContourPlot[
     q2[x]/f2[x] == q5[x]/p2[x], {x, 1*10^(14), 2*10^(15)}, {y, 1*10^(14),
       2*10^(15)}]
Improve this question
$\endgroup$
3
  • $\begingroup$ Please show your data and code. $\endgroup$
    – David G. Stork
    Mar 29, 2017 at 20:48
  • $\begingroup$ whats m,a ,b......? $\endgroup$
    – george2079
    Mar 29, 2017 at 21:19
  • $\begingroup$ those are just numbers m=1000, a = 4.8*10^(-9),b =3.33*10^(-9) $\endgroup$
    – Raj
    Mar 29, 2017 at 21:24

1 Answer 1

Reset to default
1
$\begingroup$

Edit: Simplified Code

Starting from the FindRoot argument, but with the constants rationalized for improved accuracy,

fn = Simplify[q2[x]/f2[x] - q5[x]/p2[x] /. {m -> 1000, a -> 48 10^-10, 
     b -> 333 10^-11}];

FindRoot returns the warning given in the question, along with the correct answer.

FindRoot[fn, {x, 8.9*10^14}, AccuracyGoal -> Infinity]
(* FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than MachinePrecision digits of working precision to meet these tolerances. *)
(* {x -> 1.39562*10^15 - 2.92666*10^-260 I} *)

This is true whether or not the option AccuracyGoal -> Infinity is used. In cases like this, following the advice of the warning often is a good idea. To do so, use a larger WorkingPrecision.

FindRoot[fn, {x, 8.9*10^14}, WorkingPrecision -> 30]
(* {x -> 1.39561905279431854837087407070*10^15 - 
         2.92665880860180596403069585524*10^-260 I} *)

but with no warning.

Addendum: Multiple Roots

To obtain all roots in a given range, it is helpful first to plot the function itself.

Themes`AddThemeRules["mystyle", AxesStyle -> Directive[Black, Bold, Medium], 
    AspectRatio -> 1/4, ImageSize -> 1000]; $PlotTheme = "mystyle";
Plot[Evaluate@ReIm@fn, {x, 0, 1.3 10^15}, PlotPoints -> 100]
Plot[Evaluate@ReIm@fn, {x, 1.1 10^15, 1.9 10^15}, PlotPoints -> 1000]

enter image description here

enter image description here

The function becomes oscillatory at x near 1.2 10^10, corresponding approximately to x == m/(a r), and becomes noticeably complex at x near 1.7 10^10, corresponding approximately to x == m/(b r).

Because the roots are more or less uniformly spaced and interleaved between poles, it is straightforward and efficient to compute them using FindRoot with numerous uniformly spaced initial guesses and Union to eliminate duplicates.

rts = Union[Table[(x /. FindRoot[fn, {x, x0}, WorkingPrecision -> 30]), 
    {x0, 12 10^14, 20 10^14, 2 10^12}], SameTest -> (Norm[#1 - #2] < 1 &)];
ListPlot[ReIm@rts, AxesLabel -> {Re[x], Im[x]}, PlotStyle -> PointSize[Medium]]

yields the first 169 roots in several seconds. Their real and imaginary parts are displayed below.

enter image description here

(Slight irregularities in the three figures are plotting artifacts.) Some alternative, but here less efficient, root-finding techniques are given in question 5663

Improve this answer
$\endgroup$
9
  • $\begingroup$ it worked..But when i change the starting point for example to 8.5*10^14, I still get an answer which is not right. But any starting point after 8.9*10^(14), the solution makes sense to my problem. I was wondering if it is possible to find all the possible roots in a range. $\endgroup$
    – Raj
    Mar 31, 2017 at 19:23
  • $\begingroup$ @sreekul Actually, x -> 0 is a valid answer, just not the one you want. To obtain more answers at large x, choose initial guesses greater than about 10^15. However, finding all the answer is a different matter. I believe that there are very many of them. $\endgroup$
    – bbgodfrey
    Mar 31, 2017 at 19:50
  • $\begingroup$ Is there any way to find all the possible solutions in a range. Apart from FindRoot does mathematica have any way to find solutions in a range. $\endgroup$
    – Raj
    Mar 31, 2017 at 20:59
  • $\begingroup$ "x" in general is a complex number and i am looking for the general solution. $\endgroup$
    – Raj
    Mar 31, 2017 at 21:50
  • $\begingroup$ Reduce[{Re[fnr] == 0, 1.2 10^15 <= x <= 1.4 10^15}, x] is giving me an error message "Reduce was unable to solve the system with inexact coefficients or \ the system obtained by direct rationalization of inexact numbers \ present in the system. Since many of the methods used by Reduce \ require exact input, providing Reduce with an exact version of the \ system may help. \!(*ButtonBox[\">>\", Appearance->{Automatic, None}, BaseStyle->\"Link\", ButtonData:>\"paclet:ref/Reduce\", ButtonNote->\"Reduce::inex\"])" $\endgroup$
    – Raj
    Mar 31, 2017 at 22:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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