## Abstract

A series of new methods based on the Runge-Kutta (RK) formula are proposed, which not only retain the merit of RK methods, in that the adaptive stepsize is easily implemented, but also dramatically decrease the error under the same conditions. Based on the new methods, an effective shooting algorithm is also proposed. A two-point boundary value problem for the Raman amplifier propagation equations is solved using the proposed algorithm. Our algorithm markedly increases the simulating speed for Raman amplifier propagation equations, as well as improves the accuracy, compared to the traditional algorithm.

© 2003 Optical Society of America

## 1. Introduction

Wavelength-division-multiplexing (WDM) fiber communication systems are expected to span the entire low-loss range in 1.2–1.7 µm. This huge bandwidth will only be accessible if suitable optical amplifiers are available. Then, WDM based on fiber amplifiers provides a transmission platform for optical fiber communication systems and networks [1–8]. In these WDM transmission systems with fiber amplifiers [e.g., erbium-doped fiber amplifiers (EDFA) and Raman amplifiers], which are usually bi-directional pump/signal flow, an analysis of bi-directional pump/signal propagation in fibers is essential. Therefore, the boundary value problems of ordinary differential equations must be solved in the fiber amplifier systems, and a shooting algorithm is typically adopted [5, 6]. However, this method is often divergent and is applied with difficulty by researchers [5].

Although multistep methods require the least number of function evaluations than any other method [9–15], their main drawback is that the largest number of manipulations is needed to perform between steps since the adaptive stepsize calculation must be performed each time the stepsize is changed. Generally, it is inefficient to employ a fixed stepsize when solving differential equations. A typical solution might contain segments with rapid changes where small steps are required, and segments with slow changes where larger steps are acceptable. Methods for automatic step adjustment are often based on an estimate of the local error. Fortunately, Runge-Kutta (RK) methods can be effectively used to implement the variable stepsize problem.

In this paper, we propose a series of new methods based on the RK formula. The proposed methods not only retain the merit of the RK methods, in that the adaptive stepsize is easily implemented, but also decrease the error remarkably under the same conditions. A two-point boundary value problem for the Raman amplifier propagation equations was solved by means of our proposed algorithm.

## 2. New methods

A variety of basic techniques are required for the solutions of initial value problems of ordinary differential equations. An initial value problem can be written as the form,

where ** y** represents a multi-dimensional vector containing elements

**={**

*y**y*

_{1}(

*t*),

*y*

_{2}(

*t*),

*y*

_{3}(

*t*), …} and the initial conditions are given by

Typically, methods for iterating forward in the variable given a single initial point are known as RK methods. RK methods have the general form

In Eq. (3), {*c _{i}*,

*w*,

_{i}*d*} denote all parameters. The development of a specific case entails determining the best possible values for these constants by matching the expansion of this formula with a Taylor series expansion.

_{ij}However, by substituting ** z** to ln(

**) and**

*y***(**

*g**t,*) to

**y****(**

*f**t,*)/

**y****(Here we use bold scripts only for a shortened notation for a set of equations. Hence, hereafter the division or multiplication of bold scripts do not have any meaning of vector calculus.), Eq. (1) is changed to**

*y*According to the operational process of Eq. (3), Eq. (4) can be solved as

In order to compare with the RK method, i.e., Eq. (3), Eq.(5) needs to be changed. Substituting ** y**=exp(

**) into Eq. (5), one can obtain, respectively,**

*z*The exponential change of variable ** y**=exp(

**) is advantageous in amplifier simulations because light amplification or attenuation resembles the exponential pattern with propagation distance. Thus it is intuitively obvious that such a change of variable may lead to faster convergence than some linear or polynomial approximation. Although such a technique was proposed recently [7], it was very difficult to adopt the adaptive stepsize in the work. To distinguish Eq. (6) from Eq. (3), we refer to Eq. (6) as Our-Runge-Kutta (ORK) method. Perhaps the most popular general stepping scheme used in practice is known as the fourth-order RK method [9–17]. The scheme for the fourth-order ORK method is as follows:**

*z*The fourth-order ORK method gives a local truncation error as

Hence, Eq. (8) replaces the local truncation error of the fourth-order RK method, namely

In Eqs. (8) and (9), *C* and *C*’ are constants. To compare error of fourth-order RK method with that of fourth-order ORK method, we specify one ordinary differential equation, i.e.,

The exact solution for Eq. (10) is

Figure 1 illustrates the relative errors between the numerical values and exact values. Figure 1(a) shows that the error expands considerably when *t*>~2, and the fourth-order RK method is then divergent for calculating Eq. (10). On the contrary, Fig. 1(b) shows that the error in calculating Eq. (10) based on the fourth-order ORK method is zero if the round-off error is ignored, and the result is independent of stepsize *h*. All of the results can be explained from Eqs. (8) and (11), i.e., the local truncation error *ε*=0 for each step.

Based on the Taylor expansion, it can be easily proven that the region of absolute stability for the fourth-order ORK method is (-2.785, 0), value of which is the same as the region of the fourth-order RK method.

## 3. Boundary value problems and the shooting algorithm

Many physical applications do not have specified initial conditions, but rather some given boundary conditions. A typical example of such a problem is described by first-order equations containing two elements, namely,

on *t*∊[*a, b*] with the general boundary conditions

Thus the solution is defined over a specific interval and must satisfy the relations of Eq. (13) at the end points of the interval. This type of problem is fundamentally more difficult than the initial value problem because we cannot directly apply any of the numerical integration techniques for boundary value problems and do not even have a starting point to attempt a numerical integration. The best that can be done in this situation is to make an initial guess that can then be improved by some iterative process.

To solve Eq. (12) with the constraint conditions of Eq. (13), two types of initial guesses are available. One is assumed to be *y*
_{2}(*a*)=*A*
_{1}, and the calculating process is from *t*=*a* to *t*=*b*. The other is assumed to be *y*
_{1}(*b*)=*B*, and the calculating process is from *t*=*b* to *t*=*a*. However, the region of absolute stability for the fourth-order ORK/RK method is (-2.785, 0). For some specific cases of Eq. (12), therefore, the first type of initial guess is convergent, but the second is divergent, *vice versa*. In the following, we provide a specific case of Eq. (12). That is, *f*
_{1}=*y*
_{2} and *f*
_{2}=*y*
_{1}, *a*=0.5 and *b*=2, and *y*
_{1}(*a*)=exp(0.5)+exp(-0.5) and *y*
_{2}(*b*)=exp(2)-exp(-2). Their exact solutions are *y*
_{1}(*t*)=exp(*t*)+exp(*-t*) and *y*
_{2}(*t*)=exp(*t*)-exp(*-t*). We employ the first type of initial guess to numerically calculate this specific case.

This case can be solved using a self-consistent algorithm to search for the appropriate value of *A*
_{1}, which satisfies the original problem. The basic algorithm is as follows:

- Solve the differential equation using Eq. (7) with the initial conditions
*y*_{1}(0.5)=exp(0.5)+exp(-0.5) and*y*_{2}(0.5)=*A*_{1}. - Evaluate the solution
*y*_{2}(*t*) at*t*=2, compare this value with the target value of*y*_{2}(2)=exp(2)-exp(-2), and obtain the relative error. - Adjust the value of
*A*_{1}by an interpolating method until a desired level of tolerance and accuracy is achieved. - Once the specified accuracy has been achieved, the numerical solution is complete.

To obtain a better estimate *A*
_{j+1} of the initial guess of the (*j*+1)-th iteration, in step (3), the interpolating method employs the following expression

where *A _{j}* and

*β*are the initial guess and the solution

_{j}*y*

_{2}(

*b*) in the

*j*-th iteration,

*A*

_{j-1}and

*β*

_{j-1}are the corresponding values in the (

*j*-1)-th iteration.

Details of the procedure are shown in the Appendix, where the MATLAB source codes for above case of the shooting algorithm are given. Details of the calculating procedure are demonstrated in Fig. 2, where (a) for *y*
_{1}(*t*), (b) for *y*
_{2}(*t*), and (c) for the relative error of *y*
_{1}(*t*) and *y*
_{2}(*t*) in the last iteration. In the legends of Fig. 2, the 1st-, 2nd-, and 3rd-iteration represent the first-, second-, and third-iteration in the shooting algorithm, respectively, and *exact* indicates the exact solution of *y*
_{1}(*t*) and *y*
_{2}(*t*).

Figure 2 shows that, based on the fourth-order ORK method, only a three-times iteration can satisfy the constraint conditions, and all calculating points have the relative error of <10^{-5}. The numerical results also show that it is *divergent* if the fourth-order RK method is adopted.

Based on above shooting algorithm for calculating from *t*=*a* to *t*=*b*, however, some cases are divergent. But, it is convergent when calculating from *t*=*b* to *t*=*a*. The next section will describe this in detail.

## 4. Applications for Raman amplifier propagation equations

In calculating the Raman amplifier gain profile, the major influence is pump-to-pump, signal-to-signal and pump-to-signal interactions, as well as the attenuation. In the steady state, the coupled equation is described by the following [18]

where *P _{i}*,

*v*and

_{i}*α*are the power, frequency and attenuation coefficient for the

_{i}*i*th wave, respectively;

*A*is the effective area of optical fiber; the factor of Γ accounts for randomization polarization effects, the value of which lies between 1 and 2; and

_{eff}*g*(

_{R}*v*) is the Raman gain coefficient from wave

_{j}-v_{i}*j*to

*i*. The frequency ratio

*v*/

_{i}*v*describes vibrational losses. The minus and plus signs on the left hand side describe the backward-propagating pump waves and forward-propagating signal waves, respectively. The frequencies

_{j}*v*are numerated in the decreasing order of frequency (

_{i}*i*=1, 2, …,

*m*). The terms from

*j*=1 to

*j*=

*i*-1 and from

*j*=

*i*+1 to

*j*=

*m*cause amplification and attenuation of the channel at frequency

*v*.

_{i}Apparently, Eq. (15) is a boundary value problem for backward-propagating pumps. The shooting algorithm must be adopted to solve Eq. (15). However, it is divergent when the power of the backward-propagating pump exceeds a certain value (e.g., 400 mW). By employing Eq. (7) and starting *z*=*L* instead of *z*=0, the stability of the shooting algorithm is remarkably improved.

Using Eqs. (7) and (15), a distributed multi-pump Raman amplifier is simulated. In the numerical calculations, we assumed that *A _{eff}*=80×10

^{-12}, Γ=2, the length of fiber

*L*=80 km,

*α*=0.2 and 0.35 dB/km for signals and pumps, respectively. There are 19 signal channels spaced at 200 GHz/channel, from 188.85 to 192.45 THz. The power of each channel is 1 mW, and powers of the two backward-propagating pumps are 500 mW, corresponding to wavelengths of 1455 and 1475 nm. The initial guesses of all signals at

*z*=

*L*are assumed to be 1 mW [see Fig. 4(b)].

*g*is illustrated in Fig. 3 [2].

_{R}Figure 4 shows the numerical results and the iterating procedure for shooting algorithm starting from *z*=*L*, where (a) and (b) for the procedure of the pump and signal, respectively, and (c) for all of the pumps and signals in the fourth iteration. In the legends, the 1st-, 2nd-, 3rd-, and 4th-iteration represent the first-, second-, third-, and fourth-iteration in shooting algorithm, respectively. From Fig. 4(a) and (b), it can be seen that only four iterations of the shooting algorithm can reach the constraint conditions with a relative error of <10^{-3}. Figure 4(c) shows the strong stimulated Raman scattering interactions of pump-to-pump, signal-to-signal and pump-to-signal. The numerical results also show that, when calculating from *z*=0 to *z*=*L*, it is divergent when each pump power exceeds 400 mW under our simulation parameters. To the contrary, it is convergent if starting from *z*=*L*. However, when each pump power exceeds 1 W, our proposed shooting algorithm is also divergent. The reason for this is that it is beyond the absolute stability region of the fourth-order ORK method, i.e., (-2.785, 0). Although this situation for each pump power >1 W in a Raman amplifier is infrequent, the divergent problem can still be solved by the multi-step method [19].

## 5. Conclusions

Based on the RK formula, a series of new methods for ordinary differential equations are proposed, for the first time to the authors’ knowledge. The proposed methods not only retain the merit of RK methods, in that the adaptive stepsize is easily implemented, but also decrease the error considerably under the same conditions. Using the new methods, an effective shooting algorithm is also proposed. A two-point boundary value problem for the Raman amplifier propagation equations can be solved by means of our proposed algorithm. The numerical results show that our proposed shooting algorithm can markedly increase the simulation speed for Raman amplifier propagation equations, and improve the calculation accuracy, in comparison with traditional algorithm. Simulated calculations also demonstrate that our algorithm remains stable when each pump power is <1 W under our simulating parameters, the value of which is more than twice that of the traditional algorithm. Our proposed algorithm has important applications in fiber amplifiers such as the Raman amplifier and EDFA.

## Appendix

function Ourshoot

t0=0.5; tf=2; N=250; t=linspace(t0,tf,N);

y0=ExaSol(t0); % initial value for first equation

A1=3; y0(2)=A1; % initial guesses for second equation

yf0=ExaSol(tf); % get final point value

EPSShooting=10^-5; close all; figure; figure;

col=str2mat(‘b:’,‘g:’,‘c:’,‘m:’,‘k:’,‘b--’,‘g--’,‘c--’);

ya=ExaSol(t); figure(1); plot(t,ya(:,1),‘r-‘); hold on;

figure(2); plot(t,ya(:,2),‘r-’); hold on;

for i=1:50, % maximum loop times

[y,t]=Ourrk4([t0,tf],y0,N);%y0(2) is modified each loop

if i==1,

yf1=y(end,2); %get last point value for second equation

flagEPS=abs((y(end,:)-yf0)./yf0); % relative error

A2=yf0(2).*A1./yf1; y0(2)=A2; % y0(2) is modified

else

yf2=y(end,2); %get last point value for second equation

flagEPS=abs((y(end,:)-yf0)./yf0); % relative error

A3=A1+(A2-A1)*(yf0(2)-yf1)/(yf2-yf1);

A1=A2; A2=A3; yf1=yf2; y0(2)=A2; % y0(2) is modified

end

figure(1); plot(t,y(:,1),col(i,:)); hold on;

figure(2); plot(t,y(:,2),col(i,:)); hold on;

if all(flagEPS<EPSShooting)=1, break; end

end

figure; semilogy(t,abs((y-ya)./ya));

function [y,x]=Ourrk4(x0,y0,n)

x=linspace(x0(1),x0(end),n); x=x’; h=(x(2)-x(1)); y(1,:)=y0;

for i=1:(n-1),

k1=dydxy(x(i),y(i,:)); %dy/dx/y

k2=dydxy(x(i)+h/2,y(i,:)+h*y(i,:).*k1/2);

k3=dydxy(x(i)+h/2,y(i,:)+h*y(i,:).*k2/2);

k4=dydxy(x(i)+h,y(i,:)+h*y(i,:).*k3);

y(i+1,:)=y(i,:).*exp(h*(k1+2*k2+2*k3+k4)/6);

end

function fy=dydxy(x,y) %y‘/y=f(x,y)/y

fy=[y(2)/y(1), y(1)/y(2)];

function y=ExaSol(x) % exact solution

y(1,:)=exp(x)+exp(-x); y(2,:)=exp(x)-exp(-x); y=y’;

## Acknowledgment

This work was supported by the Ministry of Education of Korea through the Brain Korea 21 (Information Technology) Program.

## References and links

**1. **T. Mizuochi, K. Kinjo, S. Kajiya, T. Tokura, and K. Motoshima, “Bidirectional unrepeatered 43 Gb/s WDM transmission with C/L band-separated Raman amplification,” J. Lightwave Technol. **20**, 2079–2085 (2002). [CrossRef]

**2. **E. M. Dianov, “Advances in Raman fibers,” J. Lightwave Technol. **20**, 1457–1462 (2002). [CrossRef]

**3. **M. Karásek and M. Menif, “Channel addition/removal response in Raman fiber amplifiers: modeling and experimentation,” J. Lightwave Technol. **20**, 1680–1687 (2002). [CrossRef]

**4. **N. Kikuchi, K. K. Wong, K. Uesaka, K. Shimizu, S. Yam, E. S. Hu, M. Marhic, and L. G. Kazovsky, “Novel in-service wavelength-band upgrade scheme for fiber Raman amplifier,” IEEE Photon. Technol. Lett. **15**, 27–29 (2003). [CrossRef]

**5. **X. M. Liu and B. Lee, “A fast and stable method for Raman amplifier propagation equations,” (to be submitted)

**6. **Q. H. Mao, J. S. Wang, X. H. Sun, and M. D. Zhang, “A theoretical analysis of amplification characteristics of bi-directional erbium-doped fiber amplifiers with single erbium-doped fiber,” Opt. Commun. **159**, 149–157 (1999). [CrossRef]

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

**8. **A. Carena, V. Curri, and P. Poggiolini, “On the optimization of hybrid Raman/erbium-doped fiber amplifiers,” IEEE Photon. Technol. Lett. **13**, 1170–1172 (2001). [CrossRef]

**9. **U. M. Ascher and L. R. Petzold, *Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations* (SIAM, Philadelphia, 1998). [CrossRef]

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

**11. **B. F. Plybon, *An Introduction to Applied Numerical Analysis* (PWS-KENT Publishing Company, Boston, 1992), pp. 428–441.

**12. **J. H. Mathews, *Numerical Methods for Mathematics, Science, and Engineering* (Second edition, Prentice Hall, New Jersey, 1992), pp. 464–475.

**13. **U. M. Ascher, R. M. M. Mattheij, and R. D. Russell, *Numerical Solution of Boundary Value Problems for Ordinary Differential Equations* (SIAM, Philadelphia, 1995). [CrossRef]

**14. **P. B. Bailey, L. F. Shampine, and P. E. Waltman, *Nonlinear Two Point Boundary Value Problems* (Academic Press, New York, 1968).

**15. **L. F. Shampine, *Numerical Solution of Ordinary Differential Equations* (Chapman & Hall, New York, 1994).

**16. **W. H. Enright and P. H. Muir, “Runge-Kutta software with defect control for boundary value ODEs,” SIAM J. Sci. Comput. **17**, 479–497 (1996). [CrossRef]

**17. **R. Weiss, “The application of implicit Runge-Kutta and collocation methods to boundary value problems,” Math. Comp. **28**, 449–464 (1974).

**18. **V. E. Perlin and H. G. Winful, “Optimal design of flat-gain wide-band fiber Raman amplifiers,” J. Lightwave Technol. **20**, 250–254 (2002). [CrossRef]

**19. **X. M. Liu and B. Lee, “A series of fast and accurate algorithms and their applications in fiber transmission systems,” (submitted to IEEE J. Quantum Electron.)