# Solving equation with Bessel functions BesselJ & BesselK

Posted 9 years ago
9661 Views
|
13 Replies
|
5 Total Likes
|
 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
13 Replies
Sort By:
Posted 9 years ago
 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 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 9 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 9 years ago
 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 9 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 9 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 9 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 9 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 9 years ago
 Thank you. But, it seems that your method fails sometimes, if the function is sophisticated enough.
Posted 9 years ago
 EVERYTHING fails sometimes. The question is whether it is useful or not.
Posted 9 years ago
 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 9 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 9 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]]