Numerical Solution of a Nonlinear Integro-Differential Equation

An algorithm for the numerical solution of a nonlinear integro-differential equation arising in the single-species annihilation reaction $A + A \rightarrow\varnothing$ modeling is discussed. Finite difference method together with the linear approximation of the unknown function is considered. For divergent integrals presented in the equation for dimension $d=2$ a regularization is used. Some numerical results are presented.


Introduction
The irreversible annihilation reaction A + A → ∅ is a fundamental model of non-equilibrium physics. The reacting A particles are assumed to perform chaotic motion due to diffusion or some external advection field such as atmospheric eddy [1]. Many reactions of this type are observed in diverse chemical, biological or physical systems [2,3].
In [1] the advection of reactive scalar using random velocity field generated by the stochastic Navier-Stokes equation, which is used for production a velocity field corresponding to thermal fluctuations [4,5] and a turbulent velocity field with the Kolmogorov scaling behavior [6] is studied by three of the present authors, and the integro-differential equation for the number density is derived. No influence of the reactant on the velocity field itself is assumed.
In this paper we present some initial experiences with the numerical solution of the integrodifferential equation mentioned above.

Problem Formulation
In [1] an integro-differential equation (72) for the mean number density a(t) of chemically active molecules in anomalous kinetics of single-species annihilation reaction A + A → ∅ which corresponds to D = uν, ∆ = 0, and Z 4 = 1.

Case d = 2
For d = 2 the singularity in the integral on the right side of Eq.
(2) at t = t is divergent. This is a consequence of the UV divergences in the model above the critical dimension d c = 2, and near the critical dimension is remedied by the UV renormalization of the model [1]. In this paper we use another approach to overcome this problem. We will use the following regularization: We rewrite Eq. (3) to the form da dt Let us consider the Eq. (4). We will use the difference method to solve it numerically. We will consider time discretization with the time step ∆t: For the discretization of the left side we will use two formula -the first and the second order finite difference approximations: For the right-side integral approximation we will use piecewise linear approximation of the function a(u): Integrals in the right side of Eq. (7) we calculate analytically using in the following way Putting approximations Eq. (6), Eq. (7), and Eq. (9) into Eq. (4) we arrive at the quadratic equations Eq. (11) and Eq. (12) in more explicit form with respect to a k : or in the standard form For the first step we have  Figure 1. Logarithms of the number density a(t) together with its fit (above), and the fitting error (below) or in the standard form The algorithm consists of the successive calculation of the values a k , k = 1, 2, . . . , K solving the quadratic equation Eq. (12) at the first step using the value a 0 = a(0), and further solving Eqs. (11) using previously determined values a i , i = 0, . . . , k − 1.

Numerical results
Below the results for a(0) = 5000, λ = 0.1, D = 0.4 and = 0.001 are presented. All calculations are done with the uniform time step ∆t.    Table 1 compares the results at selected points t for the "increasing precision" for the decreasing values ∆t from 0.01 to 0.0003125. Even for the step ∆t = 1/3200 = 0.0003125 the results may not be considered to have the sufficient precision. If the step number is increasing two times, the calculation time is increasing four times.
If we suppose, that the difference between the numerical value a ∆t (t) for the step ∆t and the "exact" value a(t) has the order p, e.g., where both C and p depends on t, then using the successive approximations for different steps ∆t, ∆t/2, and ∆t/4 we get p(t) ≈ log 2 a ∆t (t) − a ∆t/2 (t) a ∆t/2 (t) − a ∆t/4 (t) .
It seems to be reasonable to study analytically the behavior of the function a(t) for the small values of t, and start the numerical computation from some point t > 0. Also non-uniform grid could be considered.

Conclusions
Numerical results presented above could give us some basic imagination about the behavior of the number density function a(t). However, further improvements of the algorithm are necessary.
It will be also interesting to try to solve the renormalized integro-differential equation where ∆ = (d − 2)/2 and γ 0.57721 is Euler's constant, presented in [1], and compare the results. Another possibility is to study a behavior of the solution of the problem (1) for d → 2 − .