UNCONDITIONALLY POSITIVE FINITE DIFFERENCE AND STANDARD EXPLICIT FINITE DIFFERENCE SCHEMES FOR POWER FLOW EQUATION

Power flow equation for step-index glass optical fiber was solved using recently reported unconditionally-positive finite difference (UPFD) scheme. Solution obtained using UPFD scheme was compared with solution obtained using standard explicit finite difference (EFD) scheme. For accuracy testing both schemes were compared with analytical solution for steady state distribution of given fiber. The advantage of UPFD is reflected in stability of the scheme regardless of discretization step taken. Nevertheless EFD scheme has better concurrence with analytical solution than UPFD. This is due to the additional truncation-error terms in the approximations of the first and second derivatives with respect to θ.


INTRODUCTION
Partial  differential equations are used to model a wide variety of problems in natural sciences. Usually those are physical, chemical, biological and engineering problems such as transmission of light through optical fibers, heat transfer, transport and reaction of chemical species, solid state physics, adsorption of pollutants in soil and diffusion of radon and neutrons (Bear et al., 2007;Hetrick, 1971;Kevkic et al., 2019;Murray, 2002;Petrović, 2017;Shih, 1984). Often those equations are parabolic differential equations with no, or limited, analytical solution and their solvation requires various numerical techniques Savović et al., 2009;Savović et al., 2012;Savović et al., 2013a;Savović et al., 2013b). In case of three-dimensional problems solutions are often obtained by the finite element method (Urošević et al., 2003). If problems are one-dimensional they are much easier solved by finite difference method (Savović et al., 2013b). Both methods have advantages and limitations, and the choice between two methods is mostly conditioned by geometry of the problem. For neither does the standard finite discretization explicitly constrain the solutions to positive values, which may lead to numerical instabilities and oscillations of the solution. Moreover there is no unanimous opinion of choice between different finite difference methods for diverse applications.
With the development of computers most algorithms for numerical methods were implemented in a variety of programming languages. In the beginning, during 1970s and 1980s, implicit finite difference methods (IFDMs) were generally first choice. Being often unconditionally stable, the IFDM allows larger step lengths. Despite, this does not increase IFDM's computational efficiency because extremely large matrices must be manipulated at each calculation step. With further advancement of the computers this trend has been changing, shifting the emphasis to explicit finite difference methods (EFDMs). Although EFDM is not unconditionally stable we found that it is also simpler in addition to being computationally more efficient (Savović et al., 2009;Savović et al., 2012;Simovic et al., 2014).
The constant development of numerical methods in search of the best possible solution has led to the proposition of new solutions. Several authors propose new numerical solution of the parabolic differential equations that guarantees positivity of the solutions and that is independent both of the step size in z direction and mesh size (Chen-Charpentier et al., 2013;Liu et al., 2010;Quang et al., 2006). The method works with reaction terms that are the sum of a positive function and a negative function of the unknowneither or both may be zero. It is applicable to both, problems where either advection or diffusion dominates. In this work, a recently reported UPFD scheme and a standard EFD scheme for solving power flow equation are compared to analytical solution for steady-state distribution for glass optical fiber CGW-CGE-68 (Drljaca, 2011).

THEORETICAL PART
Power distribution in optical fiber can be determined by using one of the three methods: ray-tracing method, wave approach and Gloge`s power flow equation. The simplest way that can describe evolution of the power in the multimode optical fibers and account for all important characteristics of the fiber is Gloge`s power flow equation. The Gloge's power flow equation is given as (Gloge, 1972): where P(,z) is the angular power distribution, z is distance from the input end of the fiber,  is the propagation angle with respect to the core axis, D is the coupling coefficient assumed constant (Garito et al. 1998;Savovic et al., 2015) and ( (1) becomes (Djordjevich et al., 2000): Numerical solution of equation (2) is obtained by EFDM, used in our previous works (Drljaca, 2011;Savović et al. 2013a;Savovic et al., 2015) and UPFDM proposed by other authors (Chen-Charpentier et al., 2013;Liu et al., 2010;Quang et al., 2006.]. The results thus obtained are compared to the analytical solution of steady-state distribution for the given glass optical fiber: where J 0 is the Bessel function of the first kind and zero order

Standard Explicit Finite Difference Method
After applying EFD scheme (Anderson, 1995): Truncation error of the scheme is ()     Oz . The method (7) proposed by Chen-Charpentier and Kojouharov (Chen-Charpentier et al., 2013) is simultaneously explicit and unconditionally positive, but it is not unconditionally consistent. There are extra truncation error terms because the approximations of the first and second derivatives with respect to θ are evaluated at different lengths. One way to reduce this error is to choose the length step depending on the mesh size so that the inconsistent terms approach zero when the mesh is refined. Second approach is to incorporate these terms into the numerical scheme to achieve a consistent approximation of the original partial differential equation (Chen-Charpentier et al., 2013).

RESULTS
In order to present obtained far-field intensity patterns for numerical calculations we used same fiber that was previously used in different experiments and numerical simulations -CGW-CGE-68 glass optical fiber. The fiber has critical angle θ c = 7.26• measured inside the fiber and θ c = 10.6 o measured in air. Coupling coefficient for this fiber is D=7.9×10 -7 rad 2 /m at room temperature and was used for numerical calculations (Drljaca, 2011). Attenuation coefficient is calculated as 0  =2.86×10 -4 m -1 .   (Drljaca, 2011). In Figure 2. numerical solutions obtained by UPFD method are presented. As could be seen from Fig. 1 and Fig. 2 that both EFD and UPFD methods have a good agreement with analytical solution of equation (2). However it could be seen that UPFD scheme is less accurate than EFD scheme when compared to analytical solution. This is due to additional truncation error that occurs in UPFD scheme.

CONCLUSION
Standard explicit finite difference (EFD) scheme and recently reported unconditionally positive finite difference (UPFD) scheme for solving parabolic differential equations are compared to previously obtained analytical solution of power flow equation. We have shown that this scheme is less accurate, even though UPFD scheme guarantees the positivity of the solutions for arbitrary step sizes, and is stable. The reason for this is additional truncation error in the approximations to the first and second derivatives with respect to θ, which are evaluated at different lengths, which is contained in the UPFD scheme. If the accuracy of the numerical scheme over stability is needed we propose usage of EFD scheme before UPFD scheme.