0
$\begingroup$

I have a complicated function that I would like to solve for the roots. The function comprises Bessel and Hankel functions:

freq = 5*^9;
Omg = 2*Pi*freq;
rho = 15 ;
vco = 5907.5 (*core velocity*);
vcl = 5944 (*clad velocity*);
m = 0 ;
U = Sqrt[Omg^2/vco^2 - Omg^2/v];
W = Sqrt[Omg^2/vcl^2 - Omg^2/v];

NSolve[U*BesselJ[m + 1, U*rho]*HankelH1[m, W*rho] - W*HankelH1[m + 1, W*rho]*BesselJ[m, U*rho]==0, v]

However, NSolve does not evaluate. I have also tried using FindRoot but I have no idea what starting value to give. Plotting the function (with all variables equal to 1 for simplicity) is of no help for guessing starting values.

enter image description here

Improve this question
$\endgroup$
1
  • $\begingroup$ I'm getting the feeling that there are no roots. Is that possible? $\endgroup$
    – Michael E2
    Commented Nov 22, 2020 at 17:48

1 Answer 1

Reset to default
1
$\begingroup$ 
$Version

(* "12.1.1 for Mac OS X x86 (64-bit) (June 19, 2020)" *)

Clear["Global`*"]

freq = 5*^9;
Omg = 2*Pi*freq;
rho = 15;
vco = 5907.5 (*core velocity*);
vcl = 5944 (*clad velocity*);
m = 0;
U = Sqrt[Omg^2/vco^2 - Omg^2/v];
W = Sqrt[Omg^2/vcl^2 - Omg^2/v];

EDIT 2: Replacing Rationalize with Rationalize[#, 0]& to ensure conversion to rationals.

expr = U*BesselJ[m + 1, U*rho]*HankelH1[m, W*rho] - 
     W*HankelH1[m + 1, W*rho]*BesselJ[m, U*rho] // Rationalize[#, 0] & // 
   FullSimplify;

(sol = NSolve[{expr == 0, 10^7 < v < 10^8}, v, WorkingPrecision -> 50]) // N

(* {{v -> 3.533*10^7}} *)

expr /. sol // Quiet

(* {0.*10^-198964 + 0.*10^-198964 I} *)

In the definitions for U and W should the v be v^2? That provides a more reasonable value

Sqrt[v] /. sol // N

(* {5943.9} *)

EDIT: Using v^2

U = Sqrt[Omg^2/vco^2 - Omg^2/v^2];
W = Sqrt[Omg^2/vcl^2 - Omg^2/v^2];

expr = U*BesselJ[m + 1, U*rho]*HankelH1[m, W*rho] - 
     W*HankelH1[m + 1, W*rho]*BesselJ[m, U*rho] // 
    Rationalize[#, 0] & // FullSimplify;

(sol = NSolve[{expr == 0, 5000 < v < 6000}, v, 
    WorkingPrecision -> 50]) // N

(* {{v -> 5943.99}} *)
Improve this answer
$\endgroup$
4
  • $\begingroup$ Yes it should be v^2. However, I thought introducing a square root would complicate things so I just called the whole thing v. Well, two of the solution are just the parameters vco and vcl. But there is a new third solution that I didn't have before. I don't see why // Rationalize // FullSimplify would help here. It seems the real helper was the 10^7 < v < 10^8 you added. $\endgroup$
    – ngc1300
    Commented Nov 21, 2020 at 22:19
  • $\begingroup$ And I tried with v^2 and it gives different solutions, so I think using v instead really helps the numerical solver $\endgroup$
    – ngc1300
    Commented Nov 21, 2020 at 22:21
  • $\begingroup$ Also, how do you allow for complex solutions? $\endgroup$
    – ngc1300
    Commented Nov 21, 2020 at 22:35
  • $\begingroup$ Already upvoted for the method, but I think it's a harder problem than NSolve can handle. Consider Plot[RealExponent@expr, {v, 5800, 6000}, WorkingPrecision -> 100]. The gap in the graph is due to a singularity at expr /. v -> 5944. $\endgroup$
    – Michael E2
    Commented Nov 22, 2020 at 15:49

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.