APPROXIMATIONS FOR THE FIBER REFRACTIVE INDEX PROFILE

N.M. DRAGOMIR, S.S. DRAGOMIR, Y.J. CHO, AND G.W. BAXTER

Abstract. Two methods in approximating the fiber refractive index profile that have been recently obtained are reviewed. Two new methods based on the approximation of Stieltjes integral via mid-point and trapezoidal rule are also examined.

1. Introduction

Transversal interferometry of optical fibers has been established in the last twenty years as one of the most useful and accurate tools for refractive index profiling. As shown in [3], since introduction of accurate data reduction formulae by Sochacki [1], the measurement speed and precision depends mainly on the data acquisition techniques.

This, indeed, may be argued if an accurate numerical approximation of the refractive index representation formula obtained in [1]

n(u) := exp [\frac{1}{π} ∫_{u}^{1} \frac{ψ(s)}{\sqrt{s^{2} − u^{2}}} ds] , 0 ≤ u ≤ 1;

is assumed.

For some preliminary studies on the numerical approximation of the refractive index, see [4] and the references therein.

In [9], based on the Taylor’s formula with integral remainder, the authors pointed out a numerical approach in approximating the refractive index in the assumption that the derivatives up to a certain order of the deflection function ψ are known. Since in experiments only discrete values of ψ may be obtained, we consider as the actual deflection function ψ to be the best fitting curve to data (polynomial etc...), which are differentiable up to a certain order and may be obtained, for example, by classical interpolative techniques [2]. A theoretical error analysis has been performed and the uniform convergence of the numerical procedure has been established. Some numerical examples showing a very good accuracy have been provided as well.

Another approximation of the refractive index profile has been obtained in [10] using the trapezoidal rule. This rule has been successfully implemented numerically when certain values of the deflection function are known [10].

Before we present our new approximation for the refractive index, based on the mid-point and trapezoidal rule for Stieltjes integral, some preliminary facts on the phase-stepping interferometry and optical fiber profiling are necessary. To briefly describe this, we follow [3].

1.1. Basics of the Phase-Stepping Interferometry. The intensity distribution, as observed in the exit plane of the interferometer of any kind, can be generally described as

(1.1)  I = I_{0} + I_{c} cos [Δ(x, y) − Φ_{B}]

or, equivalently:

(1.2)  I = I_{0} + I_{c} cos Δ(x, y) cos Φ_{B} + I_{c} sin Δ(x, y) sin Φ_{B}

where I_{0} is the background intensity, I_{c} is the interference pattern amplitude governing the image contrast, Δ(x, y) is the phase change introduced by investigated object, and Φ_{B} is the bias (phase difference between the interfering wave fronts, a quantity characteristic for the instrument). When Φ_{B} is constant, the observed field is an homogeneous one (in absence of Δ(x, y)). When Φ_{B} is a linear function of space coordinates, a fringe pattern appears.

From (1.1) follows that the intensity I in any interferogram is a periodic function of the bias Φ_{B} with 2π period, so it can be alternately represented by a Fourier series in terms of Φ_{B} :

(1.3)  I = a_{0} + ∑_{i=1}^{∞} [a_{i} cos (i · Φ_{B}) + b_{i} sin (i · Φ_{B})]

where the Fourier coefficients are, as usual, given by

(1.4)  a_{i} = \frac{1}{2π} ∫_{0}^{2π} I (Φ_{B}) cos (i · Φ_{B}) dΦ_{B}, i = 1, 2, ...
b_{i} = \frac{1}{2π} ∫_{0}^{2π} I (Φ_{B}) sin (i · Φ_{B}) dΦ_{B}, i = 1, 2, ....

If, for simplicity, we take into consideration only the constant and the first order terms of expansion (1.3) and relating them to the formula (1.2) the exit intensity I can be finally expressed as (see for example [3]):

(1.5)  I = a_{0} + a_{1} cos (Φ_{B}) + b_{1} sin (Φ_{B})

with

(1.6)  a_{0} = \frac{1}{2π} ∫_{0}^{2π} I (Φ_{B}) dΦ_{B} = I_{0},
a_{1} = \frac{1}{2π} ∫_{0}^{2π} I (Φ_{B}) cos (Φ_{B}) dΦ_{B} = I_{1} cos Δ(x, y),
b_{1} = \frac{1}{2π} ∫_{0}^{2π} I (Φ_{B}) sin (Φ_{B}) dΦ_{B} = I_{1} sin Δ(x, y).

From this comparison follows that the ratio of the first order coefficients of the Fourier expansion gives the tangent of the phase change Δ(x, y), hence it can be obtained as (see for example [3]):

(1.7)  Δ(x, y) = arctan (\frac{b_{1}}{a_{1}} ) .

To simplify the computation of coefficients a_{1} and b_{1} the integrals in (1.6) may be replaced by finite sums according to any quadrature rule. For periodic functions the Bessel formulae are recommended. In such a case a_{1} and b_{1} may be approximated by:

(1.8)  a_{1} = \frac{1}{N} ∑_{k=0}^{2N−1} I_{k} cos (\frac{kπ}{N} ) ,

b_{1} = \frac{1}{N} ∑_{k=0}^{2N−1} I_{k} sin (\frac{kπ}{N} ) ,

where I_{k} are the exit intensity distributions observed with bias Φ_{B} = Φ_{k} = \frac{kπ}{N} , and 2N is the number of equidistant points Φ_{k} here the intensity samples are taken. They have to run over the entire period of I, so that I(Φ_{0}) = I(Φ_{2N}).

With the smallest possible number of samples (2N = 4) , the phase of interest can be retrieved from:

(1.9)  Δ(x, y) = arctan [\frac{I_{1}(x, y) − I_{3}(x, y)}{I_{0}(x, y) − I_{2}(x, y)}]

where I_{0}, I_{1}, I_{2} and I_{3} are intensity distributions as observed with Φ_{B} = 0, \frac{π}{2} , π, \frac{3π}{2} respectively. This means that by introducing a \frac{π}{2} step in the bias Φ_{B} the phase function Δ(x, y) can be unequivocally reconstructed from four interferograms.

The technique briefly described above following [3], was first applied by Brunning et al. [5] to study optical surfaces in a Twyman-Green interferometer and is commonly called the phase stepping technique. It has quickly found applications in a number of wavefront measuring devices for optical elements control, surface profiling, recognition of shape and deformation, as well as measurement of vibration amplitudes and phases (see for example [6]). Many other modifications of the original algorithm were also published, see for example the review paper by Creath [7] and the references therein.

1.2. Applications to Optical Fibers Profiling. As shown in [3], the phase stepping technique can also be applied to retrieve the phase information necessary in the refractive index profiling of optical fibers.

We will briefly describe how this can be done.

Suppose that the basic interference formula (1.1) describes the interference pattern as observed in any transverse wavefront shearing interferometer, where the fiber is investigated. Choosing the coordinate system so that the axis of the fiber image is in the x-axis, the x dependence can be eliminated from formulae (1.1) − (1.9) due to the cylindrical symmetry of the object and its homogeneous along the axis. When the wavefront shear in direction perpendicular to the fiber axis exceeds its diameter, two fiber images can be observed and the function Φ(y) = \frac{λ}{2π} Δ(y) represents the optical path length through the fiber. If the shear in the same direction is very small (comparing with fiber’s diameter), then \frac{λ}{2π} Δ(y) is proportional to the derivative of the optical path length Φ(y) in the shear direction, where the image shear s is the proportionality constant (see for example [8]):

(1.10)  \frac{dΦ(y)}{dy} = \frac{1}{s} · \frac{λ}{2π} · Δ(y).

In this case the fiber is said to be observed in differential interference contrast (DIC).

In any case the derivative (with respect to the radial position in observation plane) of the optical path length Φ(y) through the fiber can be either directly obtained or numerically constructed from phase data Δ(y) retrieved from four intensity distribution measurements as required by formula (1.9) . This derivative is essential in the refractive index profile retrieval algorithm proposed in [1]. In the first step of this algorithm it serves to calculate the deflection function ψ(ỹ)

(1.11)  ψ(ỹ) = − arcsin [\frac{dΦ(y)}{dy}]

where ỹ is the incident ray position in object plane (with illumination perpendicular to the fiber axis) and can be related to the image plane coordinate y by an appropriate mapping relation depending on the experiment conditions [1], [5].

The deflection function calculated in such of way is parametrically related to the fiber refractive index profile

(1.12)  n(u) := exp [\frac{1}{π} ∫_{u}^{1} \frac{ψ(ỹ)}{\sqrt{y^{2} − u^{2}}} dỹ]

with parameter u defined as

(1.13)  u = r/n(u), 0 ≤ u ≤ 1

and r being the radial position measured from the fiber axis, for simplicity normalized to 1 at its edge (as well as the refractive index).

From the above considerations it follows that the refractive index profile of optical fiber can be quite easily reconstructed by means of formulae (1.9) − (1.12) from four homogeneous field interferograms obtained with π/2 stepped bias Φ_{B}. Such experiment was successfully performed in [1], but we omit the details.

2. A Numerical Approximation Via Taylor’s Expansion

It is obvious that, in practice, an accurate computation of the fiber refractive index profile provided by the analytic formula:

(2.1)  n(u) := exp [\frac{1}{π} ∫_{u}^{1} \frac{ψ(s)}{\sqrt{s^{2} − u^{2}}} ds] , 0 ≤ u ≤ 1

depends on the numerical accuracy in approximating the integral

(2.2)  I(ψ, u) := \frac{1}{π} ∫_{u}^{1} \frac{ψ(s)}{\sqrt{s^{2} − u^{2}}} ds, 0 ≤ u ≤ 1.

A natural approach, if information on the derivatives of ψ are available, is to compute (2.2) by the use of Taylor’s formula

(2.3)  ψ(s) = ψ(u) + ∑_{k=1}^{n} \frac{(s − u)^{k}}{k!} ψ^{(k)}(u) + \frac{1}{n!} ∫_{u}^{s}(s − t)^{n}ψ^{(n+1)}(t)dt

for any s ∈ [u, 1] , u ∈ [0, 1] , where n ≥ 1 is a natural number. Then we get [9]

(2.4)  I(ψ, u) = \frac{1}{π} ψ(u) ∫_{u}^{1} \frac{ds}{\sqrt{s^{2} − u^{2}}}

+\frac{1}{π} ∑_{k=1}^{n} \frac{ψ^{(k)}(u)}{k!} ∫_{u}^{1}(s − u)^{k−1}\sqrt{\frac{s − u}{s + u}}ds

+\frac{1}{πn!} ∫_{u}^{1} \frac{ds}{\sqrt{s^{2} − u^{2}}} (∫_{u}^{s}(s − t)^{n}ψ^{(n+1)}(t)dt)

= \frac{1}{π} M_{0}(u)ψ(u) + \frac{1}{π} ∑_{k=1}^{n} M_{k}(u)ψ^{(k)}(u) + R_{n}(ψ, u)

where

(2.5)  M_{0}(u) : = ∫_{u}^{1} \frac{ds}{\sqrt{s^{2} − u^{2}}} ,
(2.6)  M_{k}(u) : = \frac{1}{k!} ∫_{u}^{1}(s − u)^{k−1}\sqrt{\frac{s − u}{s + u}}ds

and R_{n}(ψ, u) is the remainder in formula (2.4) , i.e.,

R_{n}(ψ, u) := \frac{1}{πn!} ∫_{u}^{1} \frac{ds}{\sqrt{s^{2} − u^{2}}} (∫_{u}^{s}(s − t)^{n}ψ^{(n+1)}(t)dt) .

We note, using Maple 6, that we may compute the functions M_{0}(u), M_{k}(u), for k = 1, 2, ... Since the expressions of M_{k}(u), for k = 2, 3, ... are quite complicate we will not present them here explicitly. We mention only the expressions for M_{0}(u) and M_{1}(u) :

M_{0}(u) = ln (1 + \sqrt{1 − u^{2}}) − ln(u),
M_{1}(u) : = [\sqrt{(1 − u) /(u + 1)}\sqrt{1 − u^{2}}
−\sqrt{(1 − u)/(u + 1)}uln(1 + \sqrt{1 − u^{2}})
+\sqrt{(1 − u)/(u + 1)}u\sqrt{1 − u^{2}}
−\sqrt{(1 − u)/(u + 1)}u^{2}ln(1 + \sqrt{1 − u^{2}})
+\sqrt{1/u}u^{3/2}ln(u)\sqrt{1 − u^{2}}]/(\sqrt{1 − u^{2}}).

The functions M_{k}(u) for k = 2, 3, ...9, which were computed by a Maple program have been used in the computer implementation of the approximation proposed for the refractive index profile n(u) (see [9]).

3. Error Analysis for Taylor’s Expansion

Using the properties of the integral, we have

(3.1)  |R_{n}(ψ, u)| ≤ \frac{1}{πn!} ∫_{u}^{1} |∫_{u}^{s}(s − t)^{n}ψ^{(n+1)}(t)dt| \frac{ds}{\sqrt{s^{2} − u^{2}}}
= : B_{n}(ψ, u).

6 N.M. DRAGOMIR, S.S. DRAGOMIR, Y.J. CHO, AND G.W. BAXTER

If ψ^{(n+1)} is bounded in the interval [0, 1] and by denoting

∥ψ^{(n+1)}∥_{[u,s],∞} := sup_{t∈[u,s]} |ψ^{(n+1)}(t)| < ∞

then we have

(3.2) |∫_{u}^{s}(s − t)^{n}ψ^{(n+1)}(t)dt| ≤ sup_{t∈[u,s]} |ψ^{(n+1)}(t)| ∫_{u}^{s}(s − t)^{n}dt
≤ \frac{(s − u)^{n+1}}{n + 1} ∥ψ^{(n+1)}∥_{[u,1],∞} .

If ψ^{(n+1)} is p-integrable (p > 1) in the interval [0, 1] and if we denote

∥ψ^{(n+1)}∥_{[u,s],p} := (∫_{u}^{s} |ψ^{(n+1)}(t)|^{p} dt)^{\frac{1}{p}} , p ≥ 1;

then by Hölder’s integral inequality for p > 1 (\frac{1}{p} + \frac{1}{q} = 1) , the following relation can be deduced |∫_{u}^{s}(s − t)^{n}ψ^{(n+1)}(t)dt| (3.3)

≤ (∫_{u}^{s}(s − t)^{nq}dt)^{\frac{1}{q}} × (∫_{u}^{s} |ψ^{(n+1)}(t)|^{p} dt)^{\frac{1}{p}}
≤ \frac{(s − u)^{n+ \frac{1}{q}}}{(nq + 1)^{1/q}} ∥ψ^{(n+1)}∥_{[u,1],p} .

If ψ^{(n+1)} is integrable in the interval [0, 1] and if we denote

∥ψ^{(n+1)}∥_{[u,s],1} := ∫_{u}^{s} |ψ^{(n+1)}(t)| dt;

then we have

(3.4) |∫_{u}^{s}(s − t)^{n}ψ^{(n+1)}(t)dt| ≤ sup_{t∈[u,s]} (s − t)^{n} ∫_{u}^{s} |ψ^{(n+1)}(t)| dt
≤ (s − u)^{n} ∥ψ^{(n+1)}∥_{[u,1],1} .

Using (3.2) − (3.4) the upper bounds for B_{n}(ψ, u), can be written as, see [9], B_{n}(ψ, u) (3.5)

≤ \frac{1}{π} × \{ \begin{matrix}\frac{1}{(n+1)!} ∫_{u}^{1}(s − u)^{n}\sqrt{\frac{s−u}{s+u}} ∥ψ^{(n+1)}∥_{[u,s],∞} dt \\ \frac{1}{n!(nq+1)^{1/q}} ∫_{u}^{1}(s − u)^{n+ \frac{1}{q} −1}\sqrt{\frac{s−u}{s+u}} ∥ψ^{(n+1)}∥_{[u,s],p} dt \\ \frac{1}{n!} ∫_{u}^{1}(s − u)^{n−1}\sqrt{\frac{s−u}{s+u}} ∥ψ^{(n+1)}∥_{[u,s],1} dt \end{matrix}

≤ \frac{1}{π} × \{ \begin{matrix}\frac{1}{(n+1)!} ∥ψ^{(n+1)}∥_{[u,1],∞} ∫_{u}^{1}(s − u)^{n}\sqrt{\frac{s−u}{s+u}}dt \\ \frac{1}{n!(nq+1)^{1/q}} ∥ψ^{(n+1)}∥_{[u,1],p} ∫_{u}^{1}(s − u)^{n+ \frac{1}{q} −1}\sqrt{\frac{s−u}{s+u}}dt \\ \frac{1}{n!} ∥ψ^{(n+1)}∥_{[u,1],1} ∫_{u}^{1}(s − u)^{n−1}\sqrt{\frac{s−u}{s+u}}dt \end{matrix} .

Since, the following type of integrals occur

N(α, u) := ∫_{u}^{1}(s − u)^{α}\sqrt{\frac{s − u}{s + u}}dt, α ≥ 0, u ∈ [0, 1]

and these type of integrals can not be exactly computed for general α (note, however, that for particular small α Maple 6 can provide the exact value), some evaluations of N(α, u) are highly desirable.

Having this in mind, we observe that, the simplest way to obtain an upper bound for N(α, u), is to consider that

0 ≤ \frac{s − u}{s + u} ≤ 1, for u ∈ [0, 1] , u ≤ s ≤ 1,

giving

N(α, u) ≤ ∫_{u}^{1}(s − u)^{α}dt = \frac{(1 − u)^{α+1}}{(α + 1)} , α ≥ 0, u ∈ [0, 1] .

Using (3.4) and (3.5) we get the following simple estimate for the absolute value of the remainder R_{n}(ψ, u) in approximating the integral I(u) with the analytic expression [9]:

(3.6)  A_{n}(ψ, u) := \frac{1}{π} M_{0}(u)ψ(u) + \frac{1}{π} ∑_{k=1}^{n} M_{k}(u)ψ^{(k)}(u)

where M_{0}(u), M_{k}(u) were defined in the previous section,

(3.7)  |R_{n}(ψ, u)| ≤ \frac{1}{π} × \{ \begin{matrix}\frac{(1−u)^{n+1}}{(n+1)!(n+1)} & ∥ψ^{(n+1)}∥ & [u,1],∞ \\ \frac{(1−u)^{n+ \frac{1}{q}}}{n!(nq+1)^{1/q+1}} & ∥ψ^{(n+1)}∥ & [u,1],p \\ \frac{(1−u)^{n}}{n!n} ∥ψ^{(n+1)}∥_{[u,1],1} . \end{matrix}

The equation (3.7) shows that the remainder R_{n}(ψ, u) is (rapidly) uniformly convergent to 0 as n → ∞, meaning that the approximation of I(u) by A_{n}(ψ, u) is accurate for large enough n.

4. A Numerical Approximation Via Trapezoidal Rule

As noted before, an accurate computation of the fiber refractive index profile provided by the analytic formula

(4.1)  n (u) := exp [\frac{1}{π} ∫_{u}^{1} \begin{matrix}ψ & (s) \\ √ & \overline{s^{2} − u^{2}} \end{matrix} ds] , 0 ≤ u ≤ 1

depends on the numerical accuracy in approximating the integral

(4.2)  I (u) := \frac{1}{π} ∫_{u}^{1} \begin{matrix}ψ & (s) \\ √ & \overline{s^{2} − u^{2}} \end{matrix} ds, 0 ≤ u ≤ 1.

If we change the variable such that τ = \sqrt{s^{2} − u^{2}}, s ∈ [u, 1] , then τ = 0 for s = u and τ = \sqrt{1 − u^{2}} for s = 1, s = \sqrt{τ^{2} + u^{2}} and

ds = \frac{τdτ}{\sqrt{τ^{2} + u^{2}}} ,

and I (u) becomes [10]

(4.3)  I (u) = \frac{1}{π} ∫_{0}^{√\overline{1−u^{2}}} \begin{matrix}ψ & (\sqrt{τ^{2} + u^{2}}) \\ √ & \overline{τ^{2} + u^{2}} \end{matrix}dτ, 0 ≤ u ≤ 1.

This integral may be approximated in many different ways by the use of classical quadrature rules including the trapezoidal rule, mid-point rule, Simpson’s rule etc...

It is well known that [2], if f : [a, b] → R is a twice differentiable function and the second derivative f′′ is bounded, this means that ∥f′′∥_{∞} := sup_{t∈[a,b]} |f′′(t)| < ∞, then for any division

I_{n} : a = x_{0} < x_{1} < ... < x_{n−1} < x_{n} = b

of the interval [a, b] with h_{i} := x_{i+1} − x_{i}(i = 0, ..., n − 1) one has

(4.4)  ∫_{a}^{b} f (x) dx = T_{n} (f; I_{n}) + R_{n} (f; I_{n})

where

T_{n} (f; I_{n}) := ∑_{i=0}^{n−1} \frac{f (x_{i}) + f (x_{i+1})}{2} · h_{i}

is the trapezoidal rule, and the error of approximating the integral in this way, denoted in (4.4) by R_{n} (f; I_{n}) , satisfies the estimate

(4.5)  |R_{n} (f; I_{n})| ≤ \frac{1}{12} ∥f′′∥_{∞} ∑_{i=0}^{n−1} h_{i}^{3}

≤ \frac{1}{12} ∥f′′∥_{∞} [ν (h)]^{2} (b − a) ,

where ν (h) := max {h_{i} := x_{i+1} − x_{i}; i = 0, ..., n − 1} is the norm of the division I_{n}. If ν (h) → 0 (as n → ∞) , then by (4.5) R_{n} (f; I_{n}) → 0 with order two of accuracy.

Now, for a given u ∈ (0, 1) , consider the division of the interval [0, \sqrt{1 − u^{2}}] given by

Δ_{n} : δ_{i} = λ_{i}\sqrt{1 − u^{2}}, 0 = λ_{0} < λ_{1} < ... < λ_{n−1} < λ_{n} = 1.

Using the trapezoidal rule, we can approximate I (u) in the following manner [10]

(4.6)  I (u) = A_{n} (ψ, λ, u) + W_{n} (ψ, λ, u) , u ∈ [0, 1] ;

where

A_{n} (ψ, λ, u) (4.7)

: = \frac{\sqrt{1 − u^{2}}}{2π} ∑_{i=0}^{n−1} (λ_{i+1} − λ_{i})

× [\begin{matrix}ψ & (\sqrt{(1 − u^{2}) λ_{i}^{2} + u^{2}}) \\ √ & \overline{(1 − u^{2}) λ_{i}^{2} + u^{2}} \end{matrix} + \begin{matrix}ψ & (\sqrt{(1 − u^{2}) λ_{i+1}^{2} + u^{2}}) \\ √ & \overline{(1 − u^{2}) λ_{i+1}^{2} + u^{2}} \end{matrix} ]

and W_{n} (ψ, λ, u) is the remainder in approximating I (u) by A_{n} (ψ, λ, u) .

5. Using Data Provided by Practical Experiments

In practice, performing accurate laboratory experiments, the function ψ may be measured in a number of (m + 1)-points, say

0 = s_{0} < s_{1} < s_{2} < ... < s_{m−1} < s_{m} = 1.

We will show now how the approximate value provided by formula (4.7) for the integral I(u) may be computed by the use of the assumed known values ψ (0) , ψ (s_{1}) , ..., ψ (s_{m−1}) , and ψ (1) .

If u ∈ (0, s_{1}], we may choose the division

λ_{0} = 0, λ_{1} = \sqrt{\frac{s_{1}^{2} − u^{2}}{1 − u^{2}}} , λ_{2} = \sqrt{\frac{s_{2}^{2} − u^{2}}{1 − u^{2}}} , ...,

λ_{m−1} = \sqrt{\frac{s_{m−1}^{2} − u^{2}}{1 − u^{2}}} , λ_{m} = 1;

which obviously satisfies the condition

0 = λ_{0} < λ_{1} < ... < λ_{m−1} < λ_{m} = 1.

Then by formula (4.7) , for n = m, we get

(5.1) A_{n} (ψ, λ, u) = \frac{1}{2π} \{\sqrt{s_{1}^{2} − u^{2}} [\frac{ψ (u)}{u} + \frac{ψ (s_{1})}{s_{1}}]

+ ∑_{i=1}^{m−2} (\sqrt{s_{i+1}^{2} − u^{2}} − \sqrt{s_{i}^{2} − u^{2}}) [\frac{ψ (s_{i+1})}{s_{i+1}} + \frac{ψ (s_{i})}{s_{i}}]

+ (\sqrt{1 − u^{2}} − \sqrt{s_{m−1}^{2} − u^{2}}) [\frac{ψ (s_{m−1})}{s_{m−1}} + ψ (1)]\} .

Since u ∈ (0, s_{1}], instead of \frac{ψ(u)}{u} in formula (5.1) we may choose either lim_{u→0+} \frac{ψ(u)}{u} or \frac{ψ(s_{1})}{s_{1}}.

For u ∈ (s_{1}, s_{2}], we may choose the division

λ_{0} = 0, λ_{1} = \sqrt{\frac{s_{2}^{2} − u^{2}}{1 − u^{2}}} , λ_{2} = \sqrt{\frac{s_{3}^{2} − u^{2}}{1 − u^{2}}} , ...,

λ_{m−2} = \sqrt{\frac{s_{m−1}^{2} − u^{2}}{1 − u^{2}}} , λ_{m−1} = 1;

which obviously satisfies

0 = λ_{0} < λ_{1} < ... < λ_{m−2} < λ_{m−1} = 1

and for n = m − 1, formula (4.7) will produce the quadrature [10] A_{m−1} (ψ, λ, u) (5.2)

= \frac{1}{2π} \{\sqrt{s_{2}^{2} − u^{2}} [\frac{ψ (u)}{u} + \frac{ψ (s_{2})}{s_{2}}]

+ ∑_{i=1}^{m−3} (\sqrt{s_{i+2}^{2} − u^{2}} − \sqrt{s_{i+1}^{2} − u^{2}}) [\frac{ψ (s_{i+2})}{s_{i+2}} + \frac{ψ (s_{i+1})}{s_{i+1}}]

+ (\sqrt{1 − u^{2}} − \sqrt{s_{m−1}^{2} − u^{2}}) [\frac{ψ (s_{m−1})}{s_{m−1}} + ψ (1)]\} .

Here we may use instead of \frac{ψ(u)}{u} in formula (5.2) the value \frac{ψ(s_{1})}{s_{1}}.

In general, if we assume that u ∈ (s_{k−1}, s_{k}], k = 1, ..., m − 1; we may choose the division

λ_{0} = 0, λ_{1} = \sqrt{\frac{s_{k}^{2} − u^{2}}{1 − u^{2}}} , λ_{2} = \sqrt{\frac{s_{k+1}^{2} − u^{2}}{1 − u^{2}}} , ...,

λ_{m−k} = \sqrt{\frac{s_{m−1}^{2} − u^{2}}{1 − u^{2}}} , λ_{m−k+1} = 1;

which satisfies the condition

0 = λ_{0} < λ_{1} < ... < λ_{m−k} < λ_{m−k+1} = 1

and for which, by applying (4.7) we get the quadrature (n = m − k + 1) [10]: A_{m−k+1} (ψ, λ, u) (5.3)

= \frac{1}{2π} \{\sqrt{s_{k}^{2} − u^{2}} [\frac{ψ (u)}{u} + \frac{ψ (s_{k})}{s_{k}}]

+ ∑_{i=1}^{m−k−1} (\sqrt{s_{k+i}^{2} − u^{2}} − \sqrt{s_{k+i−1}^{2} − u^{2}}) [\frac{ψ (s_{k+i})}{s_{k+i}} + \frac{ψ (s_{k+i−1})}{s_{k+i−1}}]

+ (\sqrt{1 − u^{2}} − \sqrt{s_{m−1}^{2} − u^{2}}) [\frac{ψ (s_{m−1})}{s_{m−1}} + ψ (1)]\} .

Here, instead of the term \frac{ψ(u)}{u} we may take \frac{ψ(s_{k−1})}{s_{k−1}}.

Now, if u ∈ (s_{m−2}, s_{m−1}], then we may choose the division

λ_{0} = 0, λ_{1} = \sqrt{\frac{s_{m−1}^{2} − u^{2}}{1 − u^{2}}} , λ_{2} = 1

and the approximation formula will be (by (4.7))

(5.4) A_{3} (ψ, λ, u) = \frac{1}{2π} \{\sqrt{s_{m−1}^{2} − u^{2}} [\frac{ψ (u)}{u} + \frac{ψ (s_{m−1})}{s_{m−1}}]

+ (\sqrt{1 − u^{2}} − \sqrt{s_{m−1}^{2} − u^{2}}) [\frac{ψ (s_{m−1})}{s_{m−1}} + ψ (1)]\}

and the term \frac{ψ(u)}{u} may be replaced by \frac{ψ(s_{m−2})}{s_{m−2}}.

Finally, for the last interval (s_{m−1}, 1], we may take the division λ_{0} = 0, λ_{1} = 1 and for n = 2, we have the quadrature

(5.5)  A_{2} (ψ, λ, u) = \frac{\sqrt{1 − u^{2}}}{2π} [\frac{ψ (u)}{u} + ψ (1)]

in which \frac{ψ(u)}{u} may be replaced by \frac{ψ(s_{m−1})}{s_{m−1}}.

To numerically implement the above quadrature rule, an equidistant partitioning of the interval [0, 1] is natural to be considered; details of this implementations may be found in [10].

6. New Approximation Using Stieltjes Integrals

Let f, g : [a, b] → R be two mappings such that the Stieltjes integral ∫_{a}^{b} f (t) dg (t) exists. It is well known that if, for instance f is continuous and g is of bounded variation, then the above Stieltjes integral exists and

(6.1)  |∫_{a}^{b} f (t) dg (t)| ≤ max_{t∈[a,b]} |f (t)| ⋁_{a}^{b} (g) ,

where ⋁_{a}^{b} (g) denotes the total variation of g. We also note that if f is Riemann

integrable on [a, b] and g is Lipschitzian on the same interval, then ∫_{a}^{b} f (t) dg (t) exists and

(6.2)  |∫_{a}^{b} f (t) dg (t)| ≤ L ∫_{a}^{b} |f (t)| dt,

where L is the Lipschitzian constant of g on [a, b] .

Now, consider I_{n} : a = x_{0} < x_{1} < ... < x_{n−1} < x_{n} = b to be a division of the interval [a, b] , h_{i} := x_{i+1}−x_{i} (i = 0, ..., n − 1) and, v (h) := max {h_{i}|i = 0, ..., n − 1} . Define the general Riemann-Stieltjes sum

(6.3)  S (f, g, I_{n}, ξ) := ∑_{i=0}^{n−1} f (ξ_{i}) [g (x_{i+1}) − g (x_{i})] ,

where ξ = (ξ_{0}, ..., ξ_{n−1}) is a sequence of intermediate points in I_{n}, i.e. ξ_{i} ∈ [x_{i}, x_{i+1}] , i = 0, ..., n − 1.

If we consider the approximation

(6.4)  ∫_{a}^{b} f (t) dg (t) = S (f, g, I_{n}, ξ) + R (f, g, I_{n}, ξ) ,

where S (f, g, I_{n}, ξ) is the Riemann-Stieltjes sum given by (6.3) and R (f, g, I_{n}, ξ) is the remainder, then it has been shown in [11] that for f Lipschitzian with the constant L and g of bounded variation, the error estimate satisfies the bound

(6.5)  |R (f, g, I_{n}, ξ)| ≤ L [\frac{1}{2}v (h) + max_{i=0,...,n−1} |ξ_{i} − \frac{x_{i} + x_{i+i}}{2} |] ⋁_{a}^{b} (g) .

It is obvious that the best inequality one can have in (6.5) is for ξ_{i} = \frac{x_{i}+x_{i+1}}{2} ; i = 0, ..., n − 1 obtaining the approximation

(6.6)  ∫_{a}^{b} f (t) dg (t) = Y (f, g, I_{n}) + R (f, g, I_{n})

where Y (f, g, I_{n}) is the mid-point formula:

(6.7)  Y (f, g, I_{n}) := ∑_{i=0}^{n−1} f (\frac{x_{i} + x_{i+1}}{2} ) [g (x_{i+1}) − g (x_{i})]

and the remainder satisfies the estimate:

(6.8)  R (f, g, I_{n}) ≤ \frac{1}{2}Lv (h) ⋁_{a}^{b} (g) .

Here the constant \frac{1}{2} is best possible in the sense that it cannot be replaced by a smaller constant.

Now, let consider the integral

(6.9)  I (Ψ, u) := \frac{1}{π} ∫_{u}^{1} \frac{Ψ (s)}{\sqrt{s^{2} − u^{2}}} ds, 0 ≤ u ≤ 1

which is connected to the fiber refractive index profile by

(6.10)  n (u) = exp [I (Ψ, u)] , 0 ≤ u ≤ 1.

One may observe that

(6.11)  I (Ψ, u) = \frac{1}{π} ∫_{u}^{1} f (s) dg (s) ,

where f (s) = \frac{Ψ(s)}{s} and g (s) = \sqrt{s^{2} − u^{2}}, s ∈ [u, 1] , 0 < u ≤ 1.

Now is we assume that Ψ is differentiable, then

(6.12)  f′ (s) = \frac{sΨ′ (s) − Ψ (s)}{s^{2}}, s ∈ [u, 1] , 0 < u ≤ 1

and if there exists a constant M > 0 such that

(6.13)  |sΨ′ (s) − Ψ (s)| ≤ Ms^{2}, s ∈ [u, 1] , 0 < u ≤ 1,

then we may conclude that f is Lipschitzian with the constant M mentioned above. Since g′ (s) = \frac{s}{\sqrt{s^{2}−u^{2}}} , s ∈ [u, 1] , 0 < u ≤ 1, then

(6.14)  ⋁_{u}^{1} (g) = ∫_{u}^{1} |g′ (s)| ds = ∫_{u}^{1} \frac{sds}{\sqrt{s^{2} − u^{2}}} = \sqrt{1 − u^{2}}.

For a given u ∈ (0, 1) , consider now the division

(6.15)  (A_{n}) : u = δ_{0} < δ_{1} < δ_{2} < ... < δ_{n−1} < δ_{n} = 1.

For this division and the above choices of the function f and g, we may consider the mid-point rule

(6.16) Y (Ψ, u, A_{n}) = ∑_{i=0}^{n−1} \frac{Ψ (\frac{δ_{i}+δ_{i+1}}{2} )}{\frac{δ_{i}+δ_{i+1}}{2}} · [\sqrt{δ_{i+1}^{2} − u^{2}} − \sqrt{δ_{i}^{2} − u^{2}}]

= 2 ∑_{i=0}^{n−1} Ψ (\frac{δ_{i} + δ_{i+1}}{2} ) · \frac{δ_{i+1} − δ_{i}}{\sqrt{δ_{i+1}^{2} − u^{2}} + \sqrt{δ_{i}^{2} − u^{2}}} .

Using the approximation (6.6) , we may conclude that the desired integral I (Ψ, u) may be approximated in the following way

(6.17)  I (Ψ, u) = \frac{1}{π} Y (Ψ, u, A_{n}) + Q (Ψ, u, A_{n}) ,

where the error Q (Ψ, u, A_{n}) , by (6.8), satisfies the estimate

(6.18)  |Q (Ψ, u, A_{n})| ≤ \frac{1}{2π} Mv (δ) \sqrt{1 − u^{2}}, 0 < u < 1.

Here M is the constant defined in (6.13) and v (δ) = max_{i=0,...,n−1} (δ_{i+1} − δ_{i}) .

In practical situations, it is useful if we consider the equidistant partitioning

δ_{i} := u + i · \frac{1 − u}{n}, i = 0, ..., n.

In this case the quadrature (6.16) becomes

Y_{n} (Ψ, u) (6.19)

= \frac{2 (1 − u)}{n} · ∑_{i=0}^{n−1} Ψ [u + (i + \frac{1}{2}) · \frac{1 − u}{n}]

×\frac{1}{\sqrt{\frac{(i+1)^{2}(1−u)^{2}}{n^{2}} + 2 (i + 1) u\frac{1−u}{n}} + \sqrt{\frac{i^{2}(1−u)^{2}}{n^{2}} + 2iu\frac{1−u}{n}}} .

Consequently we have the approximation formula

(6.20)  I (Ψ, u) = \frac{1}{π} Y_{n} (Ψ, u) + Q_{n} (Ψ, u) ,

where Y_{n} (Ψ, u) is defined by (6.19) and the error of approximation Q_{n} (Ψ, u) satisfies the estimate

(6.21)  |Q_{n} (Ψ, u)| ≤ \frac{1 − u}{2πn}M\sqrt{1 − u^{2}} ≤ \frac{M}{2πn}, 0 < u < 1, n ≥ 1.

It is obvious that Q_{n} (Ψ, u) → 0 uniformly over u ∈ [0, 1] as n → ∞.

The above inequality (6.21) provides an a priori estimate for the error in approximating I (Ψ, u) by Y_{n} (Ψ, u) , u ∈ (0, 1) when the constant M is known.

For the division I_{n} : a = x_{0} < x_{1} < ... < x_{n−1} < x_{n} = b of the interval [a, b] , consider the trapezoidal rule for Stieltjes integral

(6.22)  T (f, g, I_{n}) := ∑_{i=0}^{n−1} \frac{f (x) + f (x_{i+1})}{2} [g (x_{i+1}) − g (x_{i})] .

It has been shown that for appropriate choices of f and g one can approximate the Stieltjes integral ∫_{a}^{b} f (t) dg (t) in the following way [12, Chapter 8]

(6.23)  ∫_{a}^{b} f (t) dg (t) = T (f, g, I_{n}) + W (f, g, I_{n}) ,

where W (f, g, I_{n}) denotes the error of the approximation.

If f, g : [a, b] → R is Lipschitzian with the constant L and u : [a, b] → R is a bounded variation on [a, b] then (see [12, Chapter 8])

(6.24)  |W (f, g, I_{n})| ≤ \frac{1}{2}Lv (h) ⋁_{a}^{b} (u) ,

where v (h) := max_{i=0,...,n−1} {h_{i}} , and h_{i} := x_{i+1} − x_{i}, i ∈ {0, ..., n − 1} .

Consider now the integral from (6.11)

(6.25)  I (Ψ, u) = \frac{1}{π} ∫_{u}^{1} f (s) dg (s) ,

where f (s) = \frac{Ψ(s)}{s}, g (s) = \sqrt{s^{2} − u^{2}}, s ∈ [u, 1] , 0 < u ≤ 1, and assume that there exists a constant M > 0 such that (6.13) holds. Further, assume that u ∈ (0, 1) and utilizing the division of the interval [0, 1] as defined by (6.15) for the functions f and g, one may consider the trapezoidal rule

T (Ψ, u, A_{n}) (6.26)

= ∑_{i=0}^{n−1} \frac{\frac{Ψ(δ_{i+1})}{δ_{i+1}} + \frac{Ψ(δ_{i})}{δ_{i}}}{2} [\sqrt{δ_{i+1}^{2} − u^{2}} − \sqrt{δ_{i}^{2} − u^{2}}]

= \frac{1}{2} ∑_{i=0}^{n−1} \frac{δ_{i}Ψ (δ_{i+1}) + δ_{i+1}Ψ (δ_{i})}{δ_{i}δ_{i+1}} [\sqrt{δ_{i+1}^{2} − u^{2}} − \sqrt{δ_{i}^{2} − u^{2}}] .

Using the approximation (6.23) , the integral in (6.25) may be computed in the following way

(6.27)  I (Ψ, u) = \frac{1}{π} T (Ψ, u, A_{n}) + S (Ψ, u, A_{n}) ,

where the error S (Ψ, u, A_{n}) defined by (6.24) satisfies the estimate

(6.28)  |S (Ψ, u, A_{n})| ≤ \frac{1}{2π} Mv (δ) \sqrt{1 − u^{2}}, 0 < u < 1.

The constant M in (6.28) is defined by (6.13) and v (δ) = max_{i=0,...,n−1} (δ_{i+1} − δ_{i}) .

Usually an equidistant partitioning of the interval [u, 1] is preferred, i.e. δ_{i} = u + i · \frac{1−u}{n}, i = 0, ..., n. As such, the quadrature T (Ψ, u, A_{n}) becomes

T_{n} (Ψ, u) (6.29)

= \frac{1}{2} ∑_{i=0}^{n−1} \frac{(u + i · \frac{1−u}{n}) Ψ [u + (i + 1) \frac{1−u}{n}] + [u + (i + 1) \frac{1−u}{n}] Ψ (u + i · \frac{1−u}{n})}{(u + i · \frac{1−u}{n}) [u + (i + 1) · \frac{1−u}{n}]}

× [\sqrt{\frac{(i + 1)^{2} (1 − u)^{2}}{n^{2}} + 2 (i + 1) u\frac{1 − u}{n}} − \sqrt{\frac{i^{2} (1 − u)^{2}}{n^{2}} + 2iu\frac{1 − u}{n}}] .

In conclusion, we have the approximation formula

(6.30)  I (Ψ, u) = \frac{1}{π} T_{n} (Ψ, u) + S_{n} (Ψ, u) ,

where T_{n} (Ψ, u) is defined in (6.29) , with an error of approximation S_{n} (Ψ, u) satisfying the bound

|S_{n} (Ψ, u)| ≤ \frac{1 − u}{2πn}M\sqrt{1 − u^{2}} ≤ \frac{M}{2πn}, 0 < u < 1, n ≥ 1.

Also note that S_{n} (Ψ, u) → 0 uniformly over u ∈ [0, 1] as n → ∞.

7. Numerical experiments

Consider the function Ψ : [0, 1] → [0, ∞), Ψ (s) = s^{3}. The exact transform of Ψ is given by

(7.1)  I (Ψ, u) = (1/3\sqrt{(1 − u^{2})} + 2/3\sqrt{(1 − u^{2})}u^{2})/π, u ∈ [0, 1] .

The plot of this function is incorporated in Figure 1. The plots of the errors Q_{n} (Ψ, u) (for Stieltjes mid-point rule) and S_{n} (Ψ, u) (for Stieltjes trapezoidal rule)

Figure 1. Exact value for I (Ψ, u) , where Ψ (s) = s^{3}.

Figure 2. Comparisson of the errors as computed by the midpoint rule and trapezoidal rule for n = 10.

computed by Maple 6, for n = 10, 10^{2} and 10^{3} are incorporated in Figure 2, 3 and 4, respectively.

Figure 3. Comparisson of the errors as computed by the midpoint rule and trapezoidal rule for n = 10^{2}.

Figure 4. Comparisson of the errors as computed by the midpoint rule and trapezoidal rule for n = 10^{3}.

Acknowledgement: S. S. Dragomir and Y. J. Cho greatly acknowledge the financial support from Brain Pool Program (2002) of the Korean Federation of Science and Technology Societies. The research was performed under the Memorandum of Understanding between Victoria University and Gyeongsang National University.

References

[1] J. Sochacki, Accurate reconstruction of the refractive-index profile of fibers and preform rods from transverse interferometric data, Applied Optics, 25(19), 3473-3482, 1986.

[2] K.E. Atkinson, An Introduction to Numerical Analysis, Wiley and Sons, Second Edition, 1989.

[3] M. Sochacka, Optical fiber profiling by phase stepping transverse interferometry, SPIE, Vol. 1846, 160-175, 1994.

[4] J. Sochacki, D. Rogus and M. Sochacka, Reliability in reconstruction of the axisymmetric refractive index profiles from transverse interferograms. Accurate versus approximate methods: a comparative study, Fiber and Integrated Optics, 8(4), 279-307, 1987.

[5] J.H. Brunning, D.R. Herriott, J.E. Gallagher, D.P. Rosenfeld, A.D. White, D.J. Brangaccio, Digital wavefront measuring interferometer for testing surfaces and lenses, Applied Optics 13, 2693-2703, 1974.

[6] M. Kujawinska, Automatic Fringe Pattern Analysis in Optical Methods of Testing, in ”Prace Naukowe-Mechanika”, Warsaw Technical University, p.138, 1990.

[7] K. Creath, Phase-measurement interferometry techniques, in Progress in Optics ,Vol. XXVI, (Ed. E. Wolf), 349-393, 1988.

[8] M Pluta, Advanced Light Microscopy, (Polish Scientific Publishers), Elsevier, Amsterdam, Vol. 2, Chap. 7.1.3, 1989.

[9] N. M, Dragomir, and G. W. Baxter, An accurate approximation for the fiber refractive index profile via Taylor’s expansion formula, J. KSIAM, 6 (2), 43-52, 2002.

[10] N. M. Dragomir, and G. W. Baxter, An approximation for the fiber refractive index profile via trapezoid formula, Tamsui Oxford Journal of Mathematical Sciences, 18, 193-196, 2002.

[11] S. S. Dragomir, On the Ostrowski inequality for Riemann-Stieltjes integral ∫_{a}^{b} f (t) d (u (t)) , where f is of Hölder type and u is of bounded variation and application, J. KSIAM, 5(1), 35-45, 2001.

[12] S. S. Dragomir and Th. M. Rassias, (Eds), Ostrowski Type Inequalities and Applications in Numerical Integration, Kluwer Academic Publishers, Dordrecht, 2002.

School of Electrical Engineering, Victoria University of Technology, P.O. Box 14428, Melbourne City, MC 8001, Melbourne, Australia.

E-mail address: nicoleta.dragomir@vu.edu.au

School of Computer Science & Mathematics, Victoria University of Technology, P.O. Box 14428, Melbourne City, MC 8001, Melbourne, Australia.

E-mail address: sever.dragomir@vu.edu.au

Department of Mathematics, College of Education, Gyeongsang National University, Chinju, 660 - 701, Korea.

E-mail address: yjcho@nongae.gsnu.ac.kr

School of Electrical Engineering, Victoria University of Technology, P.O. Box 14428, Melbourne City, MC 8001, Melbourne, Australia.

E-mail address: gregory.baxter@vu.edu.au