TE mode trascendental equation problem
Show older comments
Hi,
I want to solve the characteristic equation describing the TE mode in an optical fiber (Light transmission optics, D. Marcuse):
where a is the core radius of the fiber. The key parameter to solve in this equation is the propagation constant β that is related to κ and γ as
Parameters 
and 
are constants in the core and cladding of the fiber related to the wavenumber 
as 
, being 
the refractive index of core and cladding, respectively.
and are constants in the core and cladding of the fiber related to the wavenumber as , being the refractive index of core and cladding, respectively.I know there are multiple solutions and are related to the zeros of the bessel function 
in the numerator. I tried to solve this equation using matlab functions solve() and fzero() for the values 
, 
, 
and 
but I'm stuck. I need some help to find the solution of this equation.
in the numerator. I tried to solve this equation using matlab functions solve() and fzero() for the values , , and but I'm stuck. I need some help to find the solution of this equation.Thanks for your attention.
Carolina Rickenstorff
8 Comments
Star Strider
on 17 Feb 2020
Post your code, and any necessary data.
Carolina Rickenstorff
on 17 Feb 2020
Carolina Rickenstorff
on 17 Feb 2020
Carolina Rickenstorff
on 17 Feb 2020
Carolina Rickenstorff
on 17 Feb 2020
Star Strider
on 17 Feb 2020
I usually use fsolve if fzero has problems:
X = fsolve(F,x0)
however even if I only use one of the ‘x0’ values, this throws:
Error using trustnleqn (line 28)
Objective function is returning undefined values at initial point. FSOLVE cannot
continue.
So your ‘x0’ is apparently not appropriate. I am not familiar with what you are doing, so I cannot help you choose more appropriate values that are consistent with what you want to solve. You need to experiment with those.
Carolina Rickenstorff
on 17 Feb 2020
It is very weird that there are books as Fiber Optics (F. Mitschke) where the same equation is written slightly different and I have no problem solving it.
where
and 
%%%%Te mode
syms u w
u=1.2;
%solution=solve((1/u).*besselj(0,u)./besselj(1,u)==(1/w).*besselk(0,w)./besselk(1,w), w);
solution=solve(u.*besselj(1,u)./besselj(0,u)==w.*besselk(1,w)./besselk(0,w), w);
solution=double(solution)
Star Strider
on 17 Feb 2020
In that situation, use the version you can solve!
Accepted Answer
David Goodmanson
on 18 Feb 2020
Edited: David Goodmanson
on 22 Feb 2020
Hi Carolina,
It's not a real good idea to go with an answer solely on the grounds that 'solve' can find one from an expression in a book. With books, there are definitely issues with the notation, who is using what, so you have to make sure you are solving the right equations. The E and H fields on each side of the boundary r = a are
E: G*J1(k1*a) = alpha*G*H1(k2*a)
H: G*(1/mu1)*k1*J0(k1*a) = alpha*G*(1/mu2)*k2*H0(k2*a)
% not all the constant conversion factors to get from E to H are shown.
% these cancel out when doing equality of H between the two regions.
where E is in the theta direction, around the cylinder, and H is in the z direction, along its length. (There is also a radial component of B that is matched automatically). G is an overall amplitude that doesn't matter for matching, alpha is an unknown constant to get the sizes to match up, and k1 and k2 are
k1 = sqrt(n1^2 k^2 - kz^2) k2 = sqrt(n2^2 k^2 - kz^2)
kz is your beta, which must be the same in both regions since the solution has to have the same z dependence down the length of the cylinder. k1 is your kappa, and k2 is almost your gamma. This is assuming that in your definition of k, k = 2pi / lambda, the wavelength is the free space version, lambda = c / f.
On the right, H1 and H0, which are really denoted H1(+) and H0(+) are certain Hanckel functions. (There is also an H1(-) and H0(-) ).
Dividing one equation by the other yields
mu1*J1(k1*a) / (k1*J0(k1*a)) = mu2*H1(k2*a) / (k2*H0(k2*a))
For a solution, kz has to be such that k1 is real and k2 is pure imaginary. Rewrite this as
k2 = i*g2 where g2 = +sqrt(kz^2 - n2^2 k^2)
and g2 is your gamma.
The fields in the outside region have to decay exponentially for large r, and the bessel functions K have that property. It turns out that
H0(i*g2*a) = -i*(2/pi)*K0(g2*a)
H1(i*g2*a) = -(2/pi)*K1(g2*a)
and after tossing those in we arrive at last at
mu1*J1(k1*a) / (k1*J0(k1*a)) = - mu2*K1(g2*a) / (g2*K0(g2*a))
Here the J's and K's are straight Matlab besselj and besselk.
Note the minus sign on the right. Various definitons of bessel functions include factors of (-)^n where n is the order of the bessel function. That messes with the minus signs.
This reply has gone on for a bit, but I beleive the process is absolutely necessary because of all the ways of defining bessel functions.
The larger the value of k*a, the more modes (roots) there will be. Your value of k*a is 40.54 which is large enough for several modes.
The code below shows where the roots are. I will leave it to you to use fzero or whatever method you choose to find roots.
a = 10e-6;
k = 2*pi/1.55e-6;
mu1 = 1;
n1 = 1.4;
mu2 = 1;
n2 = 1.33;
kzmax = n1*k;
kzmin = n2*k;
kz = linspace(kzmin,kzmax,100000);
kz(1) = [];
kz(end) = [];
k1 = sqrt(n1^2*k^2 - kz.^2);
g2 = sqrt(kz.^2 - n2^2*k^2);
C = mu1*besselj(1,k1*a)./(k1.*besselj(0,k1*a));
% quiet down the plot y values
Cm = mean(abs(C));
C(abs(C)>3*Cm) = nan;
D = -mu2*besselk(1,g2*a)./(g2.*besselk(0,g2*a));
Dm = mean(abs(D));
D(abs(D)>3*Dm) = nan;
figure(1)
plot(kz,C,kz,D)
grid on
kzroot = 5.4617301149e6; % first root from plot
% check, stay with besselh and k2, not besselk and g2
% besselh has an extra input to choose H(+) rather than H(-)
k1root = sqrt(n1^2*k^2 - kzroot^2);
k2root = sqrt(n2^2*k^2 - kzroot^2);
alpha = besselj(1,k1root*a)/besselh(1,1,k2root*a);
% E and H should agree on each side
E1 = besselj(1, k1root*a)
E2 = alpha*besselh(1,1,k2root*a)
% there are constant conversion factors to get from E to H, not shown here.
% these do not affect an equality check of H between regions 1and 2
H1 = (1/mu1)*k1root*besselj(0, k1root*a)
H2 = alpha*(1/mu2)*k2root*besselh(0,1,k2root*a)
10 Comments
Carolina Rickenstorff
on 18 Feb 2020
Carolina Rickenstorff
on 18 Feb 2020
Edited: Carolina Rickenstorff
on 18 Feb 2020
Hi Mr. Goodmanson,
Looking at your program I put directly in variable D the value
D = -1i.*besselh(1,1,1i.*g2*a)./(k2.*besselh(0,1,1i.*g2*a));
and plot it toguether with variable C. Despite D is complex when plotting the real part it produces the same roots as your original variable D defined in terms of besselk(). However when verifying the roots I mess it all. I just wonder if there is no way to avoid functions besselk() altoguether when finding the roots since the starting equation
is written in terms of besselj() and besselh(+).
Regards
Carolina Rickenstorff
David Goodmanson
on 18 Feb 2020
Edited: David Goodmanson
on 19 Feb 2020
Hi Carolina,
Yes you can certainly use just H0[+] and H1[+] and avoid using K. From the expressons in the answer, what you are working toward is
mu1*J1(k1*a) / (k1*J0(k1*a)) = mu2*H1(k2*a) / (k2*H0(k2*a)) [1]
where (also including k1)
k1 = sqrt(n1^2*k^2 - kz.^2)
k2 = sqrt(n2^2*k^2 - kz.^2) % imaginary
When solving Maxwell's equations in each region, k2 naturally emerges in parallel with k1. I probably should have included k2 earlier in the answer even though it wasn't used (except the verification at the end). Given the way Matlab defines besselh (that's important), [1] is real on both sides. So
D = mu2*H1(k2*a) / (k2*H0(k2*a))
works. It didn't work in your expression because you were substituting k2 = i*g2 everywhere but you left in one k2 by mistake. So your expression comes out imaginary, incorrectly.
If you are going to use the H functions, you can just get rid of g2 entirely and use the imaginary k2. The reason to use g2 is to get to the traditional K functions. People invented those functions because 150 years ago they didn't want factors of i floating around that eventually cancel to give something real. That's understandable. Instead of something like i*H(i*x), you have K(x). Real in, real out. Another reason is that K satisfies a differential equation related to the wave equation. Another reason is that K falls off exponentially and it's important to identify a solution like that. For the fiber, you want the evanescent wave in region 2 to fall off exponentially with r.
Computationally, K could be faster, I don't know, because H is working with a complex (although pure imaginary) expression.
I'll think about books at lunch. Perhaps anyone else who sees this will have some ideas, aside from AMS55 (more on that one later).
p.s. I am fine with 'David'.
Carolina Rickenstorff
on 19 Feb 2020
Edited: Carolina Rickenstorff
on 21 Feb 2020
Dear Mr. Goodmanson,
Thanks for your expertise. I'm eager to plot the radial fields along the r coordinate and see how they connect smoothly. However, I again have trouble :(. This time with respect to the amplitude factors.
Thanks to you I have the values 
. Since I'm working with the TE mode I have only the fields 
, 
and unknown amplitude factors B,D.
. Since I'm working with the TE mode I have only the fields , and unknown amplitude factors B,D.Given 
and 
m the amplitude factor D is
and m the amplitude factor D is
and I obtained a complex value
B=1;
D=-B.*besselj(0,k1root.*a)./besselh(0,1,k2root.*a); %%%% minus sign for to compensate Matlab notation of besselh, isn't?
I suppose it should be equal to the matching amplitude factor α in your original answer. I'm confused with the notation using A, B, C, D.
Sorry for the inconvenience.
Carolina R
David Goodmanson
on 20 Feb 2020
Hi Carolina,
It looks like maybe you are discovering the advantages of using K. However, H[+] will still work fine. I don't know what you are referring to with A,B,C,D but C and D were something that I just made up for the answer. There is something called the abcd method but I am not using that so I hope there is no confusion there. Here C and D are not amplitudes as in abcd but just ratios to help find the roots.
The very first set of equations in the original answer seemed to imply that that there are no additional constant conversion factors between E and H, but that is not the case. I modified the answer to point that out. The equations still work to find C and D.
E: G*J1(k1*a) = alpha*G*H1(k2*a) [*]
H: G*(1/mu1)*k1*J0(k1*a) = alpha*G*(1/mu2)*k2*H0(k2*a)
% not all the constant conversion factors to get from E to H are shown.
% these cancel out when doing equality of H between the two regions.
Looking at these, everything on the left hand side is obviously real. k2 is imaginary, and it turns out that for imaginary argument i.e. (k2*a) then
H1 is real. H0 is imaginary, so k2*H0 is real. This means that alpha is real.
The C and D I was using are C ~ J/J and D ~ H/H. Here there is no ratio like you have, J /H, which there might be in an abcd method. Anyway, divide the equations above find the same C and D as before, look for C = D as a function of kz (but see below) on the plots, find kzroot, then k1root and k2root, plug those into [*] above, then solve that equation for alpha.
Having plugged alpha, k1root, k2root ito [*] you can get E(r) by replacing a with r, then plot E for r<=a and for r>a. You should be able to see an E field that matches up at the boundary. Its derivative matches up also, but I believe that is just a coincidence for the axially symmetric case. All that is required is that E matches up.
One other thing needs mentioning. In the example I picked the smallest kzroot off the plot. However, it's the largest kzroot that gives the smallest value of k1root. That's the fundamental TE mode. It makes sense to momentarily normalize kz better and plot
kzk = kz/k;
figure(1)
plot(kzk,C,kzk,D)
grid on
In that case the for the largest root on the plot I got
kzkroot = 1.39713991855
and then kzroot = k*kzkroot and proceed from there.
The conversion factors for E to H (including an important factor of i) are not here so that might be the next thing to look at.
Note: the mu's in the equations so far are intended to be relative mu’s, i.e. mu1 = mu1_actual / mu0 etc.
p.s. I will try this again. I would be fine to go by DG instead of Mr. ...
Carolina Rickenstorff
on 21 Feb 2020
Edited: Carolina Rickenstorff
on 21 Feb 2020
David Goodmanson
on 22 Feb 2020
Hi Carolina,
you are definitely getting closer, but I will list a few comments
[1] it's never a bad idea to check out some of the intermediate details. Your list of roots has four of them listed twice, and the one with the largest kz value not listed at all. That's important because the largest kzroot leads to the smallest k1root, and the smaller k1root is, the fewer oscillations there are in J0 and J1. The smallest k1root describes the lowest, fundamental mode. I believe your plots are for the smallest kzroot --> largest k1root.
[2] I think you should retain some kind of comment that stuff like this
Cm = mean(abs(C));
C(abs(C)>3*Cm) = nan;
is for plotting purposes only. In a few months the meaning of this could seem pretty mysterious.
[3] No need to scrimp on points. N = 512 is ok, but when you start zooming in on the plots, N = 10000 is better. No need for a power of 2 either.
[4] The Hz plot works, but the other two have problems. That's because for each of those, the r<a part is imaginary and the r>a part is real. 'Plot' plots only the real part by default, so you get flat line plots for r<a, and warning messages. In this scheme of things where Hz is real, both Hr and Ephi need to be imaginary everywhere, so you need to check the form of the solutions. Once everything is imaginary then of course you will have to plot imag(whatever).
Your k2 is real as it should be to be an argument for K, but k2 started out imaginary and then you implemented k2root=imag(k2root) to make it real. I deliberately used two different variable names, k2 and g2, to avoid this notational ambiguity. Anyway if you add appropriate factors of i to the r>a parts (hopefully not just by trying +-i to see what makes the plots work!) you will get good plots.
Carolina Rickenstorff
on 25 Feb 2020
Edited: Carolina Rickenstorff
on 25 Feb 2020
Walter Roberson
on 25 Feb 2020
Thank you David for the nice breakdown!
David Goodmanson
on 27 Feb 2020
Thank you, Walter, for your comment.
Carolina,
Looks like you got this thing solved. Since you asked,
yes the largest possible order here is five, the number of roots shown on the plot. (Sometimes people start counting at 0 instead of 1, but five altogether). The limits
kzmax = n1*k; kzmin = n2*k;
are there for plotting purposes that's true, but also because outside of that range either k1 goes imaginary or (in the convention you are using), (i*k2) goes real. Either of those produces a bad nonphysical solution.
Hope you don't mind but I have a couple more comments. First of all, the code uses w*mu0 to find Ephi which is correct numerically but I think it better illustrates the physics to use
mu0*w = mu0*c*k = 120pi*k = 377*k
where mu0*c = 377 ohms is the impedence of free space. Then
Ephi1 =-1i.*b.*(120*pi)*(k/k1root).*bessellj(...
i.e., E [V/m] is a resitance times a dimensionless quantity (k/k1root) times H [A/m].
Using ind_roots=find(abs(C-D)<=0.5e-10); will work most of the time but is not the most reliable of methods. Multiple roots make it a bit inconvenient but a zero finder such as fzero is preferable.
In the definition of kz, 500000 points is way more than too much. It's only there because I temporarily upped the number in order to inspect the transition point r = a and then forgot to put it back to something better. Maybe 1000 or whatever you think, makes sense.
More Answers (2)
Carolina Rickenstorff
on 28 Feb 2020
Edited: Carolina Rickenstorff
on 28 Feb 2020
0 votes
David Goodmanson
on 2 Mar 2020
Hi Carolina,
do you have a reference for where you read that dielectric fibers have no cutoff frequency? I would be interested to see who is maintaining that. If you run your code and look at the plot of C vs D, you can see that D is always negative. D is a function of k2*a, and allowed values of k2*a depend on the details of this problem, but you can use any positive value you like for k2*a and D is always negative. That means for a solution,
C = besselj(1,k1*a)./(k1.*besselj(0,k1*a))
has to be negative. But C is only negative if J0 and J1 have opposite sign. If you plot these out, both start out positive and fiirst have opposite sign when J0 goes negative at k1*a > 2.4048, the first root of J0. The actual solution of the problem lies higher than that, but at the very least k1*a > 2.4048. Then from
k1 = sqrt(n1^2*k^2 - kz.^2);
if follows that n1*k*a > k1*a > 2.4048
and since w = c*k there is a cutoff frequency. This exact cutoff frequency turns out to be larger, but this demo does does show the existence of a cutoff frequency.
For the transmitted power, you have to integrate the Poynting vector ExH over the area of the fiber. (Let's say you don't calculate the transmitted power in the cladding since it is not really available for use). This means doing an integral of
(amplitudes)*J1(k1*r)^2 * 2*pi*r dr which is over the area of the fiber. Fortunately that integral is doable analytically and Matlab even has it. Try
syms r k1 a
int(2*pi*r*besselj(1,k1*r)^2,0,a)
2 Comments
Carolina Rickenstorff
on 2 Mar 2020
David Goodmanson
on 3 Mar 2020
Hi Carolina,
thanks very much for sending for the chapter from the book. Thanks also for the compliment, but I think I have reached the end of my knowledge of optical fibers. The solution I gave was related to m = 1 on p. 44, although his angular dependence is different. My impression is that he is playing fast and loose with some of the equations and I don't know that his boundary conditions are even correct. And he does not show any expressions for what the E and B fields actually are. Balanced against this impression is the fact that it's a Springer book in its second edition, so how incorrect could it be?
Everyone agrees that there is a transverse 01 mode with no lower frequency cutoff, which unfortunately my answer does not apply to (although I think the answer has a correct description of certain other modes). I have not yet found an expression for the fields for 01, so it looks like a trip to the university library is on the agenda.
Categories
Find more on Mathematics in Help Center and File Exchange
on 17 Feb 2020
on 3 Mar 2020
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
