On the Numerical Solution of Diffusion Problem with Singular Source Terms

Partial differential equations with singular (point) source terms arise in many different scientific and engineering applications. Singular means that within the spatial domain the source is defined by a Dirac delta function. Our interest is this type of problems is motivated by mathematical modelling of forecasting and development of new gas reservoirs [1, 2, 3, 4]. Solutions of the problems having singular source terms generally have lack of smoothness, which is an obstacle for standard numerical methods. Therefore, solving these type of problems numerically requires a great deal of attention [5, 6, 7]. In this paper we discuss the numerical solution of initial-boundary value problem with singular source terms 8>><>>: ut = D uxx +k1d (x−a1)+k2d (x−a2), 0 < x < 1, t > 0, 0 < a1 < a2 < 1, u(t,0) = uL, u(t,1) = uR, t  0, u(0,x) =j(x), 0  x  1, (1) where d (x) is a Dirac delta function. We follow the standard finite volume approach based on the integral form of (1). We consider this approach more natural than the finite difference one directly based on the differential form, since for the integral form the treatment of the Dirac delta function expression is mathematically clear. For ease of presentation, we assume that there are only two source terms. The presented material is extendable to the case with more than two source terms. Finally, this study can be readily extended to the case with time-dependent source terms.