3
$\begingroup$

I have an expression which is a combination of Bessel functions:

 eq[k_, ω_] = 0.004 k Sqrt[-k^2 + ω^2] (k^2 + ω^2 + 
k ω (-2 + 1.3*10^(-14) ω^2)) BesselJ[1, √((0.002 k - 
    0.002 ω + ω (-0.9 k + 0.0003/(
       k - ω) + ω) (-k^2 + ω^2) (-1 + 
       0.0004/(-0.9 k + ω)^2))/(0.002 - ω (-k + 
       0.0003/(k - ω) + ω) (1 + 
       0.01/(-0.001 + 
        5 (-k + ω)^2) - (0.01 (0.002 k^2 (-k + 0.0003/(
              k - ω) + ω)))/((-0.001 + (-k + ω)^2) (0.002 - ω^2 (-0.9 k + 0.0003/(
               k - ω) + ω))))))] ((k^2 - 
   0.7 ω^2) BesselJ[1, 
  2 Sqrt[-k^2 + 0.5 ω^2]] BesselY[1, 
  1.5 Sqrt[-k^2 + 0.5 ω^2]] + 
0.3 k^2 BesselJ[1, 1.5 Sqrt[-k^2 + 0.5 ω^2]] BesselY[1, 
  2 Sqrt[-k^2 + 0.5 ω^2]]) (BesselJ[1, 
  1.5 Sqrt[-k^2 + ω^2]] BesselY[1, 
  Sqrt[-k^2 + ω^2]] - 
BesselJ[1, Sqrt[-k^2 + ω^2]] BesselY[1, 
  1.5 Sqrt[-k^2 + ω^2]])

Form of equation

I want to find root of this equation for a specified value of k=1. By plotting this equation I realized that this equation acquire real and imaginary values so I tried the FindRoot& NSove & Reducecommands while restricting the range of /[Omega] as follow:

Reduce[eq[1,/[Omega]==0 && -1<Re[ω]<1 && -1<Im[ω]<1,ω]

But unfortunately it doesn't work at all. Can any one resolve this problem for me? Thanks in advance.

$\endgroup$
3
  • 1
    $\begingroup$ For me, FindRoot[eq[1, w], {w, 0}] works, and returns $$w=0.00222463$$ with some small imaginary value $\endgroup$
    – yohbs
    Commented Jun 18, 2015 at 14:29
  • 1
    $\begingroup$ Why do you mean by "it doesn't work at all". You get no answer, an error, it churns indefinitely...? Your Reduce expression contains syntax errors. Generally, I would suspect that your expression is too complicated to attempt to get a general solution with Solve and Reduce. NSolve mostly deals with polynomial equations anyway, so it won't help much here. You are left with FindRoot, which works just fine in your case, as @yohbs showed above, but won't give you all solutions. $\endgroup$
    – MarcoB
    Commented Jun 18, 2015 at 14:36
  • $\begingroup$ I know FindRoot command (with the above mentioned format) is useful to seek for a single root but I can't find all roots. $\endgroup$ Commented Jun 18, 2015 at 14:47

1 Answer 1

7
$\begingroup$

Using your definition of eq, if you try Solve even without imposing conditions on ω, you obtain the following error:

Solve[eq[1, ω] == 0, ω]

Solve::inex: Solve 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 Solve require exact input, providing Solve with an exact version of the system may help.

This is very valuable information! Mathematica is really giving you a big pointer in the right direction here.

Typically Solve will attempt to Rationalize a system with inexact coefficients, in hopes of obtaining something it can deal with more easily. However, in your case this rationalization step seems to have failed. Upon inspection of your coefficients, this is because of the value of one particularly small coefficient in your expression. In fact, see what happens if we attempt the rationalization by hand:

Rationalize[eq[1, ω]]

(1/250)*Sqrt[-1 + ω^2]*(1 + ω^2 + ω*(-2 + 1.3*^-14*ω^2))*BesselJ[1, Sqrt[(1/500 - ω/500 + (-1 + 1/(2500*(-(9/10) + ω)^2))*ω*(-(9/10) + 3/(10000*(1 - ω)) + ω)*(-1 + ω^2))/(1/500 - ω*(-1 + 3/(10000*(1 - ω)) + ω)*(1 + 1/(100*(-(1/1000) + 5*(-1 + ω)^2)) - (-1 + 3/(10000*(1 - ω)) + ω)/(50000*(-(1/1000) + (-1 + ω)^2)*(1/500 - ω^2*(-(9/10) + 3/(10000*(1 - ω)) + ω)))))]]*((1 - (7*ω^2)/10)*BesselJ[1, 2*Sqrt[-1 + ω^2/2]]*BesselY[1, (3/2)*Sqrt[-1 + ω^2/2]] + (3/10)*BesselJ[1, (3/2)*Sqrt[-1 + ω^2/2]]*BesselY[1, 2*Sqrt[-1 + ω^2/2]])*(BesselJ[1, (3/2)*Sqrt[-1 + ω^2]]*BesselY[1, Sqrt[-1 + ω^2]] - BesselJ[1, Sqrt[-1 + ω^2]]*BesselY[1, (3/2)*Sqrt[-1 + ω^2]])

Almost everything was transformed into a rational number, except the very small $1.3\times 10^{-14}$ value, which is too small for Rationalize's default tolerance (see its documentation page, under "Details").

However, we can force every number to a rational representation by setting this tolerance to $0$:

Rationalize[eq[1, ω], 0]

(1/250)*Sqrt[-1 + ω^2]*(1 + ω^2 + ω*(-2 + (13*ω^2)/1000000000000000))*BesselJ[1, Sqrt[(1/500 - ω/500 + (-1 + 1/(2500*(-(9/10) + ω)^2))*ω*(-(9/10) + 3/(10000*(1 - ω)) + ω)*(-1 + ω^2))/(1/500 - ω*(-1 + 3/(10000*(1 - ω)) + ω)*(1 + 1/(100*(-(1/1000) + 5*(-1 + ω)^2)) - (-1 + 3/(10000*(1 - ω)) + ω)/(50000*(-(1/1000) + (-1 + ω)^2)*(1/500 - ω^2*(-(9/10) + 3/(10000*(1 - ω)) + ω)))))]]*((1 - (7*ω^2)/10)*BesselJ[1, 2*Sqrt[-1 + ω^2/2]]*BesselY[1, (3/2)*Sqrt[-1 + ω^2/2]] + (3/10)*BesselJ[1, (3/2)*Sqrt[-1 + ω^2/2]]*BesselY[1, 2*Sqrt[-1 + ω^2/2]])*(BesselJ[1, (3/2)*Sqrt[-1 + ω^2]]*BesselY[1, Sqrt[-1 + ω^2]] - BesselJ[1, Sqrt[-1 + ω^2]]*BesselY[1, (3/2)*Sqrt[-1 + ω^2]])

This should help Solve for your problem:

Solve[Rationalize[eq[1, ω], 0] == 0, ω, Cubics -> False, Quartics -> False]

Not specifying a domain implicitly allows for complex solutions; this quickly evaluates to give nine solutions, expressed as Root objects. A numerical estimate of those solutions can be obtained with:

numsolns = Chop[N[solns]]

(* Out:
{ {ω -> -7.69231*10^13}, {ω -> 1. - 1.14018*10^-7 I}, 
  {ω -> 1. + 1.14018*10^-7 I}, {ω -> -1.00105}, {ω -> 0.00222463}, 
  {ω -> 0.8793}, {ω -> 0.897109}, {ω -> 0.91953}, {ω -> 1.00289} }
*)

Of course you can substitute those back into eq[1, ω] to check that they are actual roots for that expression:

Chop[eq[1, ω] /. numsolns]

(*Out: {0, 0, 0, 0, 0, 0, 0, 0, 0} *)

Finally, you can select those solutions whose real and imaginary parts are in the interval $[-1,1]$ using e.g. Select:

Select[ω /. numsolns, -1 < Re[#] < 1 && -1 < Im[#] < 1 &]

(* Out:
{1. - 1.14018*10^-7 I, 1. + 1.14018*10^-7 I, 0.00222463, 0.8793, 0.897109, 0.91953}
*)

A different take:

That very small coefficient is odd in comparison to all others in the OP's expression. I had wondered if it was just a leftover from roundoff errors. @Guesswhoitis brought up the very same point, so it is worth exploring.

We can Chop that tiny value off, Rationalize the resulting expression, and feed it to Solve:

choppedeq[k_, ω_] = Chop[eq[k, ω]];

choppedsolns = Solve[Rationalize[choppedeq[1, ω]] == 0, ω];
numchoppedsolns = ω /. N@choppedsolns

(* Out: {-1.00105, 0.00222463, 0.8793, 0.897109, 0.91953, 1.00289} *)

In this case Solve returns six Root objects, whose approximate numerical values are shown above. The solutions are all real! One wonders about the original problem...

We can also take a look at the expression itself after removing that small value:

GraphicsGrid@ With[
  {opts = Sequence[Evaluated -> True, PlotRangePadding -> Scaled[0.05], Epilog -> {AbsolutePointSize[10], Red, Point[{#, 0} & /@ numchoppedsolns]}]},
  {
   {Plot[Through[{Abs, Re, Im}[choppedeq[1, ω]]], {ω, -1.1, 1.1}, opts, AspectRatio -> Full, PlotRange -> 0.002], SpanFromLeft, SpanFromLeft, SpanFromLeft},
    Plot[Through[{Abs, Re, Im}[choppedeq[1, ω]]], {ω, #1, #2}, opts, ImageSize -> 250] & @@@ {{-1.1, -.9}, {0, 0.005}, {0.87, 0.93}, {0.98, 1.02}}
  }
]

grid of solutions

$\endgroup$
14
  • 2
    $\begingroup$ Sometimes, one just has to read error messages carefully instead of tearing your hair out after seeing them… :) But, I now wonder whether that tiny constant is actually a characteristic of the system, or mere rounding error. $\endgroup$ Commented Jun 18, 2015 at 21:25
  • 2
    $\begingroup$ @Guesswhoitis. I wondered the same thing. I had chopped that tiny constant off just for fun, and the results are wildly different, but in the interest of answering the question as faithfully as possible, I left those off. I'll dig them back up and include them in the answer. $\endgroup$
    – MarcoB
    Commented Jun 19, 2015 at 4:28
  • 1
    $\begingroup$ "the results are wildly different" - a huge red flag, in my book. It can be that OP's actual problem is badly scaled, or worse, ill-conditioned. Wilkinson's classic example comes to mind… $\endgroup$ Commented Jun 19, 2015 at 4:44
  • 1
    $\begingroup$ Here's another experiment you can try: what if you replace Rationalize[#, 0] & with SetPrecision[#, ∞] &? $\endgroup$ Commented Jun 19, 2015 at 10:14
  • 1
    $\begingroup$ @Guesswhoitis. I ran that just now, and there is no difference between the solutions found for the Rationalize[#, 0]& expression and the SetPrecision[#, Infinity]&. What was your reasoning in these cases though? Numerical instability? I would have thought that something would have been quite wrong if they had been different. $\endgroup$
    – MarcoB
    Commented Jun 19, 2015 at 11:14

Your Answer

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

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