1. Introduction

Fiber Raman amplifier (FRA) is currently considered to be an enabling technology in wide-band transmission systems for its two remarkable merits: low noise figure and flexible gain bandwidth. The increasing research interests and commercial demands of FRA urgently call for efficient simulation models to guide system design.

FRA modeling has been widely discussed in recent years and numerous models are presented for FRA simulation. In 1970s, a simplified numerical model for Raman scattering between pumps and signals as well as amplified spontaneous emission, which represents an early start of FRA research, is propounded by Auyeung and Yariv [1]. In 1999, Kidorf, Rottwitt, Nissov, Ma and Rabarijaona [2] present, for the first time, a detailed model containing almost all physical effects in FRA, which contributes a lot to the research area of FRA simulation. According to the model in [2], Min, Lee, Park [3] and Wang, Fan [4] make some simplification and give two different efficient methods to mitigate the computing time, respectively. And recently, several multi-step methods are applied by Liu, Zhang, Guo [5] and Liu, Lee [6] to enhance the accuracy and stability in the numerical solution of propagation equations for Raman amplification (PERA). But many existent FRA models are still restricted by rigorous assumptions, unacceptable time complexity and/or unsatisfied stability.

Stiffness is a subtle, difficult, and important concept in the numerical solution of ordinary differential equations. It depends on the differential equation, the initial/boundary conditions, and the numerical method. Non-stiff methods can also solve stiff problems, but they generally take a long time to do it. Thus, when we are concerned with a computational version of such problems, stiffness becomes an efficiency issue [10].

In this paper, to the best of our knowledge, the stiffness of PERA is, for the first time, dissected and a novel method for PERA, named backward differentiation formulae (BDF) method, is proposed based on this analysis. In addition, to evaluate its performance, our BDF method is compared with the reported predictor-corrector (PC) method [6], which can effectively improve the accuracy and stability compared with the one-step method and explicit multi-step method.

2. Equation performance analysis

Propagation equations governing the power evolutions of pumps and signals in FRA can be expressed as [2], [7] and [8]

where subscripts ζ and ν are optical frequencies; + and - denote the forward- and backward-propagating optical waves, respectively; α_ν , T, h, k and ε_ν are fiber loss, temperature, Planck’s constant, Boltzmann constant and Rayleigh-backscattering coefficient, respectively; g_ζν is Raman gain coefficient at frequency ν due to pump at frequency ζ (g_ζν defined here has been divided by the fiber effective area A_eff ); Γζ_ν is the polarization factor between frequencies ζ and ν. Δν and Δζ are noise spectral bandwidths for simulation.

Based on Eqs. (1a) and (1b), a Raman-amplified transmission system, which contains n forward propagating waves and m backward propagating waves, can be modeled as

where the terms with subscripts from 1 to n denote the forward-propagating waves, while the ones from n+1 to n+m denote the backward-propagating waves. In addition, to simplify the analysis, we ignore the Rayleigh-backscattering and thermal noise terms in Eq. (1). Hence, here

$F_i(z,\overrightarrow{P})=\{\begin{array}{cc}-{\alpha }_i{P}_i+P_i{\displaystyle \sum _{v_j>v_i}}\frac{g_{ji}}{{\Gamma }_{ji}}·P_j-P_i{\displaystyle \sum _{v_j<v_i}}\frac{v_i}{v_j}·\frac{g_{ij}}{{\Gamma }_{ij}}·P_j& (i=1,2,\cdots ,n)\\ {\alpha }_i{P}_i-P_i{\displaystyle \sum _{v_j>v_i}}\frac{g_{ji}}{{\Gamma }_{ji}}·P_j+P_i{\displaystyle \sum _{v_j<v_i}}\frac{v_i}{v_j}·\frac{g_{ij}}{{\Gamma }_{ij}}·P_j& (i=n+1,n+2,\cdots ,n+m)\end{array}$

By computing the Jacobi matrix of F⃗(z,P⃗) at certain points (z_c ,P⃗(z_c )), the stiffness of Raman-amplified transmission system can be obtained. A detailed introduction of the corresponding knowledge about stiffness theory is contained in the appendix.

To show the stiffness of PERA, a co-pumped transmission system is demonstrated with following configurations: fiber length is 50km; fiber loss is 0.2dB/km for signals and pumps; signal power of each channel is 1mW. There are 21 signal channels from 1530nm to 1546nm with the channel-spacing 0.8nm. The pump powers are 0.5W corresponding to the wavelengths 1445nm and 1475nm.

 

Fig. 1. The stiffness ratio of PERA along the fiber distance.

Download Full Size | PDF

Figure 1 shows that the PERA have obvious and even severe stiffness. In order to maintain the stability, an excessively small discrete step-size must be employed with respect to the smoothness of the exact solution [9]. But it seems unacceptable for simulation because smaller step-size often means larger computing time. Therefore, the stiffness of PERA is expected to be a key issue which prevents the time complexity from decreasing.

In addition, the stability of numerical methods for ordinary differential equations is deeply concerned with the error propagation. In PERA, we assume the rounding error of power P_j to be ε_j . In order to reduce the error, a valuable transformation, which has been used in the previous paper [5] and gains a good effect, can be employed here and converts Eq. (2) into

Then, ε_j becomes the rounding error of ln(P_j ), and it is easy to conclude that the corresponding error of P_j is approximately P_jε_j , which is remarkably less than ε_j because the power P_j is often in the range of [10^-4,5×10^-1] (P_j is with the scale of international standard unit: W).

To verify the above conclusion, the transmission system demonstrated for studying the stiffness of PERA is employed again here and the maximum error of ‖P⃗‖ is calculated with the variation of step-size.

 

Fig. 2. The maximum error of ‖P⃗‖ varied with the iterating step based on model (2) and (3).

Download Full Size | PDF

Figure 2 shows that within the same step-size, model Eq. (3) performs much better than Eq. (2). And in part 4, a FRA is designed based on these two models with comparison, which verifies that model Eq. (3) has a larger tolerance of the iterating step-size.

3. Numerical model

Based on the equation performance analysis, BDF method is employed for its excellent behavior in stiff problem solution [9]. The corresponding knowledge about BDF method is introduced in the appendix.

According to this method, Eq. (3) can be solved as

in which [4], [8]

$F_{\mathit{pi}}^{*}(z_k,\overrightarrow{P})=\left\{-({\alpha }_p+4{\displaystyle \sum _j}{g}_{\mathit{ij}}\frac{{\nu }_i}{{\nu }_j}{N}_{\mathit{nj}})-{\displaystyle \sum _j}{g}_{\mathit{ij}}\frac{{\nu }_i}{{\nu }_j}{P}_{\mathit{sj}}^{+}(z_k)\right\}$ $+{\displaystyle \sum _{{\nu }_n>{\nu }_i}}{g}_{\mathit{ni}}[P_{\mathit{pn}}^{+}(z_k)+P_{\mathit{pn}}^{-}(z_k)]-{\displaystyle \sum _{{\nu }_m<{\nu }_i}}{g}_{\mathit{im}}\frac{{\nu }_i}{{\nu }_m}·[P_{\mathit{pm}}^{+}(z_k)+P_{\mathit{pm}}^{-}(z_k)]\}·\Delta z$ $F_{\mathit{sj}}^{*}(z_k,\overrightarrow{P})=\{-{\alpha }_s+{\displaystyle \sum _i}{g}_{\mathit{ij}}[P_{\mathit{pi}}^{+}(z_k)+P_{\mathit{pi}}^{-}(z_k)]+{\displaystyle \sum _{{\nu }_n>{\nu }_j}}{g}_{\mathit{nj}}·P_{\mathit{sn}}^{+}(z_k)-{\displaystyle \sum _{{\nu }_m<{\nu }_j}}{g}_{\mathit{jm}}\frac{{\nu }_j}{{\nu }_m}·P_{\mathit{sm}}^{+}(z_k)\}·\Delta z$

Here, fiber distance is divided into N elemental amplifier sections for iteration [3] and subscript k denotes the corresponding discrete point; subscripts p and s denote pump and signal, respectively; i and j denote the ith pump and jth signal, respectively; Δz is the step-size between two adjacent discrete points; N_nj =hν_sj Δν is the noise power within bandwidth Δν at frequency ν_sj . It is necessary to point out that Rayleigh-scattering and the temperature dependence are not considered here. In addition, the polarization factor Γ is assumed to be 1.

However, each iteration of Eq. (4) requires estimations of F* _pi (z _k±1,P⃗) and F* _sj (z _k+1, P⃗). These need a set of approximations of $P_{\mathit{\text{pi}}}^{\pm }$ (z _k±1) and $P_{\mathit{\text{sj}}}^{+}$ (z _k+1), which can be calculated by one-step method [4], [8] or explicit multi-step method [5]. Moreover, $P_{\mathit{\text{pi}}}^{+}$ (z1), $P_{\mathit{\text{pi}}}^{+}$ (z ₂), $P_{\mathit{\text{pi}}}^{+}$ (z ₃), P^- _pi (z _N-1), P^-_pi(z _N-2), P^- _pi (z _N-3), $P_{\mathit{\text{sj}}}^{+}$ (z ₁), $P_{\mathit{\text{sj}}}^{+}$ (z ₂) and $P_{\mathit{\text{sj}}}^{+}$ (z ₃), which act as starting values, can be obtained by the semi-analytical method [4], [8] or the traditional Runge-Kutta method.

4. Simulation experiments and discussion on results

In order to estimate the performance and verify the validity of our BDF method for multiple FRA, a simulation experiment is designed. In this simulation, we assumed that the fiber length is 80km and fiber loss is 0.2dB/km. There are 101 channels located from 186THz to 196THz with the channel-spacing 100GHz. The signal power of each channel is 0.5mW, and the pump powers are 198, 145, 132, 71, 59, 47, 40 and 109mW corresponding to the wavelengths 1422.3, 1431.8, 1440.6, 1449.7, 1456.5, 1465, 1475 and 1500nm. In addition, we assumed that the polarization of pumps and signals is completely scrambled along the fiber.

The pump power evolution and net gain of FRA, shown in Figs. 3 and 4, are calculated by the commercial simulation system VPI, our BDF method and the reported PC method [6], respectively. It should be noticed that, to compare our method and PC method fairly, the same step-size is chosen in calculation.

 

Fig. 3. Pump power evolution along the fiber distance obtained from three methods.

Download Full Size | PDF

 

Fig. 4. Net gain of FRA obtained from three methods.

Download Full Size | PDF

According to the results of Tab. 1 based on Figs. 3 and 4, it can be easily seen that, within the same step-size, our method nearly has the same accuracy as PC method does, and the computational time is nearly the same (about 1.3 times) as PC method. Moreover, BDF method exhibits a good precision compared with VPI. The largest difference of pump power is 2.23×10^-3 W, which is negligible, and the largest difference of gain is about 0.28dB, which is difficult to measure experimentally.

Table 1. Comparison of the largest difference of the values calculated by three methods

View Table | View all tables in this article

Furthermore, in order to compare the stability, all of the pump powers are increased to 220mW. When iterating step extends to 1.33km, we get the following results.

 

Fig. 5. Signal power evolution along the fiber calculated by different methods, in which (a) for BDF method based on model (2), (b) for PC method, (c) for BDF method based on model (3) and (d) for VPI.

Download Full Size | PDF

Figure (5(a)) is the power evolution calculated by BDF method based on model (2), Fig. (5(b)) is obtained from PC method, Fig. (5(c)) is calculated by BDF method based on model (3) and Fig. (5(d)) is calculated by VPI. From Fig. 5, we can see that when the pump powers become larger and/or iterating step extends to a certain degree, PC method, compared with VPI, cannot convergent to the right result while BDF method still works, which shows that BDF method has a stronger stability than PC method. Moreover, Fig. (5(a)) illustrates that BDF method based on model (2) also fails to get the right result in this case, which verifies that the transformation from (2) to (3) is useful and contributes a lot to the numerical stability.

In addition, we found that, when the iteration step decreases to 0.5km, PC method can also convergent to the right result. However, its computational time is about two times longer than our BDF method with the iteration step 1.33km.

In summary, the above results of the simulation experiments show that our BDF method, compared with the PC method, can further improve the stability and simultaneously maintain a similar accuracy.

5. Conclusion

We have analyzed clearly the stiffness of PERA and find out some of the main factors which restrict the iterating step. Then a novel simulation method is presented based on this analysis. In addition, numerical experiments are employed to check the validity of our method. The results of the simulation experiments show that our BDF method performs much better in stability than the predictor-corrector method with a similar accuracy, which indicates our method to be a valuable tool for FRA design.

Appendix

It is necessary to point out that, there have been several somewhat different definitions of stiff problem and it has been agreed on that it is not possible to exactly state what it is meant by a stiff problem [9]. Thus, the following theories are based on only one alternative definition presented by Lambert [11].

First, consider a non-homogeneous linear system of ordinary differential equations (ODEs) with constant coefficients [9]

(5)$$\frac{d\overrightarrow{y}}{\mathit{dt}}=A\overrightarrow{y}(t)+\varphi (t),t\in [a,b]$$

where y⃗(t)=(y ₁(t),y ₂(t), …, y_n (t)) ^T is the solution vector, ϕ⃗(t)=(ϕ ₁(t), ϕ ₂(t),…,ϕ_n (t)) ^T is a known portion, t is the independent variable and A∈R ^n×n is the coefficient matrix. Now, assume that A has n distinct eigenvalues λ_j,j=1, 2, …,n. Then

$\overrightarrow{y}(t)={\displaystyle \sum _{j=1}^n}{C}_j{e}^{{\lambda }_{j}t}{\overrightarrow{v}}_j+\overrightarrow{\phi }(t)$

where C ₁,…,C_n are constants, {v⃗ _j } is a basis formed by the eigenvectors of A. As an initial value problem, we assume that Re(λ_j )<0 for all j to make sure that the system is a stable one, which means when t→∞, the solution tends to the particular solution φ⃗(t).

Therefore, φ(t) ${\displaystyle \sum _{j=1}^n}{C}_j{e}^{{\lambda }_{j}t}$ can be interpreted as the steady-state solution and the transient solution, respectively. If we employ a numerical scheme with a bounded region of absolute stability, the step-size is subject to a constraint that depends on the maximum module eigenvalue of A. On the other hand, the greater this module is, the shorter the corresponding time interval is needed to be. We are thus faced with a sort of paradox: the scheme is forced to employ a small integration step-size to track a component of the solution that is virtually flat for large values of t. Generally, stiffness is presented to illustrate the difference of decreasing speed corresponding to each solution portion $C_j{e}^{{\lambda }_{j}t}{\overrightarrow{v}}_j$.

Thus, a linear system Eq. (5) is called a stiff system, if

(6a)$$\mathrm{Re}({\lambda }_j)<0,j=1,2,\cdots ,n$$

(6b)$$s=\underset{1\le j\le n}{max}\mid \mathrm{Re}({\lambda }_j)\mid ⁄\underset{1\le j\le n}{min}\mid \mathrm{Re}({\lambda }_j)\mid \square 1$$

where s is called the stiffness ratio, and the system is regarded as a stiff one if s reaches the magnitude of 10 or above.

Making an extension, to a nonlinear system

(7)$$\frac{d\overrightarrow{y}}{\mathit{dt}}=f(t,\overrightarrow{y}(t)),t\in [a,b]$$

assume that (J(t) _c is the Jacobi matrix of f(t,y⃗(t)) at a certain point (t_c , y⃗(t_c )),

$J(t_c)={\frac{\partial \overrightarrow{f}(t,\overrightarrow{y})}{\partial \overrightarrow{y}}\mid }_{t=t_c}={\left[\begin{array}{cccc}\frac{\partial f_1}{\partial y_1}& \frac{\partial f_1}{\partial y_2}& \cdots & \frac{\partial f_1}{\partial y_n}\\ \frac{\partial f_2}{\partial y_1}& \frac{\partial f_2}{\partial y_2}& \cdots & \frac{\partial f_2}{\partial y_n}\\ ⋮& ⋮& \cdots & ⋮\\ \frac{\partial f_n}{\partial y_1}& \frac{\partial f_n}{\partial y_2}& \cdots & \frac{\partial f_n}{\partial y_n}\end{array}\right]\mid }_{t=t_c}$

and (λ_j =λ_j (t_c ), j=1, 2,…,n is the eigenvalues of (J(t_c ). If a nonlinear system satisfies the condition of (6), then it is regarded as a stiff one and $s(t_c)={\displaystyle \underset{1\le j\le n}{max}}\mid \mathrm{Re}({\lambda }_j)\mid ⁄{\displaystyle \underset{1\le j\le n}{min}}\mid \mathrm{Re}({\lambda }_j)\mid$ is called the local stiffness ratio of (t_c , y⃗(t_c )).

It is valuable to claim that the above definitions are for initial value problems, so the condition Eq. (6a) is necessary to ensure the stability of the systems. While in Raman-amplified transmission system, which is a boundary-value problem, the condition of Eq. (6a), according to our understanding, is unnecessary any more and the stiffness ratio can be defined as

$s={\displaystyle \underset{1\le j\le n}{max}}\mid \mathrm{Re}({\lambda }_j)\mid ⁄{\displaystyle \underset{1\le j\le n}{min}}\mid \mathrm{Re}({\lambda }_j)\mid$

to illustrate the difference of varying speed of each solution portion.

In addition, the so-called BDF method is a kind of implicit multi-step numerical scheme derived from a complementary approach to the one followed for the Adams methods [9]. The corresponding form for Eq. (7) is

$$\overrightarrow{y}(z_{k+1})={\displaystyle \sum _{i=0}^p{q}_i\overrightarrow{y}(z_{k-i})+h{b}_{-1}{\overrightarrow{f}}_{k+1}}$$

where k denotes the discrete point; h is the step-size; a_i and b _-1 are constants, whose values are shown in the following table corresponding to the order of p. The numerical scheme derived in this paper is a three-order (p=3) BDF method and the corresponding coefficients come from the following table.

Table 2. Coefficients of zero-stable BDF methods for p=0,1,…,5.

View Table | View all tables in this article

References and Links

1. J. Auyeung and A. Yariv, “Spontaneous and stimulated Raman scattering in long low loss fibers,” IEEE J. Quantum Electron. 14, 347–352 (1978). [CrossRef]  

2. H. Kidorf, K. Rottwitt, M. Nissov, M. Ma, and E. Rabarijaona, “Pump interactions in a 100-nm bandwidth Raman amplifier,” IEEE Photon. Technol. Lett. 11, 530–532 (1999). [CrossRef]  

3. B. Min, W. J. Lee, and N. Park, “Efficient formulation of Raman amplifier propagation equations with average power analysis,” IEEE Photon. Technol. Lett. 12, 1486–1488 (2000). [CrossRef]  

4. S. Wang and C. Fan, “Generalised attenuation coefficients and a novel simulation model for Raman fibre amplifiers,” IEE Proc.-Optoelectron. 148, 156–159 (2001). [CrossRef]  

5. X. Liu, H. Zhang, and Y. Guo, “A novel method for Raman amplifier propagation equations,” IEEE Photon. Technol. Lett. 15, 392–394 (2003). [CrossRef]  

6. X. Liu and B. Lee, “A fast stable method for Raman amplifier propagation equations,” Opt. Express 11, 2163–2176 (2003), http://www.opticsexpress.org/abstract.cfm?URI=OPEX-11-18-2163. [CrossRef]   [PubMed]  

7. S. Namiki and Y. Emori, “Ultrabroad-band Raman amplifiers pumped and gain-equalized by wavelength-division-multiplexed high-power laser diodes,” IEEE J. Select. Topics Quantum Electron. 7, 3–16 (2001). [CrossRef]  

8. S. Hu, H. Zhang, and Y. Guo, “Simulation model for high-speed wide-bandwidth Raman-amplified WDM system,” in Proc. APOC2003, Wuhan, China, Paper No. 5281–06.

9. A. Quarteroni, R. Sacco, and F. Saleri, Numerical Mathematics (Springer-Verlag, New York, 2000).

10. http://www.mathworks.com/company/newsletter/clevescorner/may03_cleve.shtml.

11. J. D. Lambert, Computational methods in ordinary differential equations (John Wiley & Sons Ltd., London, 1973).