See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/220656622

Eigenvalue Problems in Fiber Optic Design

Article  in  SIAM Journal on Matrix Analysis and Applications · January 2006

DOI: 10.1137/S0895479803432708 · Source: DBLP

CITATIONS READS
17 769

1 author:

Linda Kaufman
William Paterson University
66 PUBLICATIONS   8,562 CITATIONS   

SEE PROFILE

Some of the authors of this publication are also working on these related projects:

banded symmetric indefinte matrices View project

least squares View project

All content following this page was uploaded by Linda Kaufman on 03 May 2016.

The user has requested enhancement of the downloaded file.

 SIAM J. MATRIX ANAL. APPL. 
c 2006 Society for Industrial and Applied Mathematics
Vol. 28, No. 1, pp. 105–117

EIGENVALUE PROBLEMS IN FIBER OPTIC DESIGN∗

Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

LINDA KAUFMAN†

Abstract. Maxwell’s equation for modeling the guided waves in a circularly symmetric fiber
leads to a family of partial differential equation–eigenvalue systems. Incorporating the boundary
condition into a discretized system leads to an eigenvalue problem which is nonlinear in only one
element. In fiber design one would like to determine the index profile which is involved in Maxwell’s
equation so that certain optical properties, which sometimes involve derivatives of the eigenvalues,
are satisfied. This contribution discusses how to handle the nonlinear eigenvalue problem and how
to determine derivatives of the eigenvalue problem.

Key words. eigenvalue, nonlinear, Rayleigh quotient

AMS subject classifications. 15A18, 3FP30, 3FQ60

DOI. 10.1137/S0895479803432708

1. Introduction. In this paper we consider several computation issues arising

in fiber optic modeling. Assuming a fiber is perfectly straight, circular, and uniform
along its length, then Maxwell’s equations for guided waves of the fiber can be reduced
to a family in m of problems of the form
   
1 d dx m
(1.1) r + ω η (r, ω) − 2 x = β 2 x.
2 2
r dr dr r
The index of refraction profile η(r, ω) is in some regions a piecewise constant function,
as in Figure 1.1, and can be parameterized by several design parameters relating to
the widths and heights of each region. In (1.1), ω is a specified frequency and m is a
specified mode number. The finite element method converts this family of differential
equations to a family of symmetric tridiagonal eigenvalue problems in ω and m
(1.2) A(ω, m)x = μx,
where one wishes to find the positive eigenvalues and their corresponding eigenvectors.
The eigenvalue μ corresponds to β 2 in (1.1). To simplify our notation we let A
represent one member of this family. Usually the waveform is truncated at some
radius beyond the core of the fiber, and the boundary condition is expressed as the
mth order modified Bessel function of the second kind [9]. This changes the eigenvalue
problem in (1.2) to the form
(1.3) (A + s(μ)en eTn )x = μx,
where s(μ) involves the appropriate Bessel functions, A is an n × n matrix, and en
is the last column of the n × n identity matrix. In section 2, we describe several
algorithms that can be used to solve (1.3).
As explained in [6], for a particular index profile η(r, ω), one is interested in
various integrals of the modes (the eigenvectors), the dispersion
∂2β2
(1.4)
∂λ∂ω
∗ Received by the editors August 1, 2003; accepted for publication (in revised form) by B. T.

Kǎgström September 24, 2005; published electronically March 17, 2006.

http://www.siam.org/journals/simax/28-1/43270.html
† Computer Science Department, William Paterson University, Wayne, NJ 07470 (KaufmanL@

wpunj.edu).
105
 106 LINDA KAUFMAN

3
x 10
6
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

4
refractive index profile η

2
0 5 10 15 20 25
radius (microns)

Fig. 1.1. A typical index of refraction profile η(r, ω).

and the dispersion slope

∂3β2
(1.5) ,
∂λ∂ω 2
where λ is the excitation energy and 2πc = λω with the frequency ω is measured in
radians per second, and c is the speed of light. More to the point, one is interested
in finding the index profile η which might produce certain values of the dispersion
and the dispersion slope at particular values of m and ω. If the index profile η is
partially defined by a parameter p, then to use an optimization package to determine
p typically requires derivatives such as
∂4β2
(1.6) .
∂λ∂ω 2 ∂p
Numerical differentiation of these quantities has proved unsatisfactory as small changes
in a variable are not necessarily accurately reflected in the computed dispersion. Thus,
analytic differentiation is required. In section 4, we review formulas for derivatives of
eigenvalues and indicate how the amount of computation can be reduced by noticing
common subexpressions.
Note that by defining η(r, ω) we are actually specifying the chemical composition
of the radial layers for a fiber. A fiber suitable for, say, underwater transmission would
not be appropriate for a local area network or to splice into an existing network to
restore a degraded signal. Also various manufacturers tend to emphasize different op-
tical properties for a given application. When using an optimizer to determine η(r, ω),
for each function evaluation one may have to solve up to 30 nonlinear eigenvalue prob-
lems for various values of m and ω in (1.1). Because there are sometimes multiple
local minima to a particular optimization problem, the optimizer is often called sev-
eral times with different starting guesses to find a design that is manufacturable and
 EIGENVALUES IN FIBER OPTIC DESIGN 107
0
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

–0.5

S(eigenvalue)
–1

–1.5

–2

–2.5
0 0.5 1 1.5 2 2.5 3 3.5 4
eigenvalue

Fig. 2.1. The function s(μ) for varying eigenvalue μ.

more suitable. A production run may require solving thousands of nonlinear eigen-
value problems and then finding derivatives of the computed eigenvalues with respect
to the parameters that define the matrix.

2. The nonlinear problem. In this section three methods are considered for
finding the positive eigenvalues and their corresponding eigenvectors of the nonlinear
eigenvalue problem (1.3). For m = 0, s(μ) is displayed in Figure 2.1. Usually only
one or two eigenvalues are positive.
The first method is a simple Picard iteration.
Picard iteration
Let B = A + s(0)en eTn
Until convergence
Find μ, a specific positive eigenvalue of B.
Reset B to A + s(μ)en eTn
The matrix B is also tridiagonal and differs from A only in its bottom right
element. The eigenvalues of B were determined using the bisection code in EISPACK
[12] modified as in Kaufman [7] so that the inertia of the complete matrix was not
computed when the inertia of only part of the matrix sufficed. However, because there
were such good guesses for the eigenvalues and eigenvectors from previous iterations
of the optimization procedure it was much more efficient to use a Rayleigh quotient
iteration to determine μ.
Linear Rayleigh quotient algorithm for a fixed B
(in the inner loop of the Picard iteration)
Determine a lower bound μl and an upper bound μu for the desired eigenvalue μ.
Set μ0 to μl .
Set x = e , the vector of all 1’s or a good guess if available.
Set k to 0.
Until convergence
Solve (B − μk I)y = x and determine the inertia of (B − μk I)
If the inertia claims that μk is less than μ, reset μl to μk
If the inertia claims that μk is greater than μ, reset μu to μk
Set α = y T By/(y T y)
 108 LINDA KAUFMAN

Table 2.1
The growth of | s (μ) | as μ approaches zero.
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

the eigenvalue μ | s (μ) |

2.0000 0.3631
1.0000 0.5219
0.5000 0.7559
0.2500 1.1054
0.1250 1.6358
0.0625 2.4543

If α < μl , set μk+1 = .9μl + .1μu

If α > μu , set μk+1 = .9μu + .1μl
If μl ≤ α ≤ μu set μk+1 = α
set x = y/(y T y)
increment k.
In the algorithm outlined above, the inertia of the system and the solution y can
be determined using the Bunch algorithm for tridiagonal systems [1]. The choice of
.9 and .1 were arbitrary and just a heuristic.
Although the inner Rayleigh quotient iteration has cubic convergence, the outer
iteration is a fixed point scheme and is generally linear convergent. If the outer
iteration is as μ = f (μ), then f (μ) is an eigenvalue of the B matrix. The rate of
df
the convergence, and whether it converges, depends on dμ . Let x be the eigenvector
df ds
corresponding to f . Then f x = (A + s(μ)en eTn )x, which implies dμ x = T
dμ en en x,
giving
df ds 2
(2.1) = x /(xT x).
dμ dμ n

Because x2n /(xT x) < 1, convergence is guaranteed as long as | dμ

ds
| < 1. For our fiber
ds
optic problem, Table 2.1 gives dμ for various values of μ.
Thus as μ approaches zero, convergence is dubious.
In practice a major problem with the Picard iteration approach is speed. One
potentially spends time in the inner loop of the Rayleigh quotient iteration for the
wrong problem. This pitfall is especially noticeable for problems where the eigenvalue
is near zero and the function s(μ) is rather steep, as in the example at the end of this
section.
One way around the sluggishness of the Picard iteration is to mimic the derivation
of the Rayleigh quotient iteration for the nonlinear problem. In the Rayleigh quotient
iteration for a matrix A, the next iterate μk+1 is chosen so that

(2.2) (A − μI)y2

is minimized, which occurs when 2y T Ayμ + μ2 is minimized or when μ = y T Ay/y T y.

For the nonlinear Rayleigh quotient case one would like to minimize

(2.3) (A + s(μ)en en T − μI)y2 .

If z = (A + s(μ)en en T )y, then the required μ minimizes z T z − 2μz T y + μ2 y T y,

which occurs when
zT y − zT z
(2.4) μ= .
yT y − z T y
 EIGENVALUES IN FIBER OPTIC DESIGN 109

Since z  = s (μ)en en T y = s (μ)yn en , from (2.4)

z T y − s (μ)yn zn
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

(2.5) μ= .
y T y − s (μ)yn 2
With the linear Rayleigh quotient algorithm given above, the bounds on the
eigenvalue estimates ensure that the algorithm is converging to the correct eigenvalue.
The bounds can be obtained by examining the Bunch–Parlett decomposition [2], which
dictates that one can find a lower triangular matrix L and a block diagonal matrix D
with 1 × 1 blocks and 2 × 2 blocks such that

(2.6) B − μI = LDLT .

If D has i 2 × 2 blocks and j positive 1 × 1 blocks, then B has i + j eigenvalues larger

than μ. In a nonlinear Rayleigh quotient iteration one would like the same bounding
mechanism to be preserved even if the element Bn,n is changed in each iteration. The
following theorem treats the largest eigenvalue and a similar proof can be used to
handle other eigenvalues.
Theorem 2.1. Let s(μ) be a negative decreasing function. If μu is an upper
bound of the largest eigenvalue of the matrix B0 = A + s(μ0 )en en T , where μ0 is an
initial guess of the eigenvalue, then it is an upper bound on the largest eigenvalue
μmax of the problem (A + s(μ)en en T )x = μx.
Proof. If μ0 ≥ μmax , the theorem is trivially true. So assume μmax > μ0 . Let
D be the diagonal matrix in the Bunch–Parlett decomposition [2] of B0 − μu I. If μu
is an upper bound of the largest eigenvalue of B0 , then D is diagonal with negative
elements. Let B̃ = A + s(μmax )en en T , where μmax is the solution to the nonlinear
problem, and let D̃ be the diagonal matrix in the Bunch–Parlett decomposition of
B̃ − μu I.
The first n−1 diagonal elements of D and D̃ are identical and d˜n = dn +s(μmax )−
s(μ0 ). If μmax > μ0 , then s(μmax ) < s(μ0 ), so d˜n is negative and μu is an upper bound
for the nonlinear problem.
Theorem 2.2. Let s(μ) be a negative decreasing function. If μl is a lower bound
for the largest eigenvalue of the matrix B0 = (A + s(μ0 )en en T ), then it is a lower
bound for the largest eigenvalue μmax of the problem (A + s(μ)en en T )x = μx, as long
as μmax ≥ 0.
Proof. Let D be the block diagonal matrix in the Bunch–Parlett decomposition
of B0 − μl I. If there is a 2 × 2 block above the bottom two rows or if there is a
positive element on the diagonal above the bottom right corner element, then no
matter how s is changed, D will signal that μl always is a lower bound for the largest
eigenvalue. Thus, one may assume that dn is positive. Moreover, if μmax ≥ μ0 , the
theorem is also trivially true so assume μmax < μ0 . If D̃ is the diagonal matrix in
the Bunch–Parlett decomposition of A + s(μmax )en en T , where μmax is the solution
to the nonlinear problem, then d˜n = dn + s(μmax ) − s(μ0 ). If 0 ≤ μmax < μ0 , then
s(μmax ) > s(μ0 ), so d˜n will be positive and μu will be a lower bound of the nonlinear
problem.
From (2.5) one gets the following algorithm.
Nonlinear Rayleigh quotient
Determine a lower bound μl and an upper bound μu for the desired eigenvalue μ.
Set μ0 to μl .
Set x = e , the vector of all 1’s if a good guess is not available.
Set k to 0.
Until convergence
 110 LINDA KAUFMAN

Set B = A + s(μk )en eTn .

Solve (B − μk I)y = x and determine the inertia of (B − μk I)
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

If the inertia claims that μk is less than μ, reset μl to μk

If the inertia claims that μk is greater than μ, reset μu to μk
Set z = By.
T 
Set α = zyTy−s (μ)yn zn
y−s (μ)yn 2
If α < μl , set μk+1 = .9μl + .1μu
If α > μu , set μk+1 = .9μu + .1μl
If μl <= α <= μu set μk+1 = α
set x = y/(y T y)
increment k.
A third method, suggested by Cowsar [4], is designed to find a root of the function

(2.7) f (γ) = γ − s(μ(γ)),

where μ satisfies

(A + γen eTn − μI)x = 0.

Newton’s method for determining the root of (2.7) determines the new approximation
of γnew using the formula

γ − s(μ(γ))
(2.8) γnew = γ − df
.
dγ

Note that
df ds dμ
=1−
dγ dμ dγ

and if μx = (A + γen eTn )x, then

dμ
x = en eTn x.
dγ

Because x can be chosen such that xT x = 1,

dμ
= (xT en )2 .
dγ

Thus each iterate of Newton’s method for minimizing (2.7) would set the new approx-
imation of γnew to

γ − s(μ(γ))
(2.9) γnew = γ − ,
1 − dμ
ds
(xT en )2

where x is the normalized eigenvector corresponding to the eigenvalue μ of A + γen eTn .

Because Newton’s method is not globally convergent, it is wise to determine
initially an interval containing the root and to use a bisection technique if the new
iterate does not fall within the interval. If one is seeking a positive eigenvalue μ, then
0 is an upper bound for γ because f (0) > 0 since s(μ) is negative if μ is positive. A
lower bound for γ can be found if the Bunch–Parlett decomposition yields a diagonal
 EIGENVALUES IN FIBER OPTIC DESIGN 111

matrix D whose first n − 1 diagonal elements are negative but whose right corner
element dn is positive. In this case the matrix has one nonnegative eigenvalue and
if γ = −dn , then the largest eigenvalue of A + γen eTn is 0. For m = 0, s(0) = 0,
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

which means f (−dn ) < 0. If one of the first n − 1 elements of the D matrix in the
Bunch–Parlett decomposition is negative, then any very large negative number can
be used as a lower bound for γ.
In the formal definition of Newton’s method given below, to find one positive
eigenvalue of the nonlinear eigenvalue problem an upper and a lower bound are first
determined and subsequently the bounds are adjusted so that the root of f is always
bracketed. If the Newton step is outside the bound, then a step that is 90% of the
distance to the bound is taken. If the chosen step does not decrease the norm of f ,
the step size is reduced by 4. To determine several eigenvalues, the algorithm is used
first to determine the largest eigenvalue, the bounds are adjusted to find subsequent
eigenvalues, and the linear Rayleigh quotient solver is ask to find a specific numbered
eigenvalue.
Safeguarded Newton’s method
Set γu , an upper bound for γ to 0.
Find γl , a lower bound for γ as follows:
Determine D, the diagonal matrix of the Bunch–Parlett decomposition
of A.
If none of the elements of D are positive, there is no positive eigenvalue.
If one of the first n − 1 elements of D is greater than 0, set γl to a
very large negative number else set γl = −dn
If an initial eigenvalue μ0 is available, set γ0 = s(μ0 ), otherwise set γ0 = γu ;
Set k to 0.
Until convergence
Find μ > 0 and x such that (A + γk en eTn )x = μx using the linear
Rayleigh quotient algorithm
Set f = γk − s(μ); r = abs(f )
If f > 0
if γk ≤ γu set γu = γk
If f < 0
if γk ≥ γl set γl = γk
if (k > 1 and r > minr)
Set γchange = γchange /4
Set γk+1 = γk + γchange
else
Set minr to r
Set γk+1 = γk − 1−γdsk −s(μ)
(xT e )2
dμ n

If γk+1 < γl , set γk+1 = .9γl + .1γu

If γk+1 > γu , set γk+1 = .9γu + .1γl
Set γchange = γk+1 − γk
increment k

3. Computational experiments. To indicate the effectiveness of the three

algorithms outlined above, we show their effectiveness on several problems. The
problem in section 3.1 is a small example with a tunable parameter which indicates
that the Picard iteration and the nonlinear Rayleigh quotient may converge linearly
and slowly. The second example in section 3.2 is derived from a fiber optics problem
 112 LINDA KAUFMAN

Table 3.1
The course of the Picard iteration, the nonlinear Rayleigh quotient, and the Newton approach
with z = .5 in (3.1).
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Picard iteration Nonlinear Rayleigh Newton

iterate inner iterate inner derived inner
iterations iterations eigenvalue iterations
.3589622 4 .2865810 1 .3589622 4
.3589538 2 .3577595 1 .3589538 2
.3589538 1 .3589538 1 .3589538 1
.3589538 1

in the literature.
3.1. A small example. In this section we consider the following small problem:
⎧
⎨ −2 + z, 1 ≤ i ≤ 5,
(3.1) ai,i = −2 − z, i = 6, 7,
⎩
−2, 8 ≤ i ≤ 12,

and for 1 ≤ i ≤ 11

(3.2) ai+1,i = ai,i+1 = (i + .5)/ (i + 1) × (i + 2).

The function s(μ) involved computing the modified Bessel function of the second
kind calculated using ACM Algorithm 484 [3].
For z = .5, the largest eigenvalue is about .359, where s(μ) is not very steep,
and all three methods behaved rather well, as Table 3.1 indicates. The rate of
ds 2
convergence for the Picard iteration is based on dμ xn /(xT x). For this value of z,
as the algorithm converged, dμ was about −1.08, xn was about −.00242, so that
ds

−6
dμ xn /(x x) ≈ −6.36 × 10
ds 2 T
, suggesting fast convergence. Note for Newton’s method
the eigenvalue corresponding to γ is given.
For z = .081, the largest eigenvalue is about 3.2 × 10−4 , where the s function is
ds
rather steep. For this value of z, as the Picard algorithm converged, dμ was about
−3
−3.23, xn was about −.047, so that dμ xn /(x x) ≈ −5.56 × 10 , which is approx-
ds 2 T

imately the square root of the previous example. The Picard iteration algorithm
required 9 outer iterations and 42 inner iterations, which is far more than 14 required
by the nonlinear Rayleigh quotient and the 13 required by the Newton approach. The
course of the algorithms is shown in Table 3.2. Note that the first three iterations
of the nonlinear Rayleigh quotient iteration produced iterates that were negative and
then were forced to be positive.
Results like those obtained in Table 3.2 for the Picard iteration have stimulated
the search for other algorithms that might not have linear convergence. It was hoped
that the nonlinear Rayleigh quotient algorithm would be cubically convergent and
one could mimic the proof of cubic convergence given in Parlett [11], but when an
ds
eigenvalue approached zero, one could not easily bound dμ . Moreover, problems like
that given in (3.2) suggested that the algorithm was linearly convergent. Table 3.3
shows the rate of convergence of the Picard and nonlinear Rayleigh quotient iteration.
In the table μ̂ refers to the solution.
3.2. A fiber optics problem. The nonlinear Rayleigh quotient algorithm and
the Newton approach were also applied to a more realistic problem that was posed
by Lenahan and Friedrichsen [10] for m = 1 involving a core region of silicon dioxide
 EIGENVALUES IN FIBER OPTIC DESIGN 113
Table 3.2
The course of the Picard iteration, the nonlinear Rayleigh quotient, and the Newton approach
with z = .081 in (3.1).
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Picard iteration Nonlinear Rayleigh Newton

iterate inner iterate inner derived inner
iterations iterations eigenvalue iterations
.6770570×10−3 5 1.0 1 .6770570×10−3 5
.2831247×10−3 11 .1 1 .3187768×10−3 3
.3294611×10−3 3 .01 1 .3229664×10−3 3
.3219722×10−3 3 .8701964×10−4 1 .3229683×10−3 2
.3231226×10−3 2 .3682030×10−3 1
.3229445×10−3 2 .3161693×10−3 1
.3229720×10−3 2 .3240255×10−3 1
.3229677×10−3 7 .3228049×10−3 1
.3229684×10−3 2 .3229936×10−3 1
.3229683×10−3 4 .3229644×10−3 1
.3229683×10−3 1 .3229689×10−3 1
.3229682×10−3 1
.3229683×10−3 1
.3229683×10−3 1

Table 3.3
Rates of convergence with z = .081.

Picard iteration Nonlinear Rayleigh

(μk − μ̂)/(μk−1 − μ̂) (μk − μ̂)/(μk−1 − μ̂)
−.113 .100
−.163 .097
−.153 −.024
−.155 −.192
−.154 −.150
−.155 −.155
−.162 −.155
−.167 −.155
−.154
−.154
−.167

doped with germanium surrounded by a region of pure silicon dioxide. The index
profile for the core region had the form

(3.3) η(r, ω) = e(ω) + C(r)h(ω, sign(C(r)),

where e(ω)is the index of pure silicon dioxide, C(r) denotes the dopant concentration,
and h is a function of the Sellmeier coefficients [5] at ω and the reference frequency.
We started with an example where the radius of the fiber Rf and the radius of the
core Rc satisfied the relation Rf = 6Rc and C(r) was nonnegative and defined by

((1 − 2δ(r/Rc )α )/(1 − 2δ))1/2 − 1), 0 ≤ r ≤ Rc ,
(3.4) C(r) =
0, r > Rc .

We solved (3.4) for several values of λ = 2π × 4.0/(Rc ω), δ, and α. In an opti-

mization problem one may wish to vary α and δ and solve the eigenvalue problem for
about 20 values of λ. A uniform grid was used in the finite element discretization and
Rc was set at 100 units and Rf at 600 units. The matrix eigenvalue problem derived
 114 LINDA KAUFMAN

Table 3.4
Inner iteration count for the nonlinear Rayleigh quotient method and the Newton approach on
(3.4).
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

λ α δ eigenvalue nonlinear RQI Newton

0.85 25 .003 2.636370×10−4 6 7
1.1 25 .003 4.216316×10−6 10 14
0.85 20 .003 2.525251×10−4 6 7
1.1 20 .003 1.110249×10−6 14 15
0.85 25 .004 5.486164×10−4 5 6
1.1 25 .004 1.046444×10−4 5 8

Table 3.5
Inner iteration count for the nonlinear Rayleigh quotient method and the Newton approach on
the augmented model with λ = .85, α = 25, and δ = .003.

ρ σ eigenvalues nonlinear RQI Newton

−0.01 0.005 1.312023×10−3 6 7
2.824626×10−5 10 19
−0.008 0.005 1.319503×10−3 6 7
4.442786×10−5 16 20
4.807390×10−6 4 7
−0.01 0.008 2.366140×10−3 7 8
7.768107×10−4 17 18
−0.01 0.001 5.740407×10−5 11 15

from this problem involved a tridiagional matrix whose diagonal elements, ai,i , were
given by
⎧
⎨ −2 + ω 2 (η 2 − e(ω)2 ) − (1/i)2 , 1 ≤ i ≤ 100,
ai,i = −2 − (1/i)2 ,  100 < i < 600,
⎩
−2 − (1/i)2 + (1 + .5/600) × 601.5/ (600) × (601), i = 600,

and for 1 ≤ i < 600


(3.5) ai+1,i = ai,i+1 = (i + 1.5)/ (i) × (i + 1).

For the values of λ, α, and δ given in Table 3.4, there was only one positive
eigenvalue. With the smaller eigenvalues both algorithms had trouble. For λ = .85,
we used the approximation e = 1.45291 and h = 1.44943 in (3.3). For λ = 1.1, we
used e = 1.4969 and h = 1.4201.
We then augmented the model in (3.4) with a layer of fluorine doped silicon
dioxide of constant concentration ρ and then a layer of germanium doped silicon
dioxide of constant concentration σ. Each layer had the same width as the core
region. We kept the same width of the fiber, used the same uniform finite element
discretization, and set the wavelength λ at .85, α = 25, and δ = .003. As we varied
the concentrations, the number of positive eigenvalues changed. For the fluorine layer
the dopant concentration was negative and h in (3.3) was 1.456475. In Table 3.5
the number of Rayleigh quotient iterations is given for several values of σ and ρ. In
general the nonlinear Rayleigh quotient seemed to be slightly faster than the Newton
technique but the differences were rather inconsequential.
Usually when an eigenvalue was not close to zero, the nonlinear Rayleigh quotient
algorithm required the same effort as the initial iteration of the Newton approach.
In the production code a polyalgorithm was formed that initially used the Rayleigh
 EIGENVALUES IN FIBER OPTIC DESIGN 115

quotient algorithm, and if that did not work after a specified number of iterations, its
best approximation was used to determine an initial γ in the Newton procedure.
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Because we were solving a sequence of eigenvalue problems in which there were

often only small changes in the matrix, good starting approximations for eigenvalues
and eigenvectors were usually available from previous problems.
4. Analytic derivatives. The function η in (1.1) is defined by a number of pa-
rameters pk , k = 1, 2 . . . K, giving the shape of the initial region and the concentration
and width of the various layers of the fiber. Often the parameters pk are requested
∂2 β2 ∂3 β2
such that the dispersion ∂λ∂ω and the dispersion slope ∂λ∂ω 2 have prescribed values

where ω = 2πc/λ for specific value of ω and m in (1.1). In a typical optimization

problem to find an optimal shape of η, one might have to supply the dispersion for
about 30 combinations of ω and m for each function and derivative evaluation. A
signal traveling down a fiber tends to spread out over various wavelengths. Tradi-
tionally, about every 40 kilometers the signal was restored electrically. The concept
of a dispersion compensating fiber for existing fiber is to splice in a small length of
fiber with negative dispersion that is specially created to restore the signal so that
electrical restoration is not needed. A fiber optics company might also wish to lay
new fiber which has zero dispersion.
If μ = β 2 , one needs first to solve (A + s(μ)en eTn )x = μx and then to determine
∂2μ ∂μ 2πc ∂ 2 μ ∂3μ
the dispersion, 2μ11/2 ∂λ∂ω . Moreover, since ∂μ ∂λ = − ∂ω λ2 , ∂ω 2 and ∂ω 3 must be
computed to evaluate the dispersion and dispersion slope, respectively.
Let us use the shorthand A , μ , and x to denote the derivatives of the elements
of A, μ, and x with respect to ω, respectively. To simplify our notation let B =
(A + s(μ)en eTn ), φ = xT en , and σ = 1 − sμ φ2 . The following lemma provides a
formula for μ , and x .
Lemma 4.1. One can express the derivative of μ with respect to ω as

(4.1) μ = xT A x/σ,

and x satisfies xT x = 0 and

(4.2) (μI − B)x = −(μI − B) x.

Proof. Differentiating the equation Bx = μx yields

(4.3) ((μI − A) − sμ μ en eTn )x + (μI − B)x = 0,

which, when multiplied by xT , implies xT (μI − A) x − sμ μ φ2 = 0. If x is chosen such

that xT x = 1, then

(4.4) μ = xT A x/σ.

Because sμ is always negative, σ = 1 − sμ φ2 is never zero. From (4.3) and (4.1) one
gets (4.2). Unfortunately (μI − B) is singular, but if μ is not a multiple eigenvalue of
B, the vector x can be determined uniquely by using the condition xT x = 0.
Following the approach given in the proof of Lemma 4.1 it is shown in [8] that

the dispersion is given by − 2πc
λ2 μ , where

(4.5) μ = (xT A x + 2xT B  x + α)/σ,

where α = sμμ μ2 φ2 .

 116 LINDA KAUFMAN

Calculating the dispersion slope requires μ . If ρ = φφ((s (μ)μ ) + s (μ)μ μ )
2

and φ = eTn x , then one can express μ as

Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

μ = (xT A x + 6xT B  x + 6x B  x − 6μ x x + ρ)/σ.

T T
(4.6)

Since our aim is to determine the parameters in η that give a prescribed dispersion,
the derivatives of the dispersion with respect to these parameters are needed by most
optimizers. Let us use the convention that μk is the derivative of μ with respect to
the kth design parameter, Ak is the derivative of the A matrix with respect to the kth
design parameter, xk is the derivative of x with respect to its kth design parameter,
etc. Mimicking the proof of Lemma 4.1 suggests that if xT x = 1, then

(4.7) μk = xT Ak x/σ.

It is shown in [8] that if θ contains all the terms in s(μ)k φ2 except the one with
μk , then μk can be expressed several ways, including

(4.8) μk = (xT Ak x + 2(xT B  xk + xT Bk x + xT B  xk + xTk B  x ) + θ)/σ,

where xk satisfies (μI − B)xk = −(μI − B)k x − (μI − B) xk − (μI − B)k x and the
condition xT xk = −xT xk , and by the equation

μk = (xT Ak x + 4xT Bk x + 2x Bk x + 2xT Bk x + θ)/σ,

T
(4.9)

where (μI − B)x = −(μI − B) x − 2(μI − B) x and xT x = −x x .
T

The first formula for μk in (4.8) comes from differentiating (4.5) with respect to
the design parameter, but it can be the less efficient approach. The second formula
in (4.9) reverses the order of differentiation and comes from twice differentiating with
respect to ω the formula for μk in (4.7). The formula in (4.8) requires x, x , and for
the kth design parameter xk and xk . For 20 design parameters one must solve for 42
vectors. The second formula for μk requires x, x , x . Thus for 20 design parameters
the formula in (4.9) requires 3 vectors. Theoretically, (4.8) requires at least twice as
much work as (4.9).
Acknowledgments. I thank Bill Reed, formerly of Bell Labs, for bringing the
fiber optics problem to my attention and I thank Lawrence Cowsar of Bell Labs for
providing the program for determining a good discretization of the pde and the optical
properties of the fiber.

REFERENCES

[1] J. R. Bunch, Partial pivoting strategies for symmetric matrices, SIAM J. Numer. Anal., 11
(1974), pp. 521–528.
[2] J. R. Bunch and B. N. Parlett, Direct methods for solving symmetric indefinite systems of
linear systems, SIAM J. Numer. Anal., 8 (1971), pp. 539–655.
[3] K. H. Burrell, Algorithm 484: Evaluation of the modified Bessel functions K0(Z) and
K1(Z)for complex arguments, Commun. ACM, 17 (1974), pp. 524–526.
[4] L. Cowsar, personal communication, 1999.
[5] J. W. Fleming, Material dispersion in lightguide glasses, Elect. Lett., 14 (1978), pp. 326–328.
[6] S. V. Kartalopoulous, Introduction to DWDM Technology, IEEE Press, Piscataway, NJ,
2000.
[7] L. Kaufman, An observation on bisection software for the symmetric tridiagonal eigenvalue
problem, ACM Trans. Math. Software, 26 (2000), pp. 520–526.
[8] L. Kaufman, Calculating Dispersion in Fiber Optics Design, in preparation.
 EIGENVALUES IN FIBER OPTIC DESIGN 117

[9] T. A. Lenahan, Calculation of modes in an optical fiber using the finite element method and
EISPACK, Bell System Technical J., 62 (1983), pp. 2663–2694.
[10] T. A. Lenahan and H. W. Friedrichsen, Analysis of Propagation over Single-Mode Optical
Downloaded 05/03/16 to 38.98.219.157. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php

Fibers Using the Finite Element Method and EISPACK, Technical Memorandum 82-54541-
34, Bell Labs, Atlanta-Norcross, GA, Nov. 1982.
[11] B. N. Parlett, The Symmetric Eigenvalue Problem, Prentice–Hall, Englewood Cliffs, NJ,
1980.
[12] B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, Y. Ikebe, V. C. Klema, and C.
B. Moler, Matrix Eigenvalues Routines—Eispack Guide. Springer-Verlag, Berlin, 1976

View publication stats