An adaptive Huber method for nonlinear systems of Volterra integral equations with weakly singular kernels and solutions

Accepted Manuscript An adaptive huber method for nonlinear systems of volterra integral equations with weakly singular kernels and solutions L.K. Bieniasz PII: DOI: Reference:

S0377-0427(17)30188-7 http://dx.doi.org/10.1016/j.cam.2017.04.018 CAM 11097

To appear in:

Journal of Computational and Applied Mathematics

Received date: 7 February 2017 Please cite this article as: L.K. Bieniasz, An adaptive huber method for nonlinear systems of volterra integral equations with weakly singular kernels and solutions, Journal of Computational and Applied Mathematics (2017), http://dx.doi.org/10.1016/j.cam.2017.04.018 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An Adaptive Huber Method for Nonlinear Systems of Volterra Integral Equations with Weakly Singular Kernels and Solutions L.K. Bieniasz∗ Faculty of Physics, Mathematics, and Computer Science, Cracow University of Technology, ul. Warszawska 24, 31-155 Cracow, Poland

Abstract Numerical methods for solving nonlinear systems of weakly singular Volterra integral equations (VIEs) possessing weakly singular solutions appear almost nonexistent in the literature, except for a few treatments of single first kind Abel equations. To reduce this gap, an extension is presented, of the adaptive Huber method designed for VIEs with singular kernels such as K (t, τ ) = (t − τ )−1/2 and K (t, τ ) = exp[−λ(t − τ )](t − τ )−1/2 (where λ ≥ 0) and a variety of nonsingular kernels. The method was thus far restricted to bounded solutions having at least two derivatives. Under a number of assumptions specified, the extension applies to solutions Uµ (t) that can be written as sums of singular components cµ t−1/2 (with unknown coefficients cµ ), and nonsingular components U µ (t). In the solution process, factor t−1/2 is handled analytically, whereas cµ and U µ (t) are determined numerically. Computational experiments reveal that the extended method determines singular solutions equally well as the unextended method determined nonsingular solutions. The method is intended primarily for a class of VIEs encountered in electroanalytical chemistry, but it can also be of interest to other application areas. Keywords: Volterra integral equations, weakly singular kernels, weakly singular solutions, adaptive methods, product-integration, computational electrochemistry

1. Introduction The present study is inspired by the class of one-dimensional Volterra integral equations (VIEs) of the first and second kind, encountered in electroanalytical chemistry, and representing models of transient electrochemical experiments [1]. VIEs of this class typically involve integrals with weakly singular convolution ∗ Corresponding

author, tel. +48(12)6282670, http://www.cyf-kr.edu.pl/˜ nbbienia Email address: [email protected] (L.K. Bieniasz)

Preprint submitted to Journal of Computational and Applied Mathematics

April 12, 2017

2 kernels such as K (t, τ ) = (t − τ )−1/2 and K (t, τ ) = exp[−λ(t − τ )](t − τ )−1/2 with λ ≥ 0, but also simultaneously integrals with a variety of nonsingular kernels. Both single VIEs and VIE systems are of interest, and the equations can be either linear or nonlinear. In the latter case the integrands are usually linearly dependent on the unknown functions, but the resulting integrals and unknown functions form nonlinear expressions. Hence, the VIEs considered can generally be written as systems: F (t, U (t), Y (t)) = 0

(1)

T

where F (·) = [F1 (·), . . . , FM (·)] is an M −dimensional vector of functions repT resenting the individual VIEs, 0 is the zero vector, U (t) = [U1 (t), . . . , UM (t)] is an M −dimensional vector of unknown functions of independent variable t, T and Y (t) = [Y1 (t), . . . , YI (t)] is an I−dimensional vector of integrals: ˆt Yi (t) =

Kκ (t, τ ) Uµ (τ ) dτ,

(2)

0

where Kκ (t, τ ) with κ = κ(i) = 1, . . . , L is one of L possible kernels, and Uµ (t) with µ = µ(i) = 1, . . . , M represents one of the unknowns U (t) (the one associated with the ith integral). In order to solve the above VIEs numerically, an adaptive strategy based on the concept of product-integration, and on the classical discretisation of Huber [2], has been elaborated by the present author (see Refs. [3–7] for the description of the strategy, and Refs. [8–21] for its applications to a variety of electrochemical models characterised by specific kernel functions). It was demonstrated that the strategy operates in an almost automatic way, allowing one to obtain solutions with an accuracy close to the preselected error tolerance. However, the strategy has been thus far limited to VIEs possessing bounded, nonsingular solutions. In particular, the existence of at least two solution derivatives of U (t) was assumed for t belonging to the interval [0, tmax ] of interest (although the adaptive algorithm was found to operate satisfactorily also when the derivatives did not exist at t = 0). The limitation to such regular solutions is a deficiency of the strategy, since in a part of electroanalytical models one expects weak singularities of (some or all of) the solutions at t → 0: Uµ (t) ∼ t−1/2 , which can be proven in a quite general manner (see, for example, Refs. [22, 23]). If present, a solution singularity can sometimes be neglected or eliminated by reformulating the VIEs (albeit at the expense of an additional human effort). This is why accounting for the presence of the singularities was not the first priority in Refs. [3–21]. But since there are also situations when the singular components dominate in the solutions, it would be desirable to extend the strategy from Refs. [3–21], so that it can determine singular solutions numerically without any need for “manual” VIE reformulation. Such an extension is described in the present paper. We concentrate here on a general presentation of the extended strategy, and on its tests using non-electrochemical example VIEs, for which analytical solutions exist. Hence, the operation of the strategy

3

can be examined quantitatively. In addition, the results obtained may also be of interest to a wider audience of modellers from VIE application areas beyond electrochemistry. Solutions of somewhat more complicated example VIEs corresponding to electrochemical models are planned to be described separately in an electrochemical journal. Although numerical methods for solving weakly singular VIEs have been intensively investigated over the past decades (see, for example, Refs. [24–29] and references cited therein), the literature contains surprisingly few records of efforts towards numerically determining singular solutions of such VIEs. The records seem to have been limited to several proposals for the numerical treatment of single first kind Abel integral equations [30], for which the existence of singular terms in the solutions is well known and proven. Such proposals are generally hard-to-find, as the singularity of the solutions may not even be explicitely mentioned in the texts. In particular, Lubich [31] described fractional linear multistep methods theoretically applicable in such cases (although an example solved by him did not have a singular solution). Tausch [32] presented a generalized Euler-Maclaurin formula for the numerical solution of Abel-type equations, and applied it to an example having a singular solution. In their discussion of several iterative techniques, Pandey et al. [33] gave an example of deriving a singular solution, but it is debatable whether such techniques should be regarded as numerical, or they are rather analytical. Popular textbooks [24– 29] on the numerical methods for integral equations are also enigmatic about how to compute singular solutions. One exception is the book by Kythe and Puri [28], who briefly suggested that appropriate transformations of the unknowns could be used to eliminate singularities, prior to applying numerical methods suitable for nonsingular solutions (cf. pp. 175-176 in Ref. [28]). 2. The Extended Adaptive Strategy In this section we describe the essentials of the extended adaptive strategy. For generally formulated VIEs (1) it is difficult to find and/or obtain any complete mathematical treatments of the existence, uniqueness, stability, etc., of the solutions. It would also be a formidable task to derive a mathematical theory of all aspects of the present method. Therefore, an experimental route was taken, by continuing the approach adopted in previous works [3–21]. As a result, a practical algorithm has been constructed, which is likely to work, provided that a number of relatively straightforward assumptions is satisfied. The assumptions are formulated in subsect. 2.1. New details of the extended strategy are given in subsect. 2.2. Results of computational experiments, confirming the validity of the method, are provided in Sect. 3. Finally, conclusions are summarised in Sect. 4. 2.1. Assumptions The adaptive strategy from Refs. [3–21] relies on replacing the unknown solutions U (t) by piecewise linear spline functions of t, using a dynamically

2.1

Assumptions

4

selected grid of nodes tn (n = 0, 1, . . . , N ) along the t axis, with local step sizes hn = tn − tn−1 . After analytical integration of the splines one obtains systems of nonlinear algebraic equations for the approximate values of the unknowns at every successive node tn . The nonlinear algebraic systems are solved by the multi-dimensional Newton iterative method. Hence, it is assumed that Eq. (1) either has a unique solution U (t), or, if the solution is not unique, then the particular solution of interest can be reached by a proper selection of starting values for Newton iterations. This obviously imposes some requirements on the regularity of U (t), as well as F (·), which however are non-trivial to formulate theoretically, as was noted. The existence and continuity of the first partial derivatives of F (·), with respect to Uµ (t) and Yi (t) is at least desired by the Newton method. Apart from determining approximate solutions, local error estimates are calculated and used for deciding whether a discrete solution obtained is acceptable, or has to be recalculated using a modified step size hn . For all these purposes, an availability of analytical expressions, or highly accurate approximants (with errors at the level of the errors of the machine representation of real numbers) for the moment integrals: ˆt Kκ (t, τ ) τ m dτ

(3)

0

(with m = 0, 1, 2) is required for all kernels of interest and for t ≥ 0. Error estimates are additionally dependent on (proportional to) second solution derivatives. Solutions Uµ (t) that are linear functions of t are determined exactly (with machine accuracy) by the method. To incorporate the possibility of handling singular solutions, we make the following additional assumptions: (A) Unknown functions Uµ (t) (for µ = 1, . . . , M ) can be formally represented as sums of singular and nonsingular components, that is: Uµ (t) = cµ t−1/2 + U µ (t),

(4)

where cµ denotes unknown coefficients, and U µ (t) satisfies the regularity assumptions sufficient for the correct operation of the unextended atrategy from Refs. [3–21]. Obviously, some sort of mathematical insights (into the VIEs to be solved) is needed to verify that this assumption is plausible; in the electrochemical modelling this is a typical situation, as was noted. The integrals of the nonsingular components will be denoted by Y i (t) (for i = 1, . . . , I): ˆt Y i (t) =

Kκ (t, τ ) U µ (τ ) dτ, 0

(5)

2.2

New details of the strategy

5

(B) Integrals of the singular terms: Zeκ (t) =

ˆt Kκ (t, τ ) τ −1/2 dτ,

(6)

0

(for κ = 1, . . . , L) must be computable either analytically, or highly accurate approximants of quality comparable with that of integrals (3) must be available for t ≥ 0. Hence, apart from integrals (3), one more exact integral (6) is needed (per kernel function), by the present extension of the strategy, which is a relatively modest complication. Table 1 provides analytical expressions for integrals Zeκ (t) corresponding to the few simplest kernels that are of interest to electrochemical modelling. (C) In cases when the VIEs are of second kind, we require that the dependence of F (·) on U (t) vanishes in the limit t → 0. This means that lim ∂Fj /∂Uµ = 0 t→0

(with j, µ = 1, . . . , M ). Furthermore, we make this assumption stronger by requiring that lim Uµ ∂Fj /∂Uµ = 0, in order to make sure that there remains t→0

no unbounded contribution from U (t) in F (·), at t = 0. We thus consider only such VIEs that can be represented directly as finite expressions at t = 0, without applying additional variable transformations. For the same reason we exclude the possibility of other unbounded terms in F (·). Hence, for example, the VIE considered by Pandey et al. [33]: ˆt 0

(t − τ )−1/2 U (τ ) dτ + U (t) − π − t−1/2 = 0,

(7)

with the solution U (t) = t−1/2 , cannot be solved by the present method, because in this case ∂F/∂U |t=0 = 1 6= 0, and also the last term is singular.

(D) The system (1) must be solvable at t = 0, when considered as an algebraic equation system for the unknown coefficients cµ . This assumption should become clearer in the next subsection. 2.2. New details of the strategy Under assumptions A-D, Eq. (1) can be rewritten as   e (c, t) + U (t), Ye (c, t) + Y (t) = 0, F t, U

(8)

 T  T e (c, t) = where U (t) = U 1 (t), . . . , U M (t) , Y (t) = Y 1 (t), . . . , Y I (t) , U h iT T  −1/2 eκ(1) (t), . . . , cµ(I) Zeκ(I) (t) , with c1 t , . . . , cM t−1/2 , and Ye (c, t) = cµ(1) Z µ(i) and κ(i) representing indices of the unknowns and kernels associated with T any ith integral. The unknowns to be determined are c = [c1 , . . . , cM ] and U (t).

2.2

New details of the strategy

6

By analogy with Refs. [3–21] we seek to determine approximate values un = T [u1,n , . . . , uM,n ] of U (tn ) for n = 0, 1, . . . , N , and estimates est(δn ) of the T true errors δn = [δ1,n , . . . , δM,n ] of un , defined as δn = U (tn ) − un .

(9)

Vector c can, in turn, be determined as a part of the initialisation procedure associated with the calculation of u0 and u1 . Specifically, for t = t0 = 0 Eq. (8) becomes   F 0, ∗, Ye (c, 0) = 0, (10)

where (in view of assumption C) the star indicates the lack of the dependence on U (t) in the limit of t → 0, and the term Y (0) is omitted, since for nonsingular solution components U (t) Eq. (5) predicts that Y (0) = 0. Equation (10) is a system of M nonlinear algebraic equations, with (the only) M unknowns c, since appropriate exact expressions for Zeκ (0) are available (assumption B). Assumption D guarantees the solvability of Eq. (10). Note that by applying the multidimensional Newton method, Eq. (10) can be solved almost exactly (practically with machine accuracy), if a sufficient number of iterations is performed. In the case of linear VIEs Eq. (10) is also linear, so that already one iteration gives c with machine accuracy. However, although Eq. (10) suggests that the calculation of c can be decoupled from the calculation of of u0 and u1 , we argue that it is more elegant to combine these calculations. In this way one obtains a single nonlinear algebraic equation system to be solved for the unknowns at the first integration step, instead of two such systems (one for c, and one for u0 and u1 ). Consequently, the initialisation procedure from Ref. [6] is modified in the following way. Let t1/2 = t1 /2 = h1 /2, u1/2 = u1/2 (u0 , u1 ) =

u0 + u1 , 2

u0 + u1 , 2 y 1 = y 1 (u0 , u1 ) = R1,0,1 u0 + S1,0,1 u1 ,

y 1/2 = y 1/2 (u0 , u1 ) = R1/2,0,1/2 u0 + S1/2,0,1/2

(11) (12) (13)

where R1/2,0,1/2 , S1/2,0,1/2 , R1,0,1 and S1,0,1 are analytically derived matrices of kernel-dependent quadrature coefficients, defined in Ref. [6]. Using the above expressions, the nonlinear algebraic system to-be-solved at the first integration T step for the block vector [c, u0 , u1 ] of 3M unknowns, can be written as:     F 0, ∗, Ye (c, 0) = 0      e (c, t1/2 ) + u1/2 (u0 , u1 ) , Ye (c, t1/2 ) + y 1/2 (u0 , u1 ) = 0 . F t1/2 , U      e (c, t1 ) + u1 , Ye (c, t1 ) + y (u0 , u1 ) = 0  F t1 , U 1 (14) Equation (14) supersedes Eq. (13) in Ref. [6]. Due to the presence of the additional unknowns c, pre-scaling [6] of the linearised Eq. (14) by certain powers of h1 is not feasible, and was not used here.

7 e (c, t) and Ye (c, t) become known (almost) Once c is determined, both U exactly, so that Eq. (8) can be viewed as an analytically transformed VIE system (1), in which the vector of unknowns is U (t), rather than the original (total) solution U (t). Formally, the vector function F (·) can then be replaced by a new vector function F (·), such that  
e (c, t) + U (t), Ye (c, t) + Y (t) , F t, U (t), Y (t) = F t, U (15) so that the transformed VIE becomes


F t, U (t), Y (t) = 0.

(16)

Therefore, the numerical determination of u2 , u3 , . . . from Eq. (16) is analogous to the algorithm previously described in Ref. [6] for nonsingular solutions, and needs not be re-described here. The same refers to the calculation of est(δn ) for n = 0, 1, . . . , N , because all these errors are evaluated when c is already determined. In our computer implementation of the extended strategy the substitution (15) is performed internally by the code, so that no “manual” reformulation of the VIEs is needed. Having determined any of the successive vectors u0 , u1 , . . ., and their local error estimates est(δn ), the step sizes hn are adjusted automatically so that the estimated errors do not exceed a prescribed tolerance, in a way analogous to that described in Refs. [3–6]. An important property of the extended strategy is that since the coefficients c are determined almost exactly, the errors of the iT h −1/2 −1/2 are negligible compared singular solution components c1 tn , . . . , cM tn to the errors of un . Therefore, the error estimates calculated by the strategy for un can be simultaneously regarded as error estimates of the approximate total iT h −1/2 −1/2 T + un . solutions, equal to un = [u1,n , . . . , uM,n ] = c1 tn , . . . , cM tn 3. Computational experiments The extended adaptive strategy has been tested using 14 examples of single VIEs and their systems, and examination of further (electrochemical) examples is in progress, with similar results. Tables 2-4 contain 9 most interesting examples selected for the discussion in this paper. The interval of t was taken to be [0, 1] in all examples. Solutions reported below correspond to the following selection of the parameter values: λ = 1.2 (in example 2); λ = 1.3 (in example 3); λ1 = 1.5 and λ2 = 1.3 (in example 6); a = 0.4, b = 0.6, λ1 = 0.7 and λ2 = 0.8 (in example 7); k = 100 (in example 8); k = 10 (in example 9). The numerical algorithm was coded in C++, using the Bloodshed/Orwell Dev − C + + 5.7.1 environment [34, 35], and compiled using the TDM GCC 4.8.1 compiler (a 32 bit release). The resulting program was run as a sequential, 32-bit console application on a laptop computer with an Intel Centrino 2 processor (2.4 GHz), operating under Windows Vista. All calculations were done

8

in extended precision (18-19 significant digits). Special functions needed for the calculations were computed by using the present author translations (into C++) of the relevant routines from the publicly available (from netlib [36]) packages: SPECFUN of Cody [37] (functions erf, erfc, erex, daw, I0 , I1 , besei0 , besei1 , J0 , J1 ), and MISCFUN of Macleod [38] (functions H0 , H1 ). The relative errors of the special functions calculated in this way were checked not to exceed about 10−16 , and they were most often even lower. The adaptive strategy is generally found to perform very satisfactorily, when applied to the present class of example VIEs with singular solutions. Observations regarding the accuracy, convergence, reliability of the error estimation, automatism of the grid adaptation, numerical stability and efficiency, are consistent with the previous findings concerning the operation of the method for nonsingular solutions [3–21]. Figures 1 and 2 depict representative numerical solutions and their estimated errors, in comparison with exact solutions and true errors. As can be seen, numerical solutions match well exact solutions, and estimated errors match well true errors. The errors are also well controlled: their magnitude is close to the requested error tolerance parameter tol. Similar results were obtained for all examples. It is important to point out that in the presence of the initial singularities the adaptive algorithm adequately selects the first integration step h1 . Even if the starting step value hstart (chosen on input) is too large to guarantee the requested accuracy, the step is automatically decreased, in a loop of gradual reductions, until the initial errors become sufficiently small. Such successful initial step adjustments were observed in all examples. The step adaptation mechanism works correctly despite the potential numerical problems that one might expect at t ≈ 0, as a result of largely disparate magnitudes of the singular and nonsingular solution components in Eq. (4). No such or other difficulties were observed for tol & 10−6 . In some of the examples, smaller tolerances were sporadically found to be too demanding for the present method, which manifested itself by the lack of error reductions upon step reductions (but not necessarily at the first integration step), apparently as a result of an interference of machine errors enhanced by small step sizes. The quality of the adaptive error control is further illustrated by the plots showing the proportionality of the maximum solution error max |error|, to the error tolerance tol. Figure 3 presents such plots, separately for the first kind VIEs, and second kind VIEs. The data for the plots were obtained by beginning calculations with the largest tol = 10−3 , and repeating them for gradually diminished tol, until tol = 10−8 was reached, or until the aforementioned failures of the error reduction were observed. Hence, for some examples there are no data points in Fig. 3, corresponding to the lowest tol values. As can be seen, the equality max |error| ≈ tol is well obeyed in the case of first kind VIEs. A slight departure from this equality occurs for second kind VIEs, but max |error| is then still approximately a linear function of tol, so that by varying tol one is able to modify the actual error. Analogous results were previously obtained for nonsingular solutions in Refs. [3–21] (cf., in particular, discussions of Fig. 4 in Ref. [3], Fig. 4 in Ref. [5], and Fig. 3 in Ref. [6]). The true error decrease upon decreasing tol, seen in Fig. 3, confirms also the

9

convergence of the method, since with the decreasing tol the integration steps are reduced, while the number N of accepted time steps increases. From the almost linear dependence of the decimal logaritm of max |error| on the decimal logarithm of N , one can determine the practical (or average) convergence order of the method. The order is equal to the absolute value of the first derivative of the dependence. Figure 4 demonstrates that for the examples studied the convergence order is very close to 2 (for first kind VIEs it is practically equal 2; for second kind VIEs it seems to be a bit closer to 1.8). This result is again in accord with previous findings for nonsingular solutions [3–21] (cf., in particular, discussions of Fig. 5 in Ref. [3], Fig. 5 in Ref. [5], and Fig. 4 in Ref. [6]). The present work reveals no new, unexpected findings, concerning the numerical stability of the adaptive method, compared to what is already known. The stability requirements are difficult to formulate rigorously for general Eq. (1), but at least two heuristic requirements have been conjectured based on extensive experimentation [6, 7]. First, in a linearised Eq. (1), a “total” kernel of every VIE, resulting from summing up kernels of individual integrals involving a given unknown, should be a non-increasing function of t − τ . Second, in the case of second-kind VIEs that can be rearranged to the form U (t) − G (t, Y (t)) = 0

(17)

with M = I, where every unknown Uµ (t) has a corresponding integral (5) involving the same unknown, all eigenvalues of G(·) should have negative real parts. Efficiency plots, given in Fig. 5, reveal that the computing time ct needed for solving examples 1-9 varies between about 10−2 second, and about 103 seconds, when tol varies between 10−3 and 10−6 − 10−8 , depending on the examples and equation parameters. This is somewhat more than was previously observed for nonsingular solutions [3–21], but the previous calculations were on a faster computer. The present timings may also not be directly comparable with the previous ones, because of a number of modifications introduced into the adaptive code, referring to the data structures used and C++ class hierarchies. The present example 6 is solved the most efficiently, and example 8 is solved the least efficiently (partly due to the expensive calculations of several special functions involved in example 8, but mostly because of the persistent variations of U (t) over the entire interval of t).

4. Conclusions Results obtained in this study demonstrate that it is possible to extend the adaptive Huber method from Refs. [3–21], formerly restricted to weakly singular VIEs (1) possessing nonsingular solutions, onto a wider class of weakly singular VIEs having singular solutions. The extension is valid under assumptions A-D listed in Sect. 2.1. The extended method possesses practically identical characteristics (error magnitudes, convergence order, automatism of the t−grid

10

adaptation, proportionality of the errors to the error tolerance parameter, stability and efficiency), as it had before the extension. This considerable enlargement of the application scope of the method is achieved at the moderate extra price of obtaining analytical expressions for integrals (6) (or developing highly accurate approximants for these integrals). References [1] L. K. Bieniasz, Modelling Electroanalytical Experiments by the Integral Equation Method, Springer, Berlin, 2015. [2] A. Huber, Eine Näherungsmethode zur Auflösung Volterrascher Integralgleichungen, Monatsschr. Math. Phys. 47 (1939) 240–246. [3] L. K. Bieniasz, An adaptive Huber method with local error control, for the numerical solution of the first kind Abel integral equations, Computing 83 (2008) 25–39. [4] L. K. Bieniasz, Initialisation of the adaptive Huber method for solving the first kind Abel integral equation, Computing 83 (2008) 163–174. [5] L. K. Bieniasz, An adaptive Huber method for weakly singular second kind Volterra integral equations with non-linear dependencies between unknowns and their integrals, Computing 87 (2010) 35–54. [6] L. K. Bieniasz, An adaptive Huber method for non-linear systems of weakly singular second kind Volterra integral equations, Appl. Math. Comput. 217 (2011) 5622–5631. [7] L. K. Bieniasz, Extension of the adaptive Huber method for Volterra integral equations arising chemistry, to convolution kernels  in electroanalytical  exp[−α(t − τ )]erex [β(t − τ )]1/2 and exp[−α(t − τ )]daw [β(t − τ )]1/2 , J. Comput. Meth. Sci. Eng. 11 (2011) 323–338. [8] L. K. Bieniasz, Cyclic voltammetric current functions determined with a prescribed accuracy by the adaptive Huber method for Abel integral equations, Anal. Chem. 80 (2008) 9659–9665. [9] L. K. Bieniasz, Automatic simulation of cyclic voltammograms by the adaptive Huber method for weakly singular second kind Volterra integral equations, Electrochim. Acta 55 (2010) 721–728. [10] L. K. Bieniasz, Automatic simulation of cyclic voltammograms by the adaptive Huber method for systems of weakly singular Volterra integral equations, J. Electroanal. Chem. 642 (2010) 127–134. [11] L. K. Bieniasz, Extension of the adaptive Huber method for solving integral equations occurring in electroanalysis, onto kernel function representing fractional diffusion, Electroanalysis 23 (2011) 1506–1511.

11

[12] L. K. Bieniasz, A highly accurate, inexpensive procedure for computing integral transformation kernel and its moment integrals for cylindrical wire electrodes, J. Electroanal. Chem. 661 (2011) 280–286. [13] L. K. Bieniasz, Automatic simulation of electrochemical transients at cylindrical wire electrodes, by the adaptive Huber method for Volterra integral equations, J. Electroanal. Chem. 662 (2011) 371–378. [14] L. K. Bieniasz, Automatic simulation of electrochemical transients by the adaptive Huber method for Volterra integral equations involv ing kernel terms exp[−α(t − τ )]erex [β(t − τ )]1/2 and exp[−α(t −  τ )]daw [β(t − τ )]1/2 , J. Math. Chem. 50 (2012) 765–781.

[15] L. K. Bieniasz, Automatic solution of integral equations pertinent to diffusion with first order homogeneous reactions at cylindrical wire electrodes, J. Electroanal. Chem. 674 (2012) 38–47. [16] L. K. Bieniasz, Automatic simulation of electrochemical transients, assuming finite diffusion space at planar interfaces, by the adaptive Huber method for Volterra integral equations, J. Electroanal. Chem. 684 (2012) 20–31. [17] L. K. Bieniasz, Automatic solution of the Singh and Dutt integral equations for channel or tubular electrodes, by the adaptive Huber method, J. Electroanal. Chem. 693 (2013) 95–104. [18] L. K. Bieniasz, Automatic solution of integral equations describing electrochemical transients under conditions of internal spherical diffusion, J. Electroanal. Chem. 694 (2013) 104–113. [19] L. K. Bieniasz, Automatic solution of integral equations describing electrochemical transients under conditions of internal cylindrical diffusion, J. Electroanal. Chem. 700 (2013) 30–39. [20] L. K. Bieniasz, Automatic solution of integral equations describing electrochemical transients at dropping mercury electrodes, J. Electroanal. Chem. 705 (2013) 44–51. [21] L. K. Bieniasz, A new theory, and automatic computation of reversible cyclic voltammograms at a microband electrode, J. Electroanal. Chem. 767 (2016) 123–133. [22] C. G. Phillips, K. M. Jansons, The short-time transient of diffusion outside a conducting body, Proc. R. Soc. Lond. A 428 (1990) 431–449. [23] K. B. Oldham, The short-time chronoamperometric behaviour of an electrode of arbitrary shape, J. Electroanal. Chem. 297 (1991) 317–348. [24] C. T. H. Baker, The Numerical Treatment of Integral Equations, Clarendon Press, Oxford, 1978.

12

[25] P. Linz, Analytical and Numerical Methods for Volterra Equations, SIAM, Philadelphia, 1985. [26] H. Brunner, P. J. van der Houwen, The Numerical Solution of Volterra Equations, North-Holland, Amsterdam, 1986. [27] W. Hackbusch, Integral Equations, Theory and Numerical Treatment, Birkhäuser, Basel, 1995. [28] P. K. Kythe, P. Puri, Computational Methods for Linear Integral Equations, Birkhäuser, Boston, 2002. [29] H. Brunner, Collocation Methods for Volterra Integral and Related Functional Differential Equations, Cambridge University Press, Cambridge, 2004. [30] R. Gorenflo, S. Vessella, Abel integral equations, analysis and applications, Lect. Notes Math. 1461 (1991) 1–215. [31] C. Lubich, Fractional linear multistep methods for Abel-Volterra integral equations of the first kind, IMA J. Numer. Anal. 7 (1987) 97–106. [32] J. Tausch, The generalized Euler-Maclaurin formula for the numerical solution of Abel-type integral equations, J. Integral Eqs. Appl. 22 (2010) 115–140. [33] R. K. Pandey, O. P. Singh, V. K. Singh, Efficient algorithms to solve singular integral equations of Abel type, Comput. Math. Appl. 57 (2009) 664– 676. [34] Bloodshed Software, Dev-C++, http://www.bloodshed.net, Accessed 3 February 2017 (2017). [35] Orwell, Dev-C++ Blog, http://orwelldevcpp.blogspot.com, Accessed 3 February 2017 (2017). [36] AT&T Bell Laboratories, The University of Tennessee, Oak Ridge National Laboratory, Netlib Repository at UTK and ORNL, http://www.netlib.org, Accessed 3 February 2017 (2017). [37] W. J. Cody, Algorithm 715: SPECFUN - a portable FORTRAN package of special function routines and test drivers, ACM Trans. Math. Softw. 19 (1993) 22–32. [38] A. J. Macleod, Algorithm 757: MISCFUN, a software package to compute uncommon special functions, ACM Trans. Math. Softw. 22 (1996) 288–301. [39] H. E. Fettis, On the numerical solution of equations of the Abel type, Math. Comput. 18 (1964) 491–496.

}

} (

π 1/2 (λt)1/2 2

π [1 − exp(−λt)] 2λ1/2 h (λt)2 πλ1/2 t 1 − λt 2 2 + 6

2 t1/2 1 −




i

2048 (λt)11 316234143225

λt 2

(λt)14 1307674368000

+ ... +

+ besei1

2 λt 3




− ... −

+

   ( πλ )1/2h 1-erex (λt)1/2

λt 2

for λ > 0 for λ = 0

for λt < 0.4

for λt < 0.16

for λt ≥ 0.16

for λt ≥ 0.4

i

eκ (t) for several kernels Kκ (t, τ ). Symbols: erfc, daw, and erex, denote respectively: the complementary Table 1: Analytically obtainable integrals Z error function, the Dawson integral, and the expression erex(x) = exp(x2 )erfc(x). Symbols besei0 and besei1 represent expressions: besei0 (x) = exp(−x) I0 (x) and besei1 (x) = exp(−x) I1 (x), where I0 and I1 denote modified Bessel functions of the first kind and orders 0 and 1. Parameter λ is eκ (t) are ill-conditioned for small λt, so that assumed to be real and nonnegative in all cases. In the case of the last two kernels the expressions for Z well conditioned approximations based on power series expansions are also provided

1/2

daw{[λ(t − τ )]

erex{[λ(t − τ )]

1/2

(

 2 t1/2 − (πλ)1/2 t besei0




erfc{[λ(t − τ )]1/2 }

λt 2

π besei0

(t − τ )−1/2 exp[−λ(t − τ )]

2λ−1/2 daw[(λt)1/2 ] 2 t1/2

π



2 t1/2

eκ (t) Z

(t − τ )−1/2

exp[−λ(t − τ )]

1

Kκ (t, τ )

13

14

No 1

2 3

4

5

Integral Equations ´t U (τ ) dτ − π exp(t) = 0 (t−τ )1/2 0

U (τ ) (t−τ )1/2

dτ − 2λ1/2 0

´t

exp[−λ(t−τ )] U (τ ) dτ (t−τ )1/2



dτ + π2 t3/2 U (t) − f (t) = 0

t2 2

 dτ + exp [−U (t)] − f (t) = 0

π 2

+ 2t1/2

daw{[λ(t − τ )]1/2 }U (τ ) dτ − π = 0

(considered by Fettis [39]) ´t 0

´t

0

(1 + 2λt)Y1 (t) + 3(πλ)1/2 Y2 (t) − f (t) = 0 where Y1 (t) =

U (τ ) (t−τ )1/2

´t Y2 (t) = erfc{[λ(t − τ )]1/2 }U (τ ) dτ 0
π 5
1/2 f (t) = πbesei 0 λt ] 2 + 3 2 − λt erf[(λt) +(πλt)1/2 6 + (πλt)1/2 + 4λt − exp(−λt)

´t 0

U (τ ) (t−τ )1/2

where  f (t) = π 1 + t +  t ´ exp −

0

where 
  f (t) = exp −π 1 + t1/2 + t + exp − t−1/2 +




Solutions

+ 2t1/2

+ 2λt1/2

U (t) = t−1/2 + π 1/2 exp(t)erf(t1/2 )

(πλ)1/2 3

U (t) = t−1/2 + 2λt1/2 U (t) = t−1/2 +

π 2

U (t) = t−1/2 + t1/2

U (t) = t−1/2 +

Table 2: Examples of single VIEs used for numerical experiments. Symbol λ denotes a real, non-negative parameter

0

exp[−λ2 (t−τ )] U2 (τ ) dτ (t−τ )1/2

0

0 ´t

exp[−λ2 (t−τ )] U2 (τ ) dτ, (t−τ )1/2

exp[−λ1 (t − τ )] U1 (τ ) dτ

0




λ2 t 2
λ2 t 2

a11 (t) = 1 + t1/2 , a12 (t) = 1 + t, a23 (t) = 1 + t3/2 , a24 (t) = 1 + t2 b11 (t) = t4 , b12 (t) = t3 , b21 (t) = t2 , b22 (t) = t f1 (t), f2 (t) − are appropriately selected to ensure given solutions

Y4 (t) =

Y3 (t) =

´t

0

0
λ2 
λ1 t 2t f1 (t) = λ21 besei − 2 besei0 λ
 + besei1 1 2 2 
f2 (t) = λ21 besei0 λ21 t + besei1 λ21 t − λ22 besei1  a11 (t)Y1 (t) + a12 (t)Y2 (t)    +b11 (t)U1 (t) + b12 (t)U2 (t) − f1 (t) = 0 a (t)Y  23 3 (t) + a24 (t)Y4 (t)   +b21 (t)U1 (t) + b22 (t)U2 (t) − f2 (t) = 0 where ´t ´t U1 (τ ) dτ, Y (t) = U2 (τ ) dτ Y1 (t) = (t−τ 2 1/2 )

Y3 (t) =

0

´t

Integral Equations  1 λ1 Y (t) − 2π Y2 (t) − f1 (t) = 0 π 1/2 1 λ2 1 Y (t) + 1 1/2 2π Y3 (t) − f2 (t) = 0 π where ´t exp[−λ1 (t−τ )] ´t U1 (τ ) dτ, Y (t) = U2 (τ ) dτ Y1 (t) = (t−τ 2 1/2 ) (t−τ )1/2



U1 (t) = t−1/2 + t1/2 exp(−at) U2 (t) = − exp(−bt) + 1 + t1/2 t1/2

Solutions ( 2 t)−exp(−λ1 t) U1 (t) = exp(−λ2(πt 3 )1/2 −1/2 U2 (t) = t

Table 3: Examples of linear VIE systems used for numerical experiments. Symbols λ1 , λ2 , a, and b, denote real, non-negative parameters. Expression for U1 (t) in example 6 is ill-conditioned for small t and must be evaluated using power series expansions for exponential functions.

7

6

No

15

16

No 8

9

´t

0

U2 (τ ) (t−τ )1/2

dτ

− f1 (t) = 0 i2 h Y2 (t) − f2 (t) = 0 k1/2

dτ, Y2 (t) =

Integral Equations  h i h ih i  Y1 (t) 3 (t) Y2 (t)  + Yπ11/2 π 1/2 k1/2 ´t

0

U1 (τ ) (t−τ )1/2

 a (t) Y (t)2 + a (t) Y (t) −  1 1 2 1

where Y1 (t) =

a1 (t) = π4 − 4H1 (ϑ) + πH1 (ϑ)2 a2 (t) = 2π 1/2 J1 (ϑ)H0 (ϑ) [2 − πH1 (ϑ)] f1 (t) = J0 (ϑ)3 + 2J0 (ϑ)2 + πJ0 (ϑ) [J1 (ϑ)H0 (ϑ) − J0 (ϑ)H1 (ϑ)] 2 f2 (t) = − [πH0 (ϑ)J1 (ϑ)] , with ϑ = 2(kt)1/2

( Y (t) 1 3/2 + (πt) U 1 (t)U2 (t) − f1 (t) = 0 π 1/2 Y2 (t) 5/2 2 3/2 2 − (πt) U U 1 (t) + (πt) 2 (t) − 2πt U2 (t) − f2 (t) = 0 π 1/2 where i h 1/2 f1 (t) = I0 (ϑ) + (πt)1/2 cosh(ϑ) + (πt)2 sinh(2ϑ) f2 (t) = 1 + (πt)1/2 [I1 (ϑ) − 1 − πt] , with ϑ = 2(kt)1/2

U2 (t) =

1

Solutions   U (t) = 

(

U1 (t) = U2 (t) =

cos[2(kt)1/2 ] (πt)1/2 sin[2(kt)1/2 ] πt

cosh[2(kt)1/2 ] (πt)1/2   1 + sinh 2(kt)1/2 (πt)1/2

Table 4: Examples of nonlinear VIE systems used for numerical experiments. Symbols J0 and J1 denote Bessel functions of the first kind and orders 0 and 1, symbols H0 and H1 denote Struve H-functions. Symbol k denotes a real, non-negative parameter

17

a

6

U1(t), U2(t)

4 2 0 -2 -4 -6 0

0.2

0.4

0.6

0.8

1

t 1

1

b

0.5 error × 103

error × 103

0.5

c

0

-0.5

0

-0.5

-1

-1 0

0.2

0.4

0.6 t

0.8

1

0

0.2

0.4

0.6

0.8

t

Figure 1: Adaptive numerical solutions u1,n and u2,n obtained for k = 100 in exam-

ple 8 (a), and corresponding errors of u1,n (b) and u2,n (c). Notations in subfigure a: approximate solutions u1,n - (black circles); approximate solutions u2,n - (white circles); exact solutions U1 (t) and U2 (t) - (solid lines). Notations in subfigures b and c: estimated errors est(δ1,n ) or est(δ2,n ) - (white circles joined by dotted lines); true errors δ1,n or δ2,n - (black dots joined by solid lines). Parameters of the adaptive Huber method are: starting step hstart = 0.01, maximum step hmax = 0.1, error tolerance tol = 10−3 . The adapted first integration step is h1 ≈ 3.72 × 10−11

1

18

300

a

250

U1(t), U2(t)

200

150

100

50

0 0

0.2

0.4

0.6

0.8

1

t 0.2

1.5

b

0

c

1

error × 103

error × 103

-0.2 -0.4 -0.6 -0.8

0.5

0

-0.5

-1

-1 0

0.2

0.4

0.6 t

0.8

1

0

0.2

0.4

0.6

0.8

t

Figure 2: Adaptive numerical solutions u1,n and u2,n obtained for k = 10 in example 9 (a), and corresponding errors of u1,n (b) and u2,n (c). Notations and method parameters are as in Fig. 1. The adapted first integration step is h1 ≈ 4.38 × 10−8

1

19

-3

a

log( max|error|)

-4

-5

-6

-7

-8 -8

-7

-6

-5

-4

-3

-4

-3

log(tol)

-3

b

log( max|error|)

-4

-5

-6

-7

-8 -8

-7

-6

-5

log(tol)

Figure 3: Dependence of the decimal logarithm of the maximum true error max |error| (with error standing for δn , δ1,n or δ2,n ), on the decimal logarithm of the error tolerance parameter tol, observed in the case of first kind (a), and second kind (b) example VIEs. Notations in subfigure a: circles - example 1; squares - example 2; triangles with vertices up - example 3; triangles with vertices down - example 6; diamonds example 8. Notations in subfigure b: circles - example 4; squares - example 5; triangles with vertices up - example 7; triangles with vertices down - example 9. Black symbols refer to the errors of un or u1,n . White symbols refer to the errors of u2,n . All symbols are joined by straight solid lines. The dashed line represents the equality: max |error| = tol. The values of tol for the different points in the plots are: 10−3 , 5 × 10−4 , 2 × 10−4 , 10−4 , 5 × 10−5 , 2 × 10−5 , 10−5 , 5 × 10−6 , 2 × 10−6 , 10−6 , 5 × 10−7 , 2 × 10−7 , 10−7 , 5 × 10−8 , 2 × 10−8 , and 10−8 (from right to left). True errors of u2,n in example 6 are not seen among the plots, since they are much smaller than 10−8 (U2 (t) is determined with machine accuracy, because the nonsingular solution component equals zero, which is a special case of a linear function)

20

a

-3

log( max|error|)

-4

-5

-6

-7

-8 1

1.5

2

2.5

3

3.5

4

4.5

3.5

4

4.5

log(N)

b

-3

log( max|error|)

-4

-5

-6

-7

-8 1

1.5

2

2.5

3

log(N)

Figure 4: Dependence of the decimal logarithm of the maximum true error max |error|, on the decimal logarithm of the number N of integration steps, observed in the case of first kind (a), and second kind (b) example VIEs. Dashed lines are plots of the linear functions: log (max |error|) = −1.5 − 2 log(N ) and log (max |error|) = 2.1 − 2 log(N ) (in subfigure a), and log (max |error|) = −0.7 − 1.8 log(N ) and log (max |error|) = −1.8 log(N ) (in subfigure b). Other notations are as in Fig. 3

21

-3

a

log( max|error|)

-4

-5

-6

-7

-8 -2

-1

0

1

2

3

1

2

log(ct /s)

-3

b

log( max|error|)

-4

-5

-6

-7

-8 -3

-2

-1

0

log(ct /s)

Figure 5: Dependence of the decimal logarithm of the maximum true error max |error|, on the decimal logarithm of the computational time ct (in seconds), observed in the case of first kind (a), and second kind (b) example VIEs. Notations are as in Fig. 3. Computational times obtained for tol ∈ [2 × 10−5 , 10−3 ] are averaged over 10-20 program runs