Abstract

In this paper, the smooth approximation of light waves is studied for an open optical waveguide with a distinct refractive-index profile, which involves high-precision computation of the eigenmodes and corresponding eigenfunctions. During analysis, the refractive-index function is first approximated by a quadratic spline interpolation function. Since the quadratic spline function is a polynomial of degree two in every sub-interval (sub-layer), it is equivalent to a piecewise polynomial of degree two, based on which, the corresponding Sturm-Liouville eigenvalue problem of the Helmholtz operator in sub-layer can be solved analytically by the Kummer functions. Finally, the approximate dispersion equation is established to the TE case. Obviously, the approximate dispersion equations converge fast to the exact ones, as the maximum value of the sub-interval sizes tends to zero. Furthermore, eigenmodes may be obtained by Müller’s method with suitable initial values. To refine the accuracy, the equidistant partition and the non-equidistant partition are applied to divide the interval. Numerical simulations show that the eigenfunctions of the spline interpolation are much smoother than the ones with piecewise interpolation. In addition, the non-equidistant partition can help improve the accuracy and the order of convergence of general solutions reaches the third.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

The propagation characteristics (called as propagation constants) of optical waveguides play an important role in the optimal designs of optical devices, including fiber [1,2] and the computation of wave propagation with the eigenmode expansion method. Furthermore, the wave field can be decomposed into a linear combination of a few discrete guided modes and an infinite integral of a continuum of radiation modes [3]. The infinite integral in the decomposition of a field can be approximately expressed as an infinite sum of leaky modes through discretization [4,5], which requires to develop some efficient methods to obtain high-precision propagation constants [6,7]. In general, numerical methods, such as the finite difference method [8,9], the finite element method [10], the multi-domain pseudo-spectral method [11,12], and the B-spline modal method [13], can help simplify the problem into a matrix eigenvalue problem. However, for the reason that these methods all arise complex matrices, it is too complicated to calculate the related eigenvalues, let alone satisfying high precision requirements. Now, two problems are encountered. For one thing, the propagation constants of the guided modes are easier to get by common methods compared with the unreal or unbounded ones, since they are real and bounded [14]. For another, the propagation constants of the leaky modes are complex numbers [15], implying the computation difficulty, especially those with large modules. To solve above tasks, an alternative approach is proposed, that is to transform the eigen problem into a nonlinear dispersion equation, whose roots correspond to the eigenvalues or the propagation constants. Some methods can be adopted to deal with the piecewise constant slab waveguides [1416] and rib waveguides [17]. For the varying refractive-index waveguides, part of approximate dispersion equations are constructed by the Wentzel-Kramers-Brillouin (WKB) method [18] and the differential transfer matrix method (DTMM) [1922]. However, these approximations cannot guarantee the modes’ accuracy. Recently, some efficient methods are provided to compute the propagation constants, but their corresponding eigenfunctions are not smooth [23], but continuous, which changes the physical meaning of the line.

In this article, we consider the eigenmode problems of varying refractive-index waveguides in TE case. First, we approximate the refractive-index function by a quadratic spline interpolation function. Second, since the quadratic spline function is a polynomial of degree two in every sub-interval (sub-layer), the corresponding Sturm-Liouville eigenvalue problem of the Helmholtz operator in sub-layer can be solved analytically by the Kummer functions. Finally, the approximate dispersion equation is established to the TE case. The obtained dispersion equations can be solved by the secant method or Muller’s method. Numerical simulations show that our method reaches the third order accuracy. The error estimate originally proved for the real eigenvalues in [24], is still valid for the complex eigenvalues. The rest of this paper is organized as follows. In Section 2, the mathematical model and treatment are introduced, and the dispersion equation is established, where asymptotic solutions of the propagation constants (corresponding to leaky modes) are provided for TE case. In Section 3, numerical results and corresponding analysis are presented. Finally, the conclusions are given in Section 4.

2. Mathematical modelling

In this section, we consider a general two-dimensional waveguide with a refractive-index profile ${n_0}(x)$ shown in Fig. 1. For simplicity, we assume that the medium is homogeneous for $x < 0$ and $x > d$ with d being a constant. That is

$$n = \left\{ {\begin{array}{lc} {{n_1}},& {x < 0},\\ {{n_0}(x)},&{0 \le x \le d,}\\ {{n_2}},& {x > d,} \end{array}} \right.$$
where ${n_0}(x)$ is a continuous function with ${n_0}(x) > \max ({n_1},{n_2})$, and $x = 0$ and $x = d$ are two discontinuous points. The Sturm-Liouville eigenvalue equation is
$$\left\{ \begin{array}{l} \mathrm{\rho} \frac{\rm d }{{{\rm d} x}}\left( {\frac{1}{\rho }\frac{{\rm{d} \mathrm{\phi} }}{{{\rm d} x}}} \right) + \mathrm{\kappa}_0^2n_0^2(x)\mathrm{\phi} = {\mathrm{\beta}^2}\mathrm{\phi} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} - \infty < x < \infty ,\\ \mathop {\lim }\limits_{|x |\to + \infty } \mathrm{\phi} = 0. \end{array} \right.$$
with $\mathrm{\rho} = 1$ for TE case and $\mathrm{\rho} = n_0^2(x)$ for TM case, where ${\mathrm{\kappa} _0}$ is the free space wavenumber, ${\mathrm{\beta} ^2}$ is an eigenvalue, and $\mathrm{\phi} (x)$ is an eigenfunction corresponding to ${\mathrm{\beta} ^2}$. For simplification, we only study the TE case

 

Fig. 1. Subdivision diagram of a three-layers’ waveguide (Ref. [23], Fig. 1), where $d = k \cdot h, k = 3$, and $h$ is a step width.

Download Full Size | PPT Slide | PDF

In the TE case, Eq. (1) in $[0,d]$ is given by

$$\frac{{{\rm{d} ^2}\mathrm{\phi} }}{{{\rm d} {x^2}}} + \mathrm{\kappa} _0^2n_0^2(x)\mathrm{\phi} = {\mathrm{\beta} ^2}\mathrm{\phi} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} 0 < x < d,$$
$$\frac{{\rm{d} \mathrm{\phi} }}{{{\rm d} x}} ={-} i \sqrt {\mathrm{\kappa} _0^2n_1^2 - {\mathrm{\beta} ^2}} \mathrm{\phi} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} x = {0^ + },$$
$$\frac{{\rm{d} \mathrm{\phi} }}{{{\rm d} x}} = i \sqrt {\mathrm{\kappa} _0^2n_2^2 - {\mathrm{\beta} ^2}} \mathrm{\phi} ,{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} \begin{array}{{cc}} {}&{} \end{array}{\kern 1pt} x = {\rm{d}^ - }{\kern 1pt} ,{\kern 1pt} {\kern 1pt} {\kern 1pt}$$
where $i = \sqrt { - 1}$.

First, as $x \in [{x_j},{x_{j + 1}}] \subset [0,d]$, the function $\mathrm{\kappa} _0^2n_0^2(x)$ is interpolated by a polynomial of degree two, that is $\mathrm{\kappa} _0^2n_0^2(x)\textrm{ - }{{\beta} ^2} \approx {a_j}{x^2} + {b_j}x + {c_j}$. Denote ${y_j}(x)$ as the approximation of the $\mathrm{\phi} (x)$. Thus, the form of approximated equation is as following form:

$$\frac{{{\textrm{\rm{d}}^2}{y_j}}}{{\textrm{\rm{d}}{x^2}}} + ({a_j}{x^2} + {b_j}x + {c_j}){y_j} = 0,{\kern 1pt}$$
where $(j = 1,2, \cdot{\cdot} \cdot ,k)$.

2.1 Interval partitions

We adopt following two partitions as follows.

2.1.1 Equidistant partition

The interval $[0,d]$ is divided into k subintervals with width $h = d / k$, then the subintervals are denoted by ${I_j} = [{x_{j - 1}},{x_j}]$ related to ${x_j} = jh,(j = 1,2, \cdot{\cdot} \cdot ,k)$.

2.1.2 Adaptive partition

The interval $[0,d]$ is arbitrarily divided into k subintervals $[{x_1},{x_2}],[{x_2},{x_3}], \cdot{\cdot} \cdot ,[{x_k},{x_{k + 1}}]$. Assume $n_0^2(x)$ is a bounded variation function on the interval $[0,d]$, we have

$$V_{{x_i}}^{{x_{i + 1}}}({n_0^2(x)} )= \frac{1}{k}V_0^d({n_0^2(x)} ),$$
where $V_0^d({n_0^2(x)} )$ and $V_{{x_i}}^{{x_{i + 1}}}({n_0^2(x)} )$ are represented as the total variations of $n_0^2(x)$ on the intervals $[0,d]$ and $[{x_i},{x_{i + 1}}]$, respectively, $(i = 1,2, \cdot{\cdot} \cdot ,k).$

2.2 Interpolation methods

We have two forms as follows.

2.2.1 By the piecewise quadratic interpolation

The coefficients of Eq. (5) are given by

$$\left\{ \begin{array}{l} {a_j} = \frac{{2{{\kappa}_0}^2}}{{{h^2}}}[{n_0}^2({t_0}) - 2{n_0}^2({t_1}) + {n_0}^2({t_2})],\\ {b_j} = \frac{{{{\kappa}_0}^2}}{{{h^2}}}[(1 - 4j){n_0}^2({t_0}) + (8j - 4){n_0}^2({t_1}) + (3 - 4j){n_0}^2({t_2})],\\ {c_j} = {{\kappa}_0}^2[{(2{j^2} - j){n_0}^2({t_0}) + 4(j - {j^2}){n_0}^2({t_1}) + (2{j^2} - 3j + 1){n_0}^2({t_2})} ]- {{\beta}^2}. \end{array} \right.$$

Here ${t_0} = (j - 1)h$, ${t_1} = (j - 1/2)h$,and ${t_2} = jh,(j = 1,2, \cdot{\cdot} \cdot ,k)$.

2.2.2 By the quadratic spline interpolation

The coefficients of Eq. (5) are given by

$$\left\{ \begin{array}{l} {a_j} = \frac{{{M_{j + 1}} - {M_j}}}{{2({x_{j + 1}} - {x_j})}},\\ {b_j} = \frac{{{M_j}{x_{j + 1}} - {M_{j + 1}}{x_j}}}{{{x_{j + 1}} - {x_j}}},\\ {c_j} = \frac{{{M_{j + 1}}{x_j}^2 - {M_j}{x_{j + 1}}^2}}{{2({x_{j + 1}} - {x_j})}} + \frac{1}{2}{M_j}({x_{j + 1}} - {x_j}) + {\kappa}_0^2n_0^2({x_j}) - {{\beta}^2}, \end{array} \right.$$
where ${M_j} = S^{\prime}({x_j})$, $(j = 1,2, \cdots ,k)$, and $S(x)$ is the quadratic spline interpolation function.

Furthermore, ${M_j}$ and ${M_{j + 1}}$ satisfies

$${M_j} + {M_{j + 1}} = \frac{{2{{\kappa} _0}^2({n_0}^2({x_{j + 1}}) - {n_0}^2({x_j}))}}{{{x_{j + 1}} - {x_j}}}.$$

To fix ${M_j},(j = 1,2, \cdot{\cdot} \cdot ,k)$, we choose the below boundary condition: $S^{\prime\prime}(0) = {M_2} - {M_1} = 0$.

Moreover, Eq. (5) can be changed to the confluent hypergeometric equation, and a pair of linearly independent solutions $\{ {\alpha _j}(x),{{\beta} _j}(x)\}$ as $x \in [{x_j},{x_{j + 1}}]$ are given by the Kummer functions. Then the approximated solution of $\mathrm{\phi}$ as $x \in [{x_j},{x_{j + 1}}]$ can be expressed in the following form:

$${y_j}(x) = {A_j}{\alpha _j}(x) + B{}_j{{\beta} _j}(x)\;\;\;\;,\;\;\;\;\;\;(j = 1,2, \cdot \cdot \cdot ,k).$$

To obtain the constants of ${A_j}$ and ${B_j}$, $(j = 1,2, \cdots ,k)$, we require the solutions ${y_j}(x),(j = 1,2, \cdots ,k)$ of field $\mathrm{\phi} (x)$ and their first order derivatives are continuous. Thus, the interface conditions at ${x_j},(j = 1,2, \cdots ,k)$ are

$${y_j}({x_j}) = {y_{j + 1}}({x_j})\;\;\;\;,\;\;\;\;\;\;{y^{\prime}_j}({x_j}) = {y^{\prime}_{j + 1}}({x_j}).$$

And the outgoing boundary conditions at $x = 0$ and $x = d$ become

$${y^{\prime}_1}(0) ={-} i \sqrt {\mathrm{\kappa} _0^2n_1^2 - {\mathrm{\beta} ^2}} {y_1}(0), $$
$${y^{\prime}_k}(d) = i \sqrt {\mathrm{\kappa} _0^2n_2^2 - {\mathrm{\beta} ^2}} {y_k}(d). $$

According to these conditions, a linear system of ${A_j}$ and ${B_j}(j = 1,2, \cdot{\cdot} \cdot ,k)$ can be obtained.

$$\left\{ \begin{array}{l} {{\mathrm{\alpha}^{\prime}}_1}(0){A_1} + {{{\mathrm{\beta}}^{\prime}}_1}(0){B_1} ={-} i \sqrt {\mathrm{\kappa}_0^2n_1^2 - {{\mathrm{\beta}}^2}} [{{\mathrm{\alpha}_1}(0){A_1} + {{\mathrm{\beta}}_1}(0){B_1}} ],\\ {A_j}{\mathrm{\alpha}_j}({x_{j + 1}}) + {B_j}{{\mathrm{\beta}}_j}({x_{j + 1}}) = {A_{j + 1}}{\mathrm{\alpha}_{j + 1}}({x_{j + 1}}) + {B_{j + 1}}{{\mathrm{\beta}}_{j + 1}}({x_{j + 1}}),\\ {A_j}{{\mathrm{\alpha}^{\prime}}_j}({x_{j + 1}}) + {B_j}{{{\mathrm{\beta}}^{\prime}}_j}({x_{j + 1}}) = {A_{j + 1}}{{\mathrm{\alpha}^{\prime}}_{j + 1}}({x_{j + 1}}) + {B_{j + 1}}{{{\mathrm{\beta}}^{\prime}}_{j + 1}}({x_{j + 1}}),\\ {{\mathrm{\alpha}^{\prime}}_k}(d){A_k} + {{{\mathrm{\beta}}^{\prime}}_k}(d){B_k} = i \sqrt {\mathrm{\kappa}_0^2n_2^2 - {{\mathrm{\beta}}^2}} [{{\mathrm{\alpha}_k}(d){A_k} + {{\mathrm{\beta}}_k}(d){B_k}} ]. \end{array} \right.$$

In order to simplify these equations, let ${T_j} = {B_j}/{A_j},(j = 1,2, \cdot{\cdot} \cdot ,k)$, the dispersion relations of ${\beta}$ are

$$\left\{ \begin{array}{l} {T_1} ={-} \frac{{{{\alpha^{\prime}}_1}(0) + i {\alpha_1}(0){\gamma_1}}}{{{{{\beta}^{\prime}}_1}(0) + i {{\beta}_1}(0){\gamma_1}}},\\ {T_{j + 1}} ={-} \frac{{{\alpha_{j + 1}}(jh)[{{\alpha^{\prime}}_j}(jh) + {{{\beta}^{\prime}}_j}(jh){T_j}] + {{\alpha^{\prime}}_{j + 1}}(jh)[{\alpha_j}(jh) + {{\beta}_j}(jh){T_j}]}}{{{{\beta}_{j + 1}}(jh)[{{\alpha^{\prime}}_j}(jh) + {{{\beta}^{\prime}}_j}(jh){T_j}] + {{{\beta}^{\prime}}_{j + 1}}(jh)[{\alpha_j}(jh) + {{\beta}_j}(jh){T_j}]}},\\ \begin{array}{{ccc}} {}&{(j = 1,2, \cdots ,k - 1),}&{} \end{array}\\ {T_k} ={-} \frac{{{{\alpha^{\prime}}_k}(d) + i {\alpha_k}(d){\gamma_2}}}{{{{{\beta}^{\prime}}_k}(d) + i {{\beta}_k}(d){\gamma_2}}}. \end{array} \right.$$
where ${\gamma _1} = \sqrt {\mathrm{\kappa} _0^2n_1^2 - {\mathrm{\beta} ^2}}$ and ${\gamma _2} = \sqrt {\mathrm{\kappa} _0^2n_2^2 - {\mathrm{\beta} ^2}}$, and the evaluation on ${T_k}$ depends on ${T_1}$ recursively. Furthermore, ${T_k}$ is determined by ${\gamma _1}$. In addition, ${T_k}$ is given by the last equation, which depends on ${\gamma _2}$. Hence, the dispersion relation with respect to ${{\beta} ^2}$ is as follow:
$$F({\beta} ) = {T_k} + \frac{{{{\alpha ^{\prime}}_k}(d) + i {\alpha _k}(d){\gamma _2}}}{{{{{\beta} ^{\prime}}_k}(d) + i {{\beta} _k}(d){\gamma _2}}}.$$

2.3 Asymptotic solutions for the TE case

By the WKB method, we can obtain the following formulas [18].

$$\sqrt {{\kappa} _0^2n_0^2({0^ + }) - {{\beta} ^2}} \approx \frac{W}{{{a_0}}} - \frac{{a{a_2}}}{{{W^2}}} - \frac{{{a^2}{a_3}}}{{{W^3}}} - \frac{{{a^3}{a_4}}}{{{W^4}}}\textrm{ - } \cdot{\cdot} \cdot ,$$
where
$${a_0} ={-} \frac{{i d}}{2}\;\;;{a_2} = - \frac{{{\delta _0}\textrm{ + }{\delta _1}\textrm{ + }{\delta _2}}}{8} ;{a_3} = \frac{i }{{4d}}({{\delta_0}\textrm{ + }{\delta_1}\textrm{ + }{\delta_2}} );{a_4} = - \frac{{ - i d{b_4} - 2{a_3}}}{{2{a_0}}} - \frac{{{a_2}^2}}{2};$$
$${b_4} = - \frac{1}{{128}}({5{\delta_1}^2 + 5{\delta_2}^2 + {\delta_0}^2 + 2{\delta_1}{\delta_2} + 2{\delta_0}{\delta_2} + 2{\delta_0}{\delta_1}} );{\delta _0} = {\kappa} _0^2n_0^2({0^ + }) - {\kappa} _0^2n_0^2(d^-)$$
$${\delta _1} = {\kappa} _0^2n_0^2({0^ + }) - {\kappa} _0^2n_1^2\;\;\;;{\delta _2} = {\kappa} _0^2n_0^2({0^ + }) - {\kappa} _0^2n_2^2;{\delta _3} = {\delta _2} - {\delta _0},$$
$$W = {W_m} = \textrm{LambertW} (p,\frac{{ - {i ^{ - q}}}}{2} \cdot {({\delta _1}{\delta _3})^{\frac{1}{4}}}); m ={-} 4(p + 1) + q\;\;\;;(q = 0,1,2,3;p = - 1, - 2, \cdot \cdot \cdot )$$

For the sake of unification, let m be the rearrangement index (see [14]). In order to find the roots of the nonlinear equation $F({\beta} ) = 0$, we adopt Müller’s method [25] for implementation. Müller’s method requires three initial values, which are taken by the asymptotic solutions. The error analysis of the approximation is listed in the APPENDIX 1.

3. Numerical example and comparison

At last, the theory of this paper has been implemented and tested on a number of examples. Let ${n_0}(x) = 3.5 \cdot [1 - 0.08{(x - 1)^2}]$, ${n_1} = 2.1$, ${n_2} = 3.17$, $d = 2 \mu \textrm{m}$, $\lambda = 1.55 \mu \textrm{m}$, ${{\kappa} _0} = \frac{{2\pi }}{\lambda }$. The iteration stopping criteria is chosen as $|F |< {10^{ - 10}}$.

3.1 Comparison of the fitting effect by two kinds of spline interpolation methods

The fitting diagrams of the spline interpolation functions are shown in Figs. 2 and 3, respectively.

 

Fig. 2. The fitting function of the spline interpolation with $k = 4$ (marked in red), where the solid line (marked in blue) is represented as the original function (${\kappa} _0^2n_0^2(x)$).

Download Full Size | PPT Slide | PDF

 

Fig. 3. The fitting function of the spline interpolation with non-equidistant partition as $k = 4$ (marked in red), where the solid line (marked in blue) is represented as the original function (${\kappa} _0^2n_0^2(x)$).

Download Full Size | PPT Slide | PDF

It is demonstrated that the approximation effect of the non-equidistant and adaptive spline interpolation (see Fig. 3) is better than that of the equidistant spline interpolation (see Fig. 2).

We have found that two fitting curves obtained by the quadratic spline interpolation (equidistant partition) with $k = 16$ and the piecewise quadratic interpolation with $k = 2$ are almost consistent with the original curve $y = {\kappa} _0^2n_0^2(x)$, which matches the reality. The result is not a coincidence, for the original curve is a quadratic polynomial. Since the additional boundary condition at $x = 0$ of the quadratic spline interpolation is artificially assigned, rather than the boundary condition of the original curve, its spline interpolation function is not equal to the original one. However, the spline interpolation function can perfectly approximate the original function as long as the value of k is large enough.

3.2 Comparison of eigenvalues obtained by different spline interpolation methods

The diagrams of approximate eigenvalues as $k = 8$ and $k = 16$ are shown in Figs. 4 and 5, respectively. It is clear that the solutions ${\beta}$ of the spline interpolation are generally similar to the ones of piecewise interpolation except the accuracy as $k = 8$. To address this problem, we adopt non-equidistant partition to improve the accuracy, which can help achieve good results with $k = 8$ and $k = 16$ shown in Fig. 6.

 

Fig. 4. The propagation constants of the piecewise interpolation with equidistant partition, where ‘*’ stands for $\mathrm{\beta}$ obtained from the four-term WKB approximation, ‘o’ stands for results of $\mathrm{\beta}$ with $k = 8$, ‘+’ stands for results of ${\beta}$ with $k = 16$, and the horizontal axis $\textrm{Re}({\beta} /{{\kappa} _0})$ and the longitudinal axis $\textrm{Im(}{\beta} /{{\kappa} _0})$ represent the real and the imaginary parts of ${\beta} /{{\kappa} _0}$, respectively.

Download Full Size | PPT Slide | PDF

 

Fig. 5. The propagation constants of the spline interpolation with equidistant partition, where ‘*’ stands for $\mathrm{\beta}$ obtained from the four-term WKB approximation, ‘o’ stands for results of $\mathrm{\beta}$ with $k = 8$, ‘+’ stands for results of $\mathrm{\beta}$ with $k = 16$, and the horizontal axis $\textrm{Re}({\beta} /{{\kappa} _0})$ and the longitudinal axis $\textrm{Im(}{\beta} /{{\kappa} _0})$ represent the real and the imaginary parts of ${\beta} /{{\kappa} _0}$, respectively.

Download Full Size | PPT Slide | PDF

 

Fig. 6. The propagation constants of the spline interpolation with non-equidistant partition, where ‘*’ stands for $\mathrm{\beta}$ obtained from the four-term WKB approximation, ‘o’ stands for results of $\mathrm{\beta}$ with $k = 8$, ‘+’ stands for results of $\mathrm{\beta}$ with $k = 16$, and the horizontal axis $\textrm{Re}({\beta} /{{\kappa} _0})$ and the longitudinal axis $\textrm{Im(}{\beta} /{{\kappa} _0})$ represent the real and the imaginary parts of ${\beta} /{{\kappa} _0}$, respectively.

Download Full Size | PPT Slide | PDF

3.3 Numerical verification of the convergence rate

For the error analysis, we choose the absolute value of the relative error, which is denoted by ${E_r}$, and plot ${\log _2}({E_r})$ against ${\log _2}(h)$ in Figs. 78 and 9, respectively. In these figures (Figs. 79), as k increases, namely ${\log _2}(2/k)$ decreases, the slopes between the points are not less than 3. That is, our method meet the accuracy of order 3, and the accuracy can be improved by non-equidistant partition. That is to say, the order of convergence of general solutions (eigenmodes) reach the third.

 

Fig. 7. The relative error of the piecewise interpolation with equidistant partition, where ‘*’ stands for $m = 1$, ‘o’ stands for $m = 8$, ‘+’ stands for $m = 15$, the horizontal axis represents ${\log _2}(2/k)$, and the longitudinal axis represents ${\log _2}({E_r})$.

Download Full Size | PPT Slide | PDF

 

Fig. 8. The relative error of the spline interpolation with equidistant partition, where ‘*’ stands for $m = 1$, ‘o’ stands for $m = 8$, ‘+’ stands for $m = 15$, the horizontal axis represents ${\log _2}(2/k)$, and the longitudinal axis represents ${\log _2}({E_r})$.

Download Full Size | PPT Slide | PDF

 

Fig. 9. the relative error of spline interpolation with non-equidistant partition, where ‘*’ stands for $m = 1$, ‘o’ stands for $m = 8$, ‘+’ stands for $m = 15$, the horizontal axis represents ${\log _2}(2/k)$, and the longitudinal axis represents ${\log _2}({E_r})$.

Download Full Size | PPT Slide | PDF

3.4 Comparison of eigenfunctions obtained by two kinds of interpolation methods

The diagram of eigenfunctions are shown in Figs. 101112 and 13, respectively. The simulation results show that, the eigenfunction of the spline interpolation is smoother than the one of the piecewise interpolation.

 

Fig. 10. The real part of eigenfunction with $k = 4$. Where the solid line (marked in blue) stands for the spline interpolation, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Re}(\mathrm{\phi} (x))$ represents the real part of the eigenfunction function.

Download Full Size | PPT Slide | PDF

 

Fig. 11. The imaginary part of eigenfunction with $k = 4$, where the solid line (marked in blue) stands for the spline interpolation function, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Im(}\mathrm{\phi} (x))$ represents the imaginary part of the eigenfunction function.

Download Full Size | PPT Slide | PDF

 

Fig. 12. The real part of eigenfunction with $k = 16$, where the solid line (marked in blue) stands for the spline interpolation function, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Re}(\mathrm{\phi} (x))$ represents the real part of the eigenfunction function.

Download Full Size | PPT Slide | PDF

Obviously, the smoothness of the eigenfunction obtained by the spline interpolation can be well preserved regardless of the value of k. However, the piecewise interpolation cannot work like the spline interpolation, since if we apply the piecewise interpolation, there may be some jumping points. For example, in Figs. 1013, the jumping points appears near the beginning or the end of the eigenfunction. The reason is that piecewise interpolation only keeps the function continuous, but it cannot guarantee the continuity of the first derivative. Therefore, the deviation of numerical calculation appears, which in reality, is not consistent with the wave propagation.

 

Fig. 13. The imaginary part of eigenfunction with $k = 16$, where the solid line (marked in blue) stands for the spline interpolation function, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Im(}\mathrm{\phi} (x))$ represents the imaginary part of the eigenfunction function.

Download Full Size | PPT Slide | PDF

4. Conclusions

For a variety of refractive-index waveguide, the function ${\kappa} _0^2n_0^2(x)$ is approximated by a quadratic spline interpolation function. Then, it can establish to the approximate dispersion equation, only for TE case. Since the approximate equation is equivalent to the confluent hypergeometric equation, its’ solutions can be expressed analytically by the Kummer functions in each layer (sub-interval). In the numerical example, we find the roots of the dispersion equation (i.e. the propagation constants corresponding to the leaky modes) by Müller’s method, where three different asymptotic solutions play the roles of initial values. Obviously, the approximate dispersion equation converges fast to the exact one for the continuous refractive-index function, as the maximum value of the subinterval sizes tends to zero. In addition, to refine the accuracy, we applied the equidistant partition and the non-equidistant partition to divide the interval. We apply the piecewise quadratic interpolation and the quadratic spline interpolation to solve the approximate function, respectively. Numerical simulations demonstrate that the eigenfunction of the spline interpolation is smoother, compared with the one with the piecewise interpolation. Moreover, the non-equidistant partition can help promote the accuracy, and the order of convergence of general solutions meet the third.

Appendix 1

Error analysis: the error of the coefficient in Eq. (5) will inevitably lead to the error of the exact solution of the eigenvalue, and the error of eigenvalue can be obtained as follow:

$$({\hat{\lambda }\textrm{ - }\lambda } )\left( {\int_{{t_1}}^{{t_2}} {{\mathrm{\phi}^2}d x} + \sum\limits_{i = 1}^2 {\frac{{{\mathrm{\phi}^2}({t_i})}}{{2{\gamma_i}(\lambda )}}} } \right) \approx{-} \sum\limits_{i = 1}^k {\int_{{x_i}}^{{x_{i + 1}}} {R(f(x)) \cdot {\mathrm{\phi} ^2}d x} } ,$$
where ${{\beta} ^2} = \lambda$, ${\hat{{\beta} }^2} = \hat{\lambda }$, ${\gamma _i}(\lambda ) = \sqrt {\lambda - {{\kappa} _0}^2{n_i}^2}$, $R(f(x))$ is the remainder of interpolation $t{}_0 = 0 \,\, \textrm{and} \,\, t{}_1 = d$

Further,

$$|{\hat{\lambda }\textrm{ - }\lambda } |\le \frac{{2\sum\limits_{i = 1}^k {\int_{{x_i}}^{{x_{i + 1}}} {|{R(f(x))} |\cdot {{|\mathrm{\phi} |}^2}d x} } }}{{\left|{2\int_{{t_1}}^{{t_2}} {{\mathrm{\phi}^2}d x} + \frac{{{\mathrm{\phi}^2}(0)}}{{{\gamma_1}(\lambda )}} + \frac{{{\mathrm{\phi}^2}(d)}}{{{\gamma_2}(\lambda )}}} \right|}}.$$

Case 1: For the piecewise interpolation, we can obtain the remainder formula of the piecewise quadratic interpolation:

$$R(f(x))\textrm{ = }{R_1}(f(x)) = \frac{{f^{\prime\prime}{^{\prime}}(\xi )}}{6}{\omega _3}(x),$$
where $|{\omega _3}(x)|= |(x - {x_j})(x - \frac{{{x_j} + {x_{j + 1}}}}{2})(x - {x_{j + 1}})|$, which can reach the maximum value when $x = {{({x_j} + {x_{j + 1}})} / 2} + {3^{ - 1/2}}{{({x_{j + 1}} - {x_j})} / 2}$.

Let $h = \max \{ {h_1},{h_2}, \cdot{\cdot} \cdot ,{h_k}\}$, then the error estimation of the piecewise quadratic interpolation is shown as follow:

$$|{\hat{\lambda }\textrm{ - }\lambda } |\le \frac{{\sqrt 3 {h^3}}}{{108}} \cdot \frac{{\mathop {\sup }\limits_{x \in [0,d]} |{f^{\prime\prime}{^{\prime}}(x)} |\cdot \int_0^d {{{|\mathrm{\phi} |}^2}d x} }}{{\left|{2\int_0^d {{\mathrm{\phi}^2}\rm{d} x} + \frac{{{\mathrm{\phi}^2}(0)}}{{{\gamma_1}(\lambda )}} + \frac{{{\mathrm{\phi}^2}(d)}}{{{\gamma_2}(\lambda )}}} \right|}}.$$

Case 2: For the quadratic spline interpolation, we can obtain the remainder formula as follow [26]:

$$R(f(x)) = {R_2}(f(x)) = \int_0^d {K(x,t)f^{\prime\prime}{^{\prime}}(t)} d t,$$
where
$$K(x,t) = {R_2}\left\{ {\frac{{{{(x - t)}^2}_ + }}{{2!}}} \right\},\begin{array}{{cc}} {}&{\begin{array}{{cc}} {}&{} \end{array}} \end{array}$$
here ${R_2}\, (\cdot)$ is represented as the remainder of the quadratic spline interpolation.

The integral result of the function $K(x,t)$ on $x \in [{x_i},{x_{i + 1}}]$ becomes

$$\int_{{x_i}}^{{x_{i + 1}}} {|{K(x,t)} |d t} = \frac{1}{6}{(x - {x_i})^2}({x_{i + 1}} - x).$$

By the method of enlarging and reducing, we have

$$\int_0^d {|{K(x,t)} |d t} = \frac{1}{6}(x - {x_i})({x_{i + 1}} - x)(x - {x_i} + \sum\limits_{j = 1}^{i - 1} {{h_j}/{h_i}} ) < \frac{1}{{24}}{h_i}\sum\limits_{j = 1}^i {{h_j}^2} .$$

Thus, the form of remainder formula is

$$R(f(x)) = {R_2}(f(x)) = \frac{1}{{24}}\int_0^d {{h_i}\sum\limits_{j = 1}^i {{h_j}^2} |{f^{\prime\prime}{^{\prime}}(t)} |} d t.\begin{array}{{cc}} {}&{} \end{array}$$

Let $h = \max \{ {h_1},{h_2}, \cdot{\cdot} \cdot ,{h_k}\}$, we have

$$|{\hat{\lambda }\textrm{ - }\lambda } |< \frac{1}{{12}} \cdot \frac{{\mathop {\sup }\limits_{x \in [0,d]} |{f^{\prime\prime}{^{\prime}}(x)} |\cdot \sum\limits_{i = 1}^k {\int_{{x_i}}^{{x_{i + 1}}} {{h_i}\sum\limits_{j = 1}^i {{h_j}^2} } } |{{\mathrm{\phi}^2}} |d x}}{{\left|{2\int_0^d {{\mathrm{\phi}^2}d x} + \frac{{{\mathrm{\phi}^2}(0)}}{{{\gamma_1}(\lambda )}} + \frac{{{\mathrm{\phi}^2}(d)}}{{{\gamma_2}(\lambda )}}} \right|}}.$$

Particularly, if the third derivative of the original function is a bounded variation function, and the interval partition satisfies

$$\frac{1}{{{h_i}}}\sum\limits_{j = 1}^{i - 1} {|{{h_j} - {h_{j - 1}}} |} \le \frac{h}{\rm{d}},(i = 1,2, \cdot{\cdot} \cdot ,k),$$
a better residual estimate of the quadratic spline interpolation can be deduced.
$$|{{R_2}(f(x))} |\le \frac{1}{8}{h^3}\left( {\mathop {\sup }\limits_{x \in [0,d]} |{f^{\prime\prime}{^{\prime}}(x)} |+ V_0^d({f^{\prime\prime}{^{\prime}}(x)} )} \right).$$

Appendix 2

When $a \ne 0$, let’s begin to consider the second order differential equation as follow

$$y^{\prime\prime}(x) + (a{x^2} + bx + c)y(x) = 0,$$
by the change of variables $t = x + ({b / {2a}})$, we have
$$y^{\prime\prime}(t) + (a{t^2} + C)y(t) = 0,$$
where $C = c - [{{{b^2}} / {(4a)}}]$.

Further, the change of variables $z = \sqrt { - a} {t^2}$ and $\textrm{exp} (\textrm{ - }{z / 2})w(z) = y(t)$, leads to

$$z\frac{{{\textrm{\rm{d}}^2}w}}{{\textrm{\rm{d}}{z^2}}} + \left( {\frac{1}{2} - z} \right)\frac{{\textrm{\rm{d}}w}}{{\textrm{\rm{d}}z}} + (\frac{1}{4} - \frac{C}{{4\sqrt { - a} }})w(z) = 0,$$
which is in the form of the confluent hypergeometric equation
$$z\frac{{{\textrm{\rm{d}}^2}w}}{{\textrm{\rm{d}}{z^2}}} + ({B - z} )\frac{{\textrm{\rm{d}}w}}{{\textrm{\rm{d}}z}} + Aw(z) = 0,$$
with $A = \frac{1}{4}\textrm{ - }\frac{C}{{4\sqrt { - a} }}$ and $B = {1 / 2}$. Two linearly independent solutions are given by the Kummer functions $M(A,B,z)$ and $U(A,B,z)$. However, $U(A,B,z)$ shows a branch point $z = 0$. Hence, two linearly independent solutions of Eq. (30) are given by
$$\left\{ \begin{array}{l} \alpha (x) = \textrm{exp} (\textrm{ - }{{z(x)} / 2}) \cdot M\left( {A,\frac{1}{2},z(x)} \right),\\ {\beta} (x) = \textrm{exp} (\textrm{ - }{{z(x)} / 2}) \cdot {[z(x)]^{\frac{1}{2}}}M\left( {A + \frac{1}{2},\frac{3}{2},z(x)} \right) \end{array} \right.,{\kern 1pt}$$
where $z(x) = \sqrt { - a} {[{x + {b / {(2a)}}} ]^2}$, $A = \frac{C}{{4\sqrt { - a} }} = \frac{{{b^2} - 4ac}}{{16a\sqrt { - a} }}.$

Recall the following differential formulas for the Kummer function:

$$\left\{ \begin{array}{l} \alpha^{\prime}(x) ={-} \frac{1}{2}z^{\prime}(x)\alpha (x) + 2Az^{\prime}(x)\textrm{exp} (\textrm{ - }{{z(x)} / 2})M\left( {A + 1,\frac{3}{2},z(x)} \right),\\ {\beta}^{\prime}(x) ={-} \frac{{1 - z(x)}}{{2z(x)}}z^{\prime}(x){\beta} (x)\\ \begin{array}{{cc}} {\begin{array}{{cc}} {}&{} \end{array}}&{} \end{array} + \frac{{2A + 1}}{3}\sqrt {z(x)} z^{\prime}(x)\textrm{exp} [\textrm{ - }{{z(x)} / 2}]M\left( {A + \frac{3}{2},\frac{5}{2},z(x)} \right). \end{array} \right.$$
Special Case: when $a = 0$,and $b \ne 0$, Eq. (30) becomes
$$y^{\prime\prime}(x) + (bx + c)y(x) = 0.$$
and two linearly independent solutions are given by the Airy functions
$$\alpha (x) = {\rm Ai} (w(x)),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\beta} (x) = {\rm Bi} (w(x)),$$
where $w(x) ={-} \frac{{bx + c}}{{{b^{2/3}}}}$.

Therefore,

$$\alpha ^{\prime}(x) = \textrm{ - }{b^{1/3}}{\rm{A i}}^{\prime}(w(x)),{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\beta} ^{\prime}(x) = \textrm{ - }{b^{1/3}}{\rm{B i}}^{\prime}(w(x)).$$

The Airy functions and their derivatives can be evaluated in Matlab as

$$\textrm{airy} (z )= {\rm{Ai}} (z )\;\;\;,\;\;\;\textrm{airy}({1,z} )= {\rm{A i}}^{\prime}(z ),\\$$
$$\textrm{airy} ({2,z} )= {\rm{B i}} (z )\;\;,\;\;\;\textrm{airy} ({3,z} )= B i^{\prime}(z )$$

When $a = b = 0$, Eq. (30) becomes

$$y^{\prime\prime}(x) + cy(x) = 0,$$
which has two linearly independent solutions
$$\alpha (x) = {e^{ - i\sqrt c x}},{\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\kern 1pt} {\beta} (x) = {e^{i\sqrt c x}}.$$

Funding

Ministry of Science and Technology of the People's Republic of China (2018YFB1107600); Natural Science Foundation of Zhejiang Province (Q20A010007); First Class Discipline of Zhejiang - A (Zhejiang Gongshang University - Statistics).

Disclosures

The authors declare no conflicts of interest.

References

1. J. S. Bagby, D. P. Nyquist, and B. C. Drachman, “Integral formulation for analysis of integrated dielectric waveguides,” IEEE Trans. Microwave Theory Tech. MTT 33(10), 906–915 (1985). [CrossRef]  

2. R. März, Integrated Optics: Design and Modeling (Artech House, 1995).

3. J. A. DeSanto, Scalar Wave Theory (Springer-Verlag, 1992).

4. A. Ghatak, “Leaky modes in optical waveguides,” Opt. Quantum Electron. 17(5), 311–321 (1985). [CrossRef]  

5. J. Hu and C. R. Menyuk, “Understanding leaky modes: slab waveguide revisited,” Adv. Opt. Photon. 1(1), 58–106 (2009). [CrossRef]  

6. K. S. Chiang, “Review of numerical and approximate methods for the modal analysis of general optical dielectric waveguides,” Opt. Quantum Electron. 26(3), S113–S134 (1994). [CrossRef]  

7. A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014). [CrossRef]  

8. A. Snyder and J. Love, Optical Waveguide Theory (Springer, 1983).

9. C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron 11(2), 457–465 (2005). [CrossRef]  

10. Z. E. Abid, K. L. Johnson, and A. Gopinath, “Analysis of dielectric guides by vector transverse magnetic field finite elements,” J. Lightwave Technol. 11(10), 1545–1549 (1993). [CrossRef]  

11. P. J. Ciang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008). [CrossRef]  

12. S. F. Chiang, B. Y. Lin, C. H. Teng, C. Y. Wang, and S. Y. Chung, “A multidomain pseudospectral mode solver for optical waveguide analysis,” J. Lightwave Technol. 30(13), 2077–2087 (2012). [CrossRef]  

13. M. Walz, T. Zebrowski, J. Kuchenmeister, and K. Büsch, “B-spline modal method: A polynomial approach compared to the Fourier modal method,” Opt. Express 21(12), 14683–14697 (2013). [CrossRef]  

14. J. Zhu and Y. Lu, “Leaky modes of slab waveguides-asymptotic solutions,” J. Lightwave Technol. 24(3), 1619–1623 (2006). [CrossRef]  

15. L. Knockaert, H. Rogier, and D. De Zutter, “An FFT-based signal identification approach for obtaining the propagation constants of the leaky modes in layered media,” AEU Int. J. Electron. Commun. 59(4), 230–238 (2005). [CrossRef]  

16. C. C. Huang, “Numerical calculations of ARROW structures by pseudospectral approach with Mur’s absorbing boundary conditions,” Opt. Express 14(24), 11631–11652 (2006). [CrossRef]  

17. D. Song and Y. Lu, “Pesudospectral modal method for computing optical waveguide modes,” J. Lightwave Technol. 32(8), 1624–1630 (2014). [CrossRef]  

18. J. Zhu, Z. Chen, and S. Tang, “Leaky modes of optical waveguides with varied refractive index for microchip optical interconnect applications – Asymptotic solutions,” Microelecton. Reliab. 48(4), 555–562 (2008). [CrossRef]  

19. S. Khorasani and K. Mehrany, “Differential transfer-matrix method for solution of one-dimensional linearnonhomogeneous optical structures,” J. Opt. Soc. Am. B 20(1), 91–96 (2003). [CrossRef]  

20. M. Eghlidi, K. Mehrany, and B. Rashidian, “Modified differential transfer-matrix method for solution ofonedimensional linear inhomogeneous optical structures,” J. Opt. Soc. Am. B 22(7), 1521–1528 (2005). [CrossRef]  

21. J. Zhu and Z. Shen, “Dispersion relation of leaky modes in nonhomogeneous waveguides and its applications,” J. Lightwave Technol. 29(21), 3230–3236 (2011). [CrossRef]  

22. J. Zhu and J. Zheng, “Nonlinear equation of the modes in circular slab waveguides and its application,” Appl. Opt. 52(33), 8013–8023 (2013). [CrossRef]  

23. Y. Li and J. Zhu, “Efficient approximations of dispersion relations in optical waveguides with varying refractive-index profiles,” Opt. Express 23(9), 11952–11964 (2015). [CrossRef]  

24. S. Pruess, “Estimating the eigenvalues of Sturm-Liouville problems by approximating the differential equation,” SIAM J. Numer. Anal. 10(1), 55–68 (1973). [CrossRef]  

25. D. E. Muller, “A method for solving algebraic equations using an automatic computer,” Mathematical Tables and Other Aids to Computation 10(56), 208–215 (1956). [CrossRef]  

26. J. Peetre, Lecture Notes: A Theory of Interpolation of Normed Spaces (Brazilia, 1963).

References

  • View by:
  • |
  • |
  • |

  1. J. S. Bagby, D. P. Nyquist, and B. C. Drachman, “Integral formulation for analysis of integrated dielectric waveguides,” IEEE Trans. Microwave Theory Tech. MTT 33(10), 906–915 (1985).
    [Crossref]
  2. R. März, Integrated Optics: Design and Modeling (Artech House, 1995).
  3. J. A. DeSanto, Scalar Wave Theory (Springer-Verlag, 1992).
  4. A. Ghatak, “Leaky modes in optical waveguides,” Opt. Quantum Electron. 17(5), 311–321 (1985).
    [Crossref]
  5. J. Hu and C. R. Menyuk, “Understanding leaky modes: slab waveguide revisited,” Adv. Opt. Photon. 1(1), 58–106 (2009).
    [Crossref]
  6. K. S. Chiang, “Review of numerical and approximate methods for the modal analysis of general optical dielectric waveguides,” Opt. Quantum Electron. 26(3), S113–S134 (1994).
    [Crossref]
  7. A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014).
    [Crossref]
  8. A. Snyder and J. Love, Optical Waveguide Theory (Springer, 1983).
  9. C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron 11(2), 457–465 (2005).
    [Crossref]
  10. Z. E. Abid, K. L. Johnson, and A. Gopinath, “Analysis of dielectric guides by vector transverse magnetic field finite elements,” J. Lightwave Technol. 11(10), 1545–1549 (1993).
    [Crossref]
  11. P. J. Ciang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008).
    [Crossref]
  12. S. F. Chiang, B. Y. Lin, C. H. Teng, C. Y. Wang, and S. Y. Chung, “A multidomain pseudospectral mode solver for optical waveguide analysis,” J. Lightwave Technol. 30(13), 2077–2087 (2012).
    [Crossref]
  13. M. Walz, T. Zebrowski, J. Kuchenmeister, and K. Büsch, “B-spline modal method: A polynomial approach compared to the Fourier modal method,” Opt. Express 21(12), 14683–14697 (2013).
    [Crossref]
  14. J. Zhu and Y. Lu, “Leaky modes of slab waveguides-asymptotic solutions,” J. Lightwave Technol. 24(3), 1619–1623 (2006).
    [Crossref]
  15. L. Knockaert, H. Rogier, and D. De Zutter, “An FFT-based signal identification approach for obtaining the propagation constants of the leaky modes in layered media,” AEU Int. J. Electron. Commun. 59(4), 230–238 (2005).
    [Crossref]
  16. C. C. Huang, “Numerical calculations of ARROW structures by pseudospectral approach with Mur’s absorbing boundary conditions,” Opt. Express 14(24), 11631–11652 (2006).
    [Crossref]
  17. D. Song and Y. Lu, “Pesudospectral modal method for computing optical waveguide modes,” J. Lightwave Technol. 32(8), 1624–1630 (2014).
    [Crossref]
  18. J. Zhu, Z. Chen, and S. Tang, “Leaky modes of optical waveguides with varied refractive index for microchip optical interconnect applications – Asymptotic solutions,” Microelecton. Reliab. 48(4), 555–562 (2008).
    [Crossref]
  19. S. Khorasani and K. Mehrany, “Differential transfer-matrix method for solution of one-dimensional linearnonhomogeneous optical structures,” J. Opt. Soc. Am. B 20(1), 91–96 (2003).
    [Crossref]
  20. M. Eghlidi, K. Mehrany, and B. Rashidian, “Modified differential transfer-matrix method for solution ofonedimensional linear inhomogeneous optical structures,” J. Opt. Soc. Am. B 22(7), 1521–1528 (2005).
    [Crossref]
  21. J. Zhu and Z. Shen, “Dispersion relation of leaky modes in nonhomogeneous waveguides and its applications,” J. Lightwave Technol. 29(21), 3230–3236 (2011).
    [Crossref]
  22. J. Zhu and J. Zheng, “Nonlinear equation of the modes in circular slab waveguides and its application,” Appl. Opt. 52(33), 8013–8023 (2013).
    [Crossref]
  23. Y. Li and J. Zhu, “Efficient approximations of dispersion relations in optical waveguides with varying refractive-index profiles,” Opt. Express 23(9), 11952–11964 (2015).
    [Crossref]
  24. S. Pruess, “Estimating the eigenvalues of Sturm-Liouville problems by approximating the differential equation,” SIAM J. Numer. Anal. 10(1), 55–68 (1973).
    [Crossref]
  25. D. E. Muller, “A method for solving algebraic equations using an automatic computer,” Mathematical Tables and Other Aids to Computation 10(56), 208–215 (1956).
    [Crossref]
  26. J. Peetre, Lecture Notes: A Theory of Interpolation of Normed Spaces (Brazilia, 1963).

2015 (1)

2014 (2)

A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014).
[Crossref]

D. Song and Y. Lu, “Pesudospectral modal method for computing optical waveguide modes,” J. Lightwave Technol. 32(8), 1624–1630 (2014).
[Crossref]

2013 (2)

2012 (1)

2011 (1)

2009 (1)

2008 (2)

P. J. Ciang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008).
[Crossref]

J. Zhu, Z. Chen, and S. Tang, “Leaky modes of optical waveguides with varied refractive index for microchip optical interconnect applications – Asymptotic solutions,” Microelecton. Reliab. 48(4), 555–562 (2008).
[Crossref]

2006 (2)

2005 (3)

M. Eghlidi, K. Mehrany, and B. Rashidian, “Modified differential transfer-matrix method for solution ofonedimensional linear inhomogeneous optical structures,” J. Opt. Soc. Am. B 22(7), 1521–1528 (2005).
[Crossref]

L. Knockaert, H. Rogier, and D. De Zutter, “An FFT-based signal identification approach for obtaining the propagation constants of the leaky modes in layered media,” AEU Int. J. Electron. Commun. 59(4), 230–238 (2005).
[Crossref]

C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron 11(2), 457–465 (2005).
[Crossref]

2003 (1)

1994 (1)

K. S. Chiang, “Review of numerical and approximate methods for the modal analysis of general optical dielectric waveguides,” Opt. Quantum Electron. 26(3), S113–S134 (1994).
[Crossref]

1993 (1)

Z. E. Abid, K. L. Johnson, and A. Gopinath, “Analysis of dielectric guides by vector transverse magnetic field finite elements,” J. Lightwave Technol. 11(10), 1545–1549 (1993).
[Crossref]

1985 (2)

J. S. Bagby, D. P. Nyquist, and B. C. Drachman, “Integral formulation for analysis of integrated dielectric waveguides,” IEEE Trans. Microwave Theory Tech. MTT 33(10), 906–915 (1985).
[Crossref]

A. Ghatak, “Leaky modes in optical waveguides,” Opt. Quantum Electron. 17(5), 311–321 (1985).
[Crossref]

1973 (1)

S. Pruess, “Estimating the eigenvalues of Sturm-Liouville problems by approximating the differential equation,” SIAM J. Numer. Anal. 10(1), 55–68 (1973).
[Crossref]

1956 (1)

D. E. Muller, “A method for solving algebraic equations using an automatic computer,” Mathematical Tables and Other Aids to Computation 10(56), 208–215 (1956).
[Crossref]

Abid, Z. E.

Z. E. Abid, K. L. Johnson, and A. Gopinath, “Analysis of dielectric guides by vector transverse magnetic field finite elements,” J. Lightwave Technol. 11(10), 1545–1549 (1993).
[Crossref]

Alù, A.

A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014).
[Crossref]

Bagby, J. S.

J. S. Bagby, D. P. Nyquist, and B. C. Drachman, “Integral formulation for analysis of integrated dielectric waveguides,” IEEE Trans. Microwave Theory Tech. MTT 33(10), 906–915 (1985).
[Crossref]

Büsch, K.

Castaldi, G.

A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014).
[Crossref]

Chang, H. C.

P. J. Ciang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008).
[Crossref]

Chen, Z.

J. Zhu, Z. Chen, and S. Tang, “Leaky modes of optical waveguides with varied refractive index for microchip optical interconnect applications – Asymptotic solutions,” Microelecton. Reliab. 48(4), 555–562 (2008).
[Crossref]

Chiang, K. S.

K. S. Chiang, “Review of numerical and approximate methods for the modal analysis of general optical dielectric waveguides,” Opt. Quantum Electron. 26(3), S113–S134 (1994).
[Crossref]

Chiang, S. F.

Chung, S. Y.

Ciang, P. J.

P. J. Ciang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008).
[Crossref]

De Zutter, D.

L. Knockaert, H. Rogier, and D. De Zutter, “An FFT-based signal identification approach for obtaining the propagation constants of the leaky modes in layered media,” AEU Int. J. Electron. Commun. 59(4), 230–238 (2005).
[Crossref]

DeSanto, J. A.

J. A. DeSanto, Scalar Wave Theory (Springer-Verlag, 1992).

Drachman, B. C.

J. S. Bagby, D. P. Nyquist, and B. C. Drachman, “Integral formulation for analysis of integrated dielectric waveguides,” IEEE Trans. Microwave Theory Tech. MTT 33(10), 906–915 (1985).
[Crossref]

Eghlidi, M.

Engheta, N.

A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014).
[Crossref]

Galdi, V.

A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014).
[Crossref]

Ghatak, A.

A. Ghatak, “Leaky modes in optical waveguides,” Opt. Quantum Electron. 17(5), 311–321 (1985).
[Crossref]

Gopinath, A.

Z. E. Abid, K. L. Johnson, and A. Gopinath, “Analysis of dielectric guides by vector transverse magnetic field finite elements,” J. Lightwave Technol. 11(10), 1545–1549 (1993).
[Crossref]

Hu, J.

Huang, C. C.

C. C. Huang, “Numerical calculations of ARROW structures by pseudospectral approach with Mur’s absorbing boundary conditions,” Opt. Express 14(24), 11631–11652 (2006).
[Crossref]

C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron 11(2), 457–465 (2005).
[Crossref]

C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron 11(2), 457–465 (2005).
[Crossref]

Johnson, K. L.

Z. E. Abid, K. L. Johnson, and A. Gopinath, “Analysis of dielectric guides by vector transverse magnetic field finite elements,” J. Lightwave Technol. 11(10), 1545–1549 (1993).
[Crossref]

Khorasani, S.

Knockaert, L.

L. Knockaert, H. Rogier, and D. De Zutter, “An FFT-based signal identification approach for obtaining the propagation constants of the leaky modes in layered media,” AEU Int. J. Electron. Commun. 59(4), 230–238 (2005).
[Crossref]

Kuchenmeister, J.

Li, Y.

Lin, B. Y.

Love, J.

A. Snyder and J. Love, Optical Waveguide Theory (Springer, 1983).

Lu, Y.

März, R.

R. März, Integrated Optics: Design and Modeling (Artech House, 1995).

Mehrany, K.

Menyuk, C. R.

Monticone, F.

A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014).
[Crossref]

Muller, D. E.

D. E. Muller, “A method for solving algebraic equations using an automatic computer,” Mathematical Tables and Other Aids to Computation 10(56), 208–215 (1956).
[Crossref]

Nyquist, D. P.

J. S. Bagby, D. P. Nyquist, and B. C. Drachman, “Integral formulation for analysis of integrated dielectric waveguides,” IEEE Trans. Microwave Theory Tech. MTT 33(10), 906–915 (1985).
[Crossref]

Peetre, J.

J. Peetre, Lecture Notes: A Theory of Interpolation of Normed Spaces (Brazilia, 1963).

Pruess, S.

S. Pruess, “Estimating the eigenvalues of Sturm-Liouville problems by approximating the differential equation,” SIAM J. Numer. Anal. 10(1), 55–68 (1973).
[Crossref]

Rashidian, B.

Rogier, H.

L. Knockaert, H. Rogier, and D. De Zutter, “An FFT-based signal identification approach for obtaining the propagation constants of the leaky modes in layered media,” AEU Int. J. Electron. Commun. 59(4), 230–238 (2005).
[Crossref]

Shen, Z.

Silva, A.

A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014).
[Crossref]

Snyder, A.

A. Snyder and J. Love, Optical Waveguide Theory (Springer, 1983).

Song, D.

Tang, S.

J. Zhu, Z. Chen, and S. Tang, “Leaky modes of optical waveguides with varied refractive index for microchip optical interconnect applications – Asymptotic solutions,” Microelecton. Reliab. 48(4), 555–562 (2008).
[Crossref]

Teng, C. H.

S. F. Chiang, B. Y. Lin, C. H. Teng, C. Y. Wang, and S. Y. Chung, “A multidomain pseudospectral mode solver for optical waveguide analysis,” J. Lightwave Technol. 30(13), 2077–2087 (2012).
[Crossref]

P. J. Ciang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008).
[Crossref]

Walz, M.

Wang, C. Y.

Wu, C. L.

P. J. Ciang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008).
[Crossref]

Yang, C. S.

P. J. Ciang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008).
[Crossref]

Yang, J. Y.

C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron 11(2), 457–465 (2005).
[Crossref]

Zebrowski, T.

Zheng, J.

Zhu, J.

Adv. Opt. Photon. (1)

AEU Int. J. Electron. Commun. (1)

L. Knockaert, H. Rogier, and D. De Zutter, “An FFT-based signal identification approach for obtaining the propagation constants of the leaky modes in layered media,” AEU Int. J. Electron. Commun. 59(4), 230–238 (2005).
[Crossref]

Appl. Opt. (1)

IEEE J. Quantum Electron. (1)

P. J. Ciang, C. L. Wu, C. H. Teng, C. S. Yang, and H. C. Chang, “Full-vectorial optical waveguide mode solvers using multidomain pseudospectral frequency-domain (PSFD) formulations,” IEEE J. Quantum Electron. 44(1), 56–66 (2008).
[Crossref]

IEEE J. Sel. Top. Quantum Electron (1)

C. C. Huang, C. C. Huang, and J. Y. Yang, “A full-vectorial pseudospectral modal analysis of dielectric optical waveguides with stepped refractive index profiles,” IEEE J. Sel. Top. Quantum Electron 11(2), 457–465 (2005).
[Crossref]

IEEE Trans. Microwave Theory Tech. MTT (1)

J. S. Bagby, D. P. Nyquist, and B. C. Drachman, “Integral formulation for analysis of integrated dielectric waveguides,” IEEE Trans. Microwave Theory Tech. MTT 33(10), 906–915 (1985).
[Crossref]

J. Lightwave Technol. (5)

J. Opt. Soc. Am. B (2)

Mathematical Tables and Other Aids to Computation (1)

D. E. Muller, “A method for solving algebraic equations using an automatic computer,” Mathematical Tables and Other Aids to Computation 10(56), 208–215 (1956).
[Crossref]

Microelecton. Reliab. (1)

J. Zhu, Z. Chen, and S. Tang, “Leaky modes of optical waveguides with varied refractive index for microchip optical interconnect applications – Asymptotic solutions,” Microelecton. Reliab. 48(4), 555–562 (2008).
[Crossref]

Opt. Express (3)

Opt. Quantum Electron. (2)

A. Ghatak, “Leaky modes in optical waveguides,” Opt. Quantum Electron. 17(5), 311–321 (1985).
[Crossref]

K. S. Chiang, “Review of numerical and approximate methods for the modal analysis of general optical dielectric waveguides,” Opt. Quantum Electron. 26(3), S113–S134 (1994).
[Crossref]

Science (1)

A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, “Performing mathematical operations with metamaterials,” Science 343(6167), 160–163 (2014).
[Crossref]

SIAM J. Numer. Anal. (1)

S. Pruess, “Estimating the eigenvalues of Sturm-Liouville problems by approximating the differential equation,” SIAM J. Numer. Anal. 10(1), 55–68 (1973).
[Crossref]

Other (4)

J. Peetre, Lecture Notes: A Theory of Interpolation of Normed Spaces (Brazilia, 1963).

A. Snyder and J. Love, Optical Waveguide Theory (Springer, 1983).

R. März, Integrated Optics: Design and Modeling (Artech House, 1995).

J. A. DeSanto, Scalar Wave Theory (Springer-Verlag, 1992).

Cited By

OSA participates in Crossref's Cited-By Linking service. Citing articles from OSA journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (13)

Fig. 1.
Fig. 1. Subdivision diagram of a three-layers’ waveguide (Ref. [23], Fig. 1), where $d = k \cdot h, k = 3$, and $h$ is a step width.
Fig. 2.
Fig. 2. The fitting function of the spline interpolation with $k = 4$ (marked in red), where the solid line (marked in blue) is represented as the original function (${\kappa} _0^2n_0^2(x)$).
Fig. 3.
Fig. 3. The fitting function of the spline interpolation with non-equidistant partition as $k = 4$ (marked in red), where the solid line (marked in blue) is represented as the original function (${\kappa} _0^2n_0^2(x)$).
Fig. 4.
Fig. 4. The propagation constants of the piecewise interpolation with equidistant partition, where ‘*’ stands for $\mathrm{\beta}$ obtained from the four-term WKB approximation, ‘o’ stands for results of $\mathrm{\beta}$ with $k = 8$, ‘+’ stands for results of ${\beta}$ with $k = 16$, and the horizontal axis $\textrm{Re}({\beta} /{{\kappa} _0})$ and the longitudinal axis $\textrm{Im(}{\beta} /{{\kappa} _0})$ represent the real and the imaginary parts of ${\beta} /{{\kappa} _0}$, respectively.
Fig. 5.
Fig. 5. The propagation constants of the spline interpolation with equidistant partition, where ‘*’ stands for $\mathrm{\beta}$ obtained from the four-term WKB approximation, ‘o’ stands for results of $\mathrm{\beta}$ with $k = 8$, ‘+’ stands for results of $\mathrm{\beta}$ with $k = 16$, and the horizontal axis $\textrm{Re}({\beta} /{{\kappa} _0})$ and the longitudinal axis $\textrm{Im(}{\beta} /{{\kappa} _0})$ represent the real and the imaginary parts of ${\beta} /{{\kappa} _0}$, respectively.
Fig. 6.
Fig. 6. The propagation constants of the spline interpolation with non-equidistant partition, where ‘*’ stands for $\mathrm{\beta}$ obtained from the four-term WKB approximation, ‘o’ stands for results of $\mathrm{\beta}$ with $k = 8$, ‘+’ stands for results of $\mathrm{\beta}$ with $k = 16$, and the horizontal axis $\textrm{Re}({\beta} /{{\kappa} _0})$ and the longitudinal axis $\textrm{Im(}{\beta} /{{\kappa} _0})$ represent the real and the imaginary parts of ${\beta} /{{\kappa} _0}$, respectively.
Fig. 7.
Fig. 7. The relative error of the piecewise interpolation with equidistant partition, where ‘*’ stands for $m = 1$, ‘o’ stands for $m = 8$, ‘+’ stands for $m = 15$, the horizontal axis represents ${\log _2}(2/k)$, and the longitudinal axis represents ${\log _2}({E_r})$.
Fig. 8.
Fig. 8. The relative error of the spline interpolation with equidistant partition, where ‘*’ stands for $m = 1$, ‘o’ stands for $m = 8$, ‘+’ stands for $m = 15$, the horizontal axis represents ${\log _2}(2/k)$, and the longitudinal axis represents ${\log _2}({E_r})$.
Fig. 9.
Fig. 9. the relative error of spline interpolation with non-equidistant partition, where ‘*’ stands for $m = 1$, ‘o’ stands for $m = 8$, ‘+’ stands for $m = 15$, the horizontal axis represents ${\log _2}(2/k)$, and the longitudinal axis represents ${\log _2}({E_r})$.
Fig. 10.
Fig. 10. The real part of eigenfunction with $k = 4$. Where the solid line (marked in blue) stands for the spline interpolation, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Re}(\mathrm{\phi} (x))$ represents the real part of the eigenfunction function.
Fig. 11.
Fig. 11. The imaginary part of eigenfunction with $k = 4$, where the solid line (marked in blue) stands for the spline interpolation function, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Im(}\mathrm{\phi} (x))$ represents the imaginary part of the eigenfunction function.
Fig. 12.
Fig. 12. The real part of eigenfunction with $k = 16$, where the solid line (marked in blue) stands for the spline interpolation function, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Re}(\mathrm{\phi} (x))$ represents the real part of the eigenfunction function.
Fig. 13.
Fig. 13. The imaginary part of eigenfunction with $k = 16$, where the solid line (marked in blue) stands for the spline interpolation function, the dotted line (marked in red) stands for the piecewise interpolation function, the horizontal axis x represents the position of the point, and the longitudinal axis $\textrm{Im(}\mathrm{\phi} (x))$ represents the imaginary part of the eigenfunction function.

Equations (47)

Equations on this page are rendered with MathJax. Learn more.

n = { n 1 , x < 0 , n 0 ( x ) , 0 x d , n 2 , x > d ,
{ ρ d d x ( 1 ρ d ϕ d x ) + κ 0 2 n 0 2 ( x ) ϕ = β 2 ϕ , < x < , lim | x | + ϕ = 0.
d 2 ϕ d x 2 + κ 0 2 n 0 2 ( x ) ϕ = β 2 ϕ , 0 < x < d ,
d ϕ d x = i κ 0 2 n 1 2 β 2 ϕ , x = 0 + ,
d ϕ d x = i κ 0 2 n 2 2 β 2 ϕ , x = d ,
\rm{d} 2 y j \rm{d} x 2 + ( a j x 2 + b j x + c j ) y j = 0 ,
V x i x i + 1 ( n 0 2 ( x ) ) = 1 k V 0 d ( n 0 2 ( x ) ) ,
{ a j = 2 κ 0 2 h 2 [ n 0 2 ( t 0 ) 2 n 0 2 ( t 1 ) + n 0 2 ( t 2 ) ] , b j = κ 0 2 h 2 [ ( 1 4 j ) n 0 2 ( t 0 ) + ( 8 j 4 ) n 0 2 ( t 1 ) + ( 3 4 j ) n 0 2 ( t 2 ) ] , c j = κ 0 2 [ ( 2 j 2 j ) n 0 2 ( t 0 ) + 4 ( j j 2 ) n 0 2 ( t 1 ) + ( 2 j 2 3 j + 1 ) n 0 2 ( t 2 ) ] β 2 .
{ a j = M j + 1 M j 2 ( x j + 1 x j ) , b j = M j x j + 1 M j + 1 x j x j + 1 x j , c j = M j + 1 x j 2 M j x j + 1 2 2 ( x j + 1 x j ) + 1 2 M j ( x j + 1 x j ) + κ 0 2 n 0 2 ( x j ) β 2 ,
M j + M j + 1 = 2 κ 0 2 ( n 0 2 ( x j + 1 ) n 0 2 ( x j ) ) x j + 1 x j .
y j ( x ) = A j α j ( x ) + B j β j ( x ) , ( j = 1 , 2 , , k ) .
y j ( x j ) = y j + 1 ( x j ) , y j ( x j ) = y j + 1 ( x j ) .
y 1 ( 0 ) = i κ 0 2 n 1 2 β 2 y 1 ( 0 ) ,
y k ( d ) = i κ 0 2 n 2 2 β 2 y k ( d ) .
{ α 1 ( 0 ) A 1 + β 1 ( 0 ) B 1 = i κ 0 2 n 1 2 β 2 [ α 1 ( 0 ) A 1 + β 1 ( 0 ) B 1 ] , A j α j ( x j + 1 ) + B j β j ( x j + 1 ) = A j + 1 α j + 1 ( x j + 1 ) + B j + 1 β j + 1 ( x j + 1 ) , A j α j ( x j + 1 ) + B j β j ( x j + 1 ) = A j + 1 α j + 1 ( x j + 1 ) + B j + 1 β j + 1 ( x j + 1 ) , α k ( d ) A k + β k ( d ) B k = i κ 0 2 n 2 2 β 2 [ α k ( d ) A k + β k ( d ) B k ] .
{ T 1 = α 1 ( 0 ) + i α 1 ( 0 ) γ 1 β 1 ( 0 ) + i β 1 ( 0 ) γ 1 , T j + 1 = α j + 1 ( j h ) [ α j ( j h ) + β j ( j h ) T j ] + α j + 1 ( j h ) [ α j ( j h ) + β j ( j h ) T j ] β j + 1 ( j h ) [ α j ( j h ) + β j ( j h ) T j ] + β j + 1 ( j h ) [ α j ( j h ) + β j ( j h ) T j ] , ( j = 1 , 2 , , k 1 ) , T k = α k ( d ) + i α k ( d ) γ 2 β k ( d ) + i β k ( d ) γ 2 .
F ( β ) = T k + α k ( d ) + i α k ( d ) γ 2 β k ( d ) + i β k ( d ) γ 2 .
κ 0 2 n 0 2 ( 0 + ) β 2 W a 0 a a 2 W 2 a 2 a 3 W 3 a 3 a 4 W 4  -  ,
a 0 = i d 2 ; a 2 = δ 0  +  δ 1  +  δ 2 8 ; a 3 = i 4 d ( δ 0  +  δ 1  +  δ 2 ) ; a 4 = i d b 4 2 a 3 2 a 0 a 2 2 2 ;
b 4 = 1 128 ( 5 δ 1 2 + 5 δ 2 2 + δ 0 2 + 2 δ 1 δ 2 + 2 δ 0 δ 2 + 2 δ 0 δ 1 ) ; δ 0 = κ 0 2 n 0 2 ( 0 + ) κ 0 2 n 0 2 ( d )
δ 1 = κ 0 2 n 0 2 ( 0 + ) κ 0 2 n 1 2 ; δ 2 = κ 0 2 n 0 2 ( 0 + ) κ 0 2 n 2 2 ; δ 3 = δ 2 δ 0 ,
W = W m = LambertW ( p , i q 2 ( δ 1 δ 3 ) 1 4 ) ; m = 4 ( p + 1 ) + q ; ( q = 0 , 1 , 2 , 3 ; p = 1 , 2 , )
( λ ^  -  λ ) ( t 1 t 2 ϕ 2 d x + i = 1 2 ϕ 2 ( t i ) 2 γ i ( λ ) ) i = 1 k x i x i + 1 R ( f ( x ) ) ϕ 2 d x ,
| λ ^  -  λ | 2 i = 1 k x i x i + 1 | R ( f ( x ) ) | | ϕ | 2 d x | 2 t 1 t 2 ϕ 2 d x + ϕ 2 ( 0 ) γ 1 ( λ ) + ϕ 2 ( d ) γ 2 ( λ ) | .
R ( f ( x ) )  =  R 1 ( f ( x ) ) = f ( ξ ) 6 ω 3 ( x ) ,
| λ ^  -  λ | 3 h 3 108 sup x [ 0 , d ] | f ( x ) | 0 d | ϕ | 2 d x | 2 0 d ϕ 2 d x + ϕ 2 ( 0 ) γ 1 ( λ ) + ϕ 2 ( d ) γ 2 ( λ ) | .
R ( f ( x ) ) = R 2 ( f ( x ) ) = 0 d K ( x , t ) f ( t ) d t ,
K ( x , t ) = R 2 { ( x t ) 2 + 2 ! } ,
x i x i + 1 | K ( x , t ) | d t = 1 6 ( x x i ) 2 ( x i + 1 x ) .
0 d | K ( x , t ) | d t = 1 6 ( x x i ) ( x i + 1 x ) ( x x i + j = 1 i 1 h j / h i ) < 1 24 h i j = 1 i h j 2 .
R ( f ( x ) ) = R 2 ( f ( x ) ) = 1 24 0 d h i j = 1 i h j 2 | f ( t ) | d t .
| λ ^  -  λ | < 1 12 sup x [ 0 , d ] | f ( x ) | i = 1 k x i x i + 1 h i j = 1 i h j 2 | ϕ 2 | d x | 2 0 d ϕ 2 d x + ϕ 2 ( 0 ) γ 1 ( λ ) + ϕ 2 ( d ) γ 2 ( λ ) | .
1 h i j = 1 i 1 | h j h j 1 | h d , ( i = 1 , 2 , , k ) ,
| R 2 ( f ( x ) ) | 1 8 h 3 ( sup x [ 0 , d ] | f ( x ) | + V 0 d ( f ( x ) ) ) .
y ( x ) + ( a x 2 + b x + c ) y ( x ) = 0 ,
y ( t ) + ( a t 2 + C ) y ( t ) = 0 ,
z \rm{d} 2 w \rm{d} z 2 + ( 1 2 z ) \rm{d} w \rm{d} z + ( 1 4 C 4 a ) w ( z ) = 0 ,
z \rm{d} 2 w \rm{d} z 2 + ( B z ) \rm{d} w \rm{d} z + A w ( z ) = 0 ,
{ α ( x ) = exp (  -  z ( x ) / 2 ) M ( A , 1 2 , z ( x ) ) , β ( x ) = exp (  -  z ( x ) / 2 ) [ z ( x ) ] 1 2 M ( A + 1 2 , 3 2 , z ( x ) ) ,
{ α ( x ) = 1 2 z ( x ) α ( x ) + 2 A z ( x ) exp (  -  z ( x ) / 2 ) M ( A + 1 , 3 2 , z ( x ) ) , β ( x ) = 1 z ( x ) 2 z ( x ) z ( x ) β ( x ) + 2 A + 1 3 z ( x ) z ( x ) exp [  -  z ( x ) / 2 ] M ( A + 3 2 , 5 2 , z ( x ) ) .
y ( x ) + ( b x + c ) y ( x ) = 0.
α ( x ) = A i ( w ( x ) ) , β ( x ) = B i ( w ( x ) ) ,
α ( x ) =  -  b 1 / 3 A i ( w ( x ) ) , β ( x ) =  -  b 1 / 3 B i ( w ( x ) ) .
airy ( z ) = A i ( z ) , airy ( 1 , z ) = A i ( z ) ,
airy ( 2 , z ) = B i ( z ) , airy ( 3 , z ) = B i ( z )
y ( x ) + c y ( x ) = 0 ,
α ( x ) = e i c x , β ( x ) = e i c x .