Numerical simulations of stochastic biochemical oxygen demand equations via RBF method

ABSTRACT This paper proposes a new approach to numerical simulations of stochastic biochemical oxygen demand equations. The analysis of biochemical oxygen demand is necessary in evaluating the effects of water pollution. Here, we consider radial basis functions. This approximation process can also be interpreted as a simple kind of neural network. An illustrative example is included to demonstrate the validity and applicability of the approach. The numerical experiments show that the method performs well.


Introduction
White noise is a widespread and common phenomenon in many engineering, biological, economics and chemical models including some form of error prediction. In these types of applications, the deterministic model is often already available and stochastic models are an active area of research.
In recent decades, mathematical models describing wastewater treatment have been popular. Biochemical Oxygen Demand (BOD) is a measure of the dissolved oxygen consumed by microorganisms during the oxidation of reduced substances in waters and wastes. The analysis and forecast of BOD are vital and essential in assessing the effects of water pollution. The purpose of this paper is to extend the Radial Basis Functions (RBFs) method on stochastic models for BOD in a more stochastic and therefore more realistic approach. The final model will be used to predict the amount of BOD at any point on the stream [1][2][3][4].
The biochemical oxygen demand model presented in this paper is often used in water quality determination of river and estuarine systems and is an extension of the so-called Streeter Phelps model. In this model, the concentration levels of Carbonaceous Biochemical Oxygen Demand (CBOD), Dissolved Oxygen (DO), and Nitrogenous Biochemical Oxygen Demand (NBOD) are related to each other by a set of differential equations. Although the concentrations are time dependent, the purpose behind the model is often to monitor the concentration levels within a fixed volume of water, flowing downstream a river. An underlying assumption is that the velocity of the water flow is constant and thus time and distance are linearly related [5][6][7]. This paper establishes a direct method for solving a system of stochastic differential equations via a set of RBFs. The approximation using RBF is an extremely powerful method to interpolate functions. Sums of radial basis functions are typically used to approximate given functions. The RBF methodology was first introduced in 1971 and has been variously studied [8][9][10]. Recently, these functions have also been used in the numerical solution of stochastic equations. They have also been used to solve equations in financial mathematics [11][12][13][14]. This method has been applied for solving different kinds of equations, problems and models such as PDEs, IEs and ODEs in physics, biology and neural network etc.
This paper is organised as follows. The next section reviews radial basis functions and the approach in solving stochastic Volterra integral equations. Section 3 introduces a model that we want to solve in this paper and our numerical findings are reported. Section 4 discusses the error of the RBFs method for the linear integral equations. Finally, Section 5 is the conclusions.

RBF
The goal of this section is to recall notations and definitions of RBFs that are used for approximating the solutions of BOD equations. This section also describes applying the method for solving stochastic integral equations. Definition 2.1. A function ΦðtÞ : R n ! R is called radial if there exists a one variable function φ : ð0; 1Þ ! R such that ΦðtÞ ¼ φð t k kÞ, where : k k is the Euclidean norm. A radial basis function φðrÞ is a univariate continuous real-valued function that depends on the distance from the origin (or any other fixed centre point).
1Þ ! R is defined as follows: Where, φ i : R n ! R and t i for i ¼ 1; 2; . . . ; N is the centre of RBFs. Additionally,t; t i are belong to R n .
RBFs are mostly identified on the basis of smoothness. Some functions are infinitely smooth and some are piecewise smooth. Gaussian Function (GS), Multiquadric (MQ), Inverse Multi-quadric (IMQ) and Inverse quadric (IQ) are some examples of infinitely smooth RBFs whereas Thin Plate Spline (TPS) and Linear radial function (LR) are piecewise smooth RBFs. For infinitely smooth RBFs, there exists a free parameter called the shape parameter that controls the shape of RBF. The RBF become flat if the shape parameter is closer to 0.
Let t 1 ; t 2 ; . . . ; t N f g � R n be a given set of distinct nodal points. To approximate a function yðtÞ using the radial function ΦðtÞ ¼ φð t k kÞ we can give the following linear combination: where the coefficients λ 1 ; λ 2 ; . . . ; λ N are determined by the interpolation conditions Therefore, the solution of the interpolation problem, based on the extended expansion (1), reduces to the solution of a system of linear equations of the matrix form where the pieces are given by

Description of method on stochastic Volterra integral equation
The nonlinear stochastic Volterra integral equation takes the following form: where u t ð Þ; u 0 t ð Þ; k 1 s; t ð Þ and k 2 s; t ð Þ; for s; t 2 0; T ½ Þ; are the stochastic processes defined on the same probability space Ω; F; P ð Þ with a filtration fF t ; t � 0} that is increasing and right continuous and F 0 contains all P-null sets. u t ð Þ is unknown random function and B t ð Þ is a standard Brownian motion defined on the probability space and ð t 0 k 2 ðs; tÞuðsÞdBðsÞ is the Itˆo integral. Let us approximate the function u t ð Þ in terms of radial basis functions,ϕ t ð Þ; as follows: uðtÞ ' Then, from substituting Equation (3) into Equation (2) we have Substituting the collocation points t j N j¼1 into Eq.(4), we obtain: In the above equation, we let s ¼ λ i φðk s À t i kÞdBðsÞ; t 2 ½0; T�; j ¼ 1; . . . ; N: For computing or approximating the first integral, several methods can be used. Now, by applying Legendre-Gauss-Lobatto integration formula, we approximate the first integral of Equation (6) as follows λ i φðk s À t i kÞdBðsÞ; t 2 ½0; T�; j ¼ 1; . . . ; N: where w k ; x k ; k ¼ 1; . . . ; N: are weights and nodes of the integration rule, respectively.
For approximating the Itˆo integral, let 0 may be approximated by the Riemann sum where the discrete points s j ¼ jδt: Indeed, the integral may be defined by taking the limit δt ! 0 in (8). In a similar way, we may consider a sum of the form X NÀ 1 j¼0 f ðs j ÞðBðs jþ1 Þ À Bðs j ÞÞ; Which, by analogy with (8), may be regarded as an approximation to a stochastic integral f ðs j ÞðBðs jþ1 Þ À Bðs j ÞÞ; Here, we are integrating f with respect to Brownian motion. In Equation (6) for Itˆo integral, let 0 ¼ s 0 < s 1 < . . . < s m ¼ t 0 . So, achieved λ i φðk s k À t i kÞ Bðs kþ1 Þ À Bðs k Þ ½ �; t 2 ½0; T�; j ¼ 1; . . . ; N: We have a system of equations that can be solved by mathematical software for the unknowns vector λ. By computing that, unknown function u t ð Þ can be approximated.

Stochastic biochemical oxygen demand model
The biochemical oxygen demand (BOD) model shows the concentration of oxygen needed by aerobic microorganisms. BOD is essentially a measure of the concentration of oxygen required to remove waste organic matter from water in the process of decomposition by aerobic bacteria. The water quality in bodies of water is generally measured by BOD and oxygen concentration. Let b, o and n be the concentration levels, measured in mg/l, for Carbonaceous Biochemical Oxygen Demand (CBOD), Dissolved Oxygen (DO) and Nitrogenous Biochemical Oxygen Demand (NBOD) respectively [6,15]. The CBOD, DO and NBOD are defined by the following differential equation where k b ¼ k c þ k 3 ,s � 2 ¼ k 2 d s þ p 0 À r 0 þ s 2 , k c is the reaction rate coefficient, k 2 is the reaeration rate coefficient, k 3 is the sedimentation and adsorption loss rate for CBO, k n is the decay rate of NBOD, p 0 is photosynthesis of oxygen, r 0 is the respiration of oxygen, d s is the saturation concentration of oxygen, s 1 is the independent source for CBOD, s 2 is the independent source for OD, and s 3 is the independent source for NBOD. As for many environmental models, BOD is also subject to various uncertainties. These uncertainties can be incorporated into the model by adding white noise processes to each of the three equations of the model. For this equation assumes that we can model uncertainties in the source terms (s 1 , s 2 , s 3 ) by adding a random noise factor to each of three equations. The resulting linear SDE is elaborated by Nsengiyumva [15]: where B i (t) are independent standard Brownian motions and σ b , σ o , and σ n are diffusion coefficients in the noises of CBOD, OD and NBOD respectively and represent the intensities of B i (t), i= 1, 2, 3. Usually stochasticity is introduced into the model by making a simple assumption about the random nature of coefficient. If, we can assume that the coefficients k b , k 2 and k n have random nature, the equation can be written as follows: For more information about this model see Berg et al [6]. By integration of Equation (13) the method presented in Section 2 is applied for the following system of equations: We approximate the functions based on RBF method and after solving the system, b t ð Þ; o t ð Þ and n t ð Þ are approximated. Table 1 defines the parameter values used in the simulations. In order to conform the results above, initial value b 0 ð Þ ¼ 15; o 0 ð Þ ¼ 8:5 and n 0 ð Þ ¼ 5 are chosen and parameters corresponding to Table 1. Then b t ð Þ; o t ð Þ and n t ð Þ are shown in the Figure 1.

Conclusion
Some systems of stochastic differential equations cannot be solved analytically; in this paper we present a new technique for solving them numerically. We have used Radial Basis Functions for approximating the solution of the BOD equations. Such nonlinear models occur in financial epidemiology, and biology. For solving the Riemann integral part, the Legendre-Gauss-Lobatto integration formula was applied, then the problem was discretized. The proposed method reduces an integral equation to a system of equations. Implementation of our method is easy and the accuracy has been verified by test examples. Efficiency of this method and a good degree of accuracy was confirmed by a numerical example.