Message Boards Message Boards

0
|
10935 Views
|
13 Replies
|
5 Total Likes
View groups...
Share
Share this post:

Solving equation with Bessel functions BesselJ & BesselK

Posted 10 years ago

Hi, I am trying to solve the following equation using Mathematica. NSolve[] seems not to be capable of finding the solution. NSolve[] returns back the same equation after "thinking" for a while. Solve[] does nearly the same thing. Is there anything, I can change in the equation to get a better result.

a = 25 10^-6; n1 = 1.458 + 5 10^-3; n2 = 1.458 ; [Lambda] = 
 1.5 10^-6; k0 = (2 [Pi])/[Lambda]; p = Sqrt[n1^2 k0^2 - [Beta]^2];
 q = Sqrt[[Beta]^2 - n2^2 k0^2];
DJm[z_] := !(
*SubscriptBox[([PartialD]), (z)](BesselJ[m, z])) /. m -> 1;
DKm[z_] := !(
*SubscriptBox[([PartialD]), (z)](BesselK[m, z])) /. m -> 1;




NSolve[((1/2 (BesselJ[0, p a] - BesselJ[2, p a]))/(
     p BesselJ[1, p a]) + (1/2 (-BesselK[0, q a] - BesselK[2, q a]))/(
     q BesselK[1, q a])) ((1/2 (BesselJ[0, p a] - BesselJ[2, p a]))/(
     p BesselJ[1, p a]) + 
     n2^2/n1^2 (1/2 (-BesselK[0, q a] - BesselK[2, q a]))/(
      q BesselK[1, q a])) == (1)^2/
   a^2 (1/p^2 + 1/q^2) (1/p^2 + n2^2/n1^2 1/q^2), [Beta]]

Thanks in advance

POSTED BY: Saf Al
13 Replies

Note that:

Reduce[
 n1^2 k0^2 - [Beta]^2 > 0 && [Beta]^2 - n2^2 k0^2 > 0]

gives

 -6.1282  10^6 < [Beta] < -6.10726  10^6 || 
 6.10726 10^6 < [Beta] < 6.1282  10^6

so any places where p and q are both real are located in one of these ranges.

So go ahead and plot the full expression (i.e., the left hand side minus the right hand side) in, for example, the second range:

Plot[((1/2 (BesselJ[0, p a] - BesselJ[2, p a]))/(p BesselJ[1, 
         p a]) + (1/
         2 (-BesselK[0, q a] - BesselK[2, q a]))/(q BesselK[1, 
         q a])) ((1/2 (BesselJ[0, p a] - BesselJ[2, p a]))/(p BesselJ[
         1, p a]) + 
     n2^2/n1^2 (1/2 (-BesselK[0, q a] - BesselK[2, q a]))/(q BesselK[
          1, q a])) - (1)^2/
    a^2 (1/p^2 + 1/q^2) (1/p^2 + n2^2/n1^2 1/q^2), {[Beta], 
  6.107256118578557`*^6, 6.1282000696024895`*^6}]

giving

enter image description here

So there appear to be some zeros in this region. And one can use this plot to select a starting point for the FindRoot function as in:

FindRoot[((1/2 (BesselJ[0, p a] - BesselJ[2, p a]))/(p BesselJ[1, 
          p a]) + (1/
          2 (-BesselK[0, q a] - BesselK[2, q a]))/(q BesselK[1, 
          q a])) ((1/
          2 (BesselJ[0, p a] - BesselJ[2, p a]))/(p BesselJ[1, p a]) +
       n2^2/n1^2 (1/
           2 (-BesselK[0, q a] - BesselK[2, q a]))/(q BesselK[1, 
           q a])) - (1)^2/
     a^2 (1/p^2 + 1/q^2) (1/p^2 + n2^2/n1^2 1/q^2) == 0, {[Beta], 
  6.12 10^6}]

giving

{[Beta] -> 6.11989[CenterDot]10^6}

But there are several solutions around here (to the accuracy of the FindRoot computation), some possibly elsewhere on the complex plane.

using 6.125 10^6 as the seed value gives a root at

{[Beta] -> 6.1248[CenterDot]10^6}
POSTED BY: David Reiss
Posted 10 years ago

Perfect! I am trying to calculate the propagation constant (Beta) for modes in an optical fiber according to a well known equation. The method, you used, works just fine. However, you assumed that p and q are both real. I was wondering what led to that assumption. I checked the reference, I am reading for the time being. It does not say anything about p and q being purely real. Of course, if they are complex, then I can't use the REDUCE[] command to find the areas, where the solutions exist. Any advice?

Thanks in advance

POSTED BY: Saf Al

Actually I didn't assume that p and q were real. What I did was ask what happens int he region where they are real as a way of tracking down a place where I could potentially plot the function and get real values. Given that all the other parameters were real, I then suspected that the value of the function would be real in those areas. The hunch was correct in that there was a plot generated and that indeed the plot seemed to cross the axis--indicating some zeros in that region. There may well be other zeros (or curves of zeros) elsewhere in the complex plane--I haven;t thought much more about it. But the key lesson is to explore the behavior of the function and learn as much about it as possible. Then one goes ahead and uses appropriate Mathematica functions to try to extract the locations of zeros. Some functions and algorithms can be automatically plugged into a high level Mathematica function which then goes on to find solutions on its own. In messy (i.e., many real world) cases though one needs to use Mathematica as a hands-on tool to explore the behavior of the objects that one wants to extract information from (e.g., like zeros). So I think that the best way to take what I wrote was as a bit of an example of that sort of exploratory play.

POSTED BY: David Reiss
Posted 10 years ago

Yet, your solution is neat. Which makes me wonder whether q and p have, in fact, real values. I guess, this will be my next exploration step. But regardless of this, I will keep this way of function-analyzing in mind whenever Solve[] fails. Thanks. you've been very helpful.

POSTED BY: Saf Al
Posted 10 years ago

For any WorkingPrecision between 48 and 128 NMinimize converges on a beta about 153746 where the resulting minimum is significantly different from zero and does not get closer to zero with increasing Working precision.

In[7]:= a = 25 10^-6; n1 = 1458/1000 + 5 10^-3; n2 = 1458/1000; \[Lambda] = 15/10 10^-6;
k0 = (2 \[Pi])/\[Lambda]; p = Sqrt[n1^2 k0^2 - \[Beta]^2];q = Sqrt[\[Beta]^2 - n2^2 k0^2];
NMinimize[ Norm[((1/2 (BesselJ[0, p a] - BesselJ[2, p a]))/(p BesselJ[1, p a]) +
 (1/2 (-BesselK[0, q a] - BesselK[2, q a]))/(q BesselK[1, q a])) ((1/2 (BesselJ[0, p a] -
 BesselJ[2, p a]))/(p BesselJ[1, p a]) +n2^2/n1^2 (1/2 (-BesselK[0, q a] - BesselK[2, q a]))/
 (q BesselK[1, q a])) - (1)^2/a^2 (1/p^2 + 1/q^2) (1/p^2 + n2^2/n1^2 1/q^2)], \[Beta], 
 WorkingPrecision -> 128]

Out[9]= {2.6644908885382207535682631873522244614076638003775608946182982\
 185058485936623112823132496185424062412703162736336930371650049643*10^-14,
 {\[Beta] -> 153746.742756413486195499781367574310403493602174002003504819\
 78661459273051117012575097394837519002121150764177004071175454428111}}

Because of the way that minimum varies with WorkingPrecision I'm a little skeptical whether that is precisely correct or not.

POSTED BY: Bill Simpson
Posted 10 years ago

Thanks for the reply. I am curious about the how using the Norm[] command would help finding the roots of the mentioned equation. It would be kind if you'd explain your method briefly. Regards.

POSTED BY: Saf Al
Posted 10 years ago

Sometimes Solve, NSolve, Reduce, FindRoot, etc, etc can have great difficulty and take a long time looking for a zero of a complicated function. NMinimize can sometimes be faster, but you want it to find where your complicated function is near zero, whether it is positive, negative or even Complex, Norm will give a measure of the distance from zero for all those cases. If it finds the minimum of the norm is nonzero then that may indicate there is no zero but, as was shown, it is sometimes possible that it just did not search far enough to find a root.

POSTED BY: Bill Simpson
Posted 10 years ago

Thank you. But, it seems that your method fails sometimes, if the function is sophisticated enough.

POSTED BY: Saf Al
Posted 10 years ago

EVERYTHING fails sometimes. The question is whether it is useful or not.

POSTED BY: Bill Simpson

If I plot the left-hand side of the equation minus the right-hand side, as a function of beta, I basically get 0 for any value of beta

POSTED BY: Frank Kampas
Posted 10 years ago

Thanks for the reply. The difficulty with Plotting is probably a result of P or/and q being complex depending on the value of (Beta). Regards.

POSTED BY: Saf Al

Your code is corrupted. Please edit your post to fix it or attach a notebook.

POSTED BY: Sam Carrettie
Posted 10 years ago

Hi, I don't know why the code is corrupted. I copied the code as INLINE from Mathematica notebook then I paste it in the Sample Code field using "CTR k". However, I included another code snippet in this reply. I double checked it to make sure that it is not corrupted. Thanks

a = 25 10^-6; n1 = 1.458 + 5 10^-3; n2 = 1.458 ; \[Lambda] = 
 1.5 10^-6; k0 = (2 \[Pi])/\[Lambda]; p = Sqrt[n1^2 k0^2 - \[Beta]^2];
 q = Sqrt[\[Beta]^2 - n2^2 k0^2];


NSolve[((1/2 (BesselJ[0, p a] - BesselJ[2, p a]))/(
     p BesselJ[1, p a]) + (1/2 (-BesselK[0, q a] - BesselK[2, q a]))/(
     q BesselK[1, q a])) ((1/2 (BesselJ[0, p a] - BesselJ[2, p a]))/(
     p BesselJ[1, p a]) + 
     n2^2/n1^2 (1/2 (-BesselK[0, q a] - BesselK[2, q a]))/(
      q BesselK[1, q a])) == (1)^2/
   a^2 (1/p^2 + 1/q^2) (1/p^2 + n2^2/n1^2 1/q^2), \[Beta]]
POSTED BY: Saf Al
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract