Next Article in Journal
Fuzzy Logic Prediction of Hypertensive Disorders in Pregnancy Using the Takagi–Sugeno and C-Means Algorithms
Next Article in Special Issue
Comprehensive Numerical Analysis of Time-Fractional Reaction–Diffusion Models with Applications to Chemical and Biological Phenomena
Previous Article in Journal
A Self-Training-Based System for Die Defect Classification
Previous Article in Special Issue
Inequalities for Riemann–Liouville-Type Fractional Derivatives of Convex Lyapunov Functions and Applications to Stability Theory
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Solution of Linear Second-Kind Convolution Volterra Integral Equations Using the First-Order Recursive Filters Method

Laboratoire Angevin de Mécanique, Procédés et innovAtion (LAMPA), Arts et Métiers ParisTech, 49035 Angers, France
Mathematics 2024, 12(15), 2416; https://doi.org/10.3390/math12152416
Submission received: 28 June 2024 / Revised: 25 July 2024 / Accepted: 1 August 2024 / Published: 3 August 2024

Abstract

:
A new numerical method for solving Volterra linear convolution integral equations (CVIEs) of the second kind is presented in this work. This new approach uses first-order infinite impulse response digital filters method (IIRFM). Three convolutive kernels were analyzed, the unit kernel and two singular kernels: the logarithmic and generalized Abel kernels. The IIRFM is based on the combined use of the Laplace transformation, a first-order decomposition, and a bilinear transformation. This approach often leads to simple analytical expressions of the approximate solutions, enabling efficient numerical calculation, even using single-precision floating-point numbers. When compared with the method of homotopic perturbations with Laplace transformation (HPM-L), the IIRFM approach does not present, in linear cases, the convergence difficulties inherent to iterative approaches. Unlike most solution methods based on the Laplace transform, the IIRFM has the dual advantage of not requiring the calculation of the Laplace transform of the source function, and of not requiring the systematic calculation of inverse Laplace transforms.

1. Introduction

In the present work, a new method is presented for the numerical solution of the convolutive Volterra integral equations (CVIEs) of the second kind, written in the following form:
u λ , f ( t ) = f ( t ) + λ 0 t K ( t , t ) a L u λ , f ( t ) + a NL NL u λ , f ( t ) d t ,
where K ( t , t ) is the kernel of the integral equation; f ( t ) is a source function assumed to be continuous; λ , a L , a NL R are parameters; NL [ u λ , f ( t ) ] is a non-linear operator of the desired solution u λ , f ( t ) . These integral equations are used in a wide variety of scientific areas, such as point mechanics (Abel’s problem), fluid mechanics (water waves scattering, [1]), electrochemistry (voltammetry, [2]), as well as heat transfer (see Section 4).
The linear situation, which is the only one studied in this first paper, corresponds to a NL = 0 , in which case we find Volterra’s linear integral equations of the second type, defined, for  t 0 and a L = 1 [3], by:
u λ , f ( t ) = f ( t ) + λ 0 t K ( t , t ) u λ , f ( t ) d t .
The present paper studies a particularly important class of Volterra integral equations of the second type, which involves convolutional kernels, depending only upon the difference t t , and written in the general form K ( t , t ) = K ( t t ) . Three kinds of convolution kernels are studied in the present work: the unit kernel K ( t , t ) = 1 ; generalized Abel’s kernel K ( t , t ) = 1 / ( t t ) α with 0 < α < 1 and logarithmic kernel K ( t , t ) = ln t t . The last two convolution kernels share the property of being singular at t = t .
Note that using the generalized Abel kernel in Volterra’s integral Equation (2) leads to the generalized Abel integral equation of the second kind, as follows:
u λ , α , f ( t ) = f ( t ) + λ 0 t u λ , α , f ( t ) t t α d t .
A general analytical solution exists in closed form for Abel’s integral equation, involving the calculation (explicit or numerical) of an integral and a series of functions (see [4] (p. 142)):
u λ , α , f ( t ) = f ( t ) + 0 t R ( t t ) f ( t ) d t , with R ( t ) = n = 1 λ Γ ( 1 α ) t 1 α n t Γ n ( 1 α ) ,
where Γ ( z ) = 0 e x x z 1 d x is the Euler gamma function. Result (4) will provide analytical reference solutions for some of the numerical results that will be calculated in this work.
Several approaches based on the Laplace transformation have already been proposed to solve analytically or numerically Volterra (CVIE) or Abel (CAIE) second-kind linear convolution integral equations. In [5,6], the Homotopic Perturbation Method (HPM) and the Laplace transformation were used to numerically determine the solution of various convolution Abel integral equations of the second kind. In [7], systems of coupled convolutive Abel integral equations have been solved numerically using the homotopic analysis method (HAM) introduced by Lio [8,9], coupled with the Laplace transformation. In [10], Bairwa et al. used the q-homotopic analysis and Laplace transformation method (q-HATM) to numerically solve various CAIEs.
All these approaches generally assume that there is an analytical expression for the source function f ( t ) . There are, however, situations where it is not possible to obtain this expression; for example, when f ( t ) is a noisy signal or the result of an experimental measurement. In such cases, the IIR filters method presented here can be very useful for determining the numerical solution of the linear convolution Volterra equations.
The main objective of the present study is therefore to apply the first-order infinite impulse response filters method (IIRFM) proposed by Heyd et al. [11,12] to numerically solve linear Volterra convolutive integral equations of the second kind. The use of infinite impulse response (IIR) digital filters is usually limited to linear time-invariant systems (LTISs), whose dynamics are described by linear ordinary differential equations with constant coefficients, or equivalently by Laplace transfer functions expressed as polynomial fractions. The results presented in this paper make it possible, for the first time, to extend the field of application of IIR filters to linear systems described by Volterra convolutive integral equations of the second kind, which are involved in a wide range of scientific fields, including heat transfer, for example.
This new approach has the great advantage of being non-iterative, which means it has none of the convergence problems inherent in iterative methods. In addition, for some Volterra integral equations, the IIRFM can lead to an analytical expression of the approximate numerical solution, making the calculation very efficient compared with iterative numerical methods.
The results obtained with the IIRFM are compared in this paper with those provided by other numerical approaches, including the homotopic perturbation method with Laplace transformation (HPM-L), the Daubechies wavelet method of Mouley et al. [13], and the polynomial approximation method of Singha et al. [14].
Several new results have been obtained in this work using the IIRFM, such as a very efficient numerical calculation of the function g λ ( t ) = 1 + erf ( λ π t ) when λ < 0 (see Section 3.2.1) and this even using single-precision floating-point numbers. This new approach to calculating the erf ( λ π t ) function could be very usefully incorporated into the free codes currently available, Python, Julia, Octave, and Scilab.
The remainder of this paper is organized as follows. The foundations of the IIRFM and HPM-L methods are outlined in Section 2. The mechanism of the IIRFM is illustrated in Section 3 by numerically solving various linear CVIEs of the second kind. The results obtained by the IIRFM are also compared with those provided by different numerical approaches. Finally, two examples of linear CVIEs used in heat transfers are solved in Section 4 using the IIRFM method.

2. IIRFM and HPM-L Numerical Methods

The two numerical approaches used in this work to solve the integral Equation (2) are outlined in this paragraph, and the main steps involved in their implementation are recalled or introduced here.

2.1. First-Order IIR Filters Method (IIRFM)

The IIRFM consists of numerically solving Equation (2) in discrete time, using sufficiently precise recurrence relations, which are deduced from successive applications of the Laplace transformation, a first-order partial fractions decomposition and of a bilinear transformation.
This method is particularly well suited to solving convolutional linear integral equations (CVIEs), but also ordinary differential equations (ODEs), integro-differential equations (IDEs), and parabolic partial differential equations (PDEs). Each of these differential or integral problems is assumed to be linear and to have constant coefficients.

2.1.1. Principle of the IIRFM Method

The IIRFM consists schematically of replacing each of the differential or integral problems mentioned above, by a set of M independent, first-order, exact, or approximate ordinary differential problems.
Let ( Σ ) be a differential (or integral), linear and invariant problem of the previously mentioned kinds. Let f ( t ) be the excitation (or source function) and u λ , f ( t ) the response (or solution) of the problem ( Σ ) to the excitation f ( t ) , for a set of constant λ coefficients (see Figure 1).
Let us symbolically denote F ( u λ , f , u λ , f , u ˙ λ , f , f , λ , t ) = 0 the equation that governs the dynamics of ( Σ ) in direct space (see Figure 1a). To simplify the presentation, we assume at this stage that the initial condition on u λ , f ( t ) is null: u λ , f ( 0 ) = 0 . Since the functional relationship F = 0 is linear and has constant coefficients, the Laplace transformation L can be applied to it. This gives the algebraic relationship u ¯ λ , f ( s ) = H ( s ) f ¯ ( s ) , which governs the dynamics of ( Σ ) in the Laplace space, where H ( s ) is the Laplace transfer function of the system (see Figure 1b), with  u ¯ λ , f ( s ) = L [ u λ , t ( t ) ] and f ¯ ( s ) = L [ f ( t ) ] .
To transform the equation F ( u λ , f , u λ , f , u ˙ λ , f , f , λ , t ) = 0 into M independent, linear, first-order ODEs with constant coefficients, the function H ( s ) C 0 , with  C 0 as a constant, is first approximated by a rational fraction RF M ( s ) of the following form:
RF M ( s ) = n = 0 M 1 h n ( λ ) s n m = 0 M m ( λ ) s m ,
where { h n , m } is a set of 2 M + 1 constant real coefficients, determined by interpolation (by least squares or Chebyshev for example) or approximation (of Padé for example) of the H ( s ) C 0 function. Note that the value of the constant C 0 depends on the problem under consideration. In the case of linear convolutional Volterra integral equations of the second kind (CVIE), which will be examined in detail later in this work, the value C 0 = 1 is chosen.
The rational fraction RF M ( s ) is then decomposed, by partial fractions decomposition, in the usual form, as follows:
RF M ( s ) = i = 1 N 1 r i s p i + j = 1 N 2 c j + d j s e j + f j s + g j s 2 , = i = 1 N 1 H i 1 ( s ) + j = 1 N 2 H j 2 ( s ) ,
where N 1 + 2 N 2 = M and c j , d j , e j , f j , g j , r i , p i are real constants. Note that it is assumed here that the multiplicity ν i of each pole p i is unity. This property is usually guaranteed by the fact that the poles are approximated and are calculated numerically here. Finally, each second-order transfer function H j 2 is decomposed into single complex elements, in the following form:
H j 2 ( s ) = α j + i β j s z ̲ j + α j i β j s z ̲ j * ,
where z ̲ j = x j + i y j and z ̲ j * = x j i y j are the complex conjugate roots (with y j 0 ) of the polynomial equation e j + f j s + g j s 2 = 0 . The real constants α j and β j are α j = d j / 2 and β j = 1 2 ( c j + d j x j ) / y j , respectively. The previous results can be combined by writing RF M ( s ) as a combination of M first-order transfer functions with complex coefficients, as follows:
RF M ( s ) = i = 1 M a M , i ( λ ) b M , i ( λ ) + s ,
where a M , i ( λ ) , b M , i ( λ ) Z . Finally, the exact H ( s ) transfer function of the system is replaced by a transfer function H M ( s ) written in the following form:
H M ( s ) = C 0 + i = 1 M a M , i ( λ ) b M , i ( λ ) + s .
This is a founding result of the IIRFM. Note that in certain situations, the transfer function H M ( s ) may coincide with the exact transfer function (see for example Section 3.1).
Using result (9), the approximate solution u ¯ λ , f IIR , M ( s ) can be expressed in the Laplace space as
u ¯ λ , f IIR , M ( s ) = H M ( s ) f ¯ ( s ) = C 0 f ¯ ( s ) + i = 1 M a M , i ( λ ) b M , i ( λ ) + s f ¯ ( s ) ,
which can also be written as u ¯ λ , f IIR , M ( s ) = u ¯ λ , f , 0 IIR , M ( s ) + i = 1 M u ¯ λ , f , i > 0 IIR , M ( s ) , with
u ¯ λ , f , 0 IIR , M ( s ) = C 0 f ¯ ( s ) and u ¯ λ , f , i > 0 IIR , M ( s ) = a M , i ( λ ) b M , i ( λ ) + s f ¯ ( s ) ,
which are the approximate partial solutions of the problem in the Laplace space.
The approximate solution u λ , f IIR , M ( t ) of the problem in the direct space remains to be determined. Several approaches are possible (without being exhaustive):
  • It is possible to consider that each first-order partial transfer function a M , i ( λ ) / [ b M , i ( λ ) + s ] corresponds to the Laplace transform of a first-order linear ODE with constant coefficients and zero initial condition, for  i = 1 , 2 , , M :
    a M , i ( λ ) b M , i ( λ ) + s = u ¯ λ , f , i > 0 IIR , M ( s ) f ¯ ( s ) u ˙ λ , f , i > 0 IIR , M ( t ) + b M , i ( λ ) u λ , f , i > 0 IIR , M ( t ) = a M , i ( λ ) f ( t ) ,
    where u λ , f , i > 0 IIR , M ( 0 ) = 0 . The partial solutions u λ , f , i > 0 IIR , M ( t ) are obtained by applying the method of variation of constants, giving
    u λ , f IIR , M ( t ) = C 0 f ( t ) + i = 1 M a M , i ( λ ) 0 t f ( t ) e b M , i ( λ ) ( t t ) d t .
    Depending on the expression of f ( t ) , the M convolutional integrals figuring in (13) will be calculated analytically or numerically.
  • Noting that L 1 [ a M , i ( λ ) b M , i ( λ ) + s ] = a M , i ( λ ) e b M , i ( λ ) t , the result (13) can be found directly by applying the convolution theorem.
  • Finally, it is also possible to transform ordinary differential Equation (12) into partial recurrence relations by applying a bilinear transformation s 2 T e z 1 z + 1 [12] to each partial transfer function a M , i ( λ ) / [ b M , i ( λ ) + s ] . It follows that the partial solution u λ , f , i > 0 IIR , M ( t k ) verifies the following recurrence relation for t k = k T e , k 1 , i = 1 , , M and u λ , f , i > 0 IIR , M ( t 0 ) = 0 :
    u λ , f , 0 IIR , M ( t k ) = C 0 f ( t k ) and u λ , f , 0 IIR , M ( t 0 ) = 0 , u λ , f , i > 0 IIR , M ( t k ) = 2 T e b M , i ( λ ) 2 + T e b M , i ( λ ) u λ , f , i > 0 IIR , M ( t k 1 ) + a M , i ( λ ) 2 + T e b M , i ( λ ) f ( t k ) + f ( t k 1 ) ,
    where T e > 0 is the calculation step.
This latter approach will be discussed in detail in the remainder of this paper, which focuses on the numerical solution of linear convolutional Volterra integral equations of the second type (CVIE), for which it is recalled that C 0 = 1 .

2.1.2. Application of the IIRFM to the Case of Second-Kind CVIEs

The different steps of the IIRFM are schematically as follows:
Step 1.
The Laplace transformation ( L ) is applied to Equation (1):
L ( 1 ) u ¯ λ , f ( s ) = f ¯ ( s ) + λ L 0 t K ( t t ) a L u λ , f ( t ) + a NL NL u λ , f ( t ) d t ,
with u ¯ λ , f ( s ) = L [ u λ , f ( t ) ] and f ¯ ( s ) = L [ f ( t ) ] . Using the convolution theorem, the Laplace transform of the previous integral can be written as follows:
L 0 t K ( t t ) a L u λ , f ( t ) + a NL NL u λ , f ( t ) d t = K ¯ ( s ) × a L u ¯ λ , f ( s ) + a NL NL ¯ u λ , f ( s ) ,
where K ¯ ( s ) = L [ K ( t , 0 ) ] and NL ¯ [ u λ , f ] ( s ) = L { NL [ u λ , f ( t ) ] } , of which the analytical expression is usually unknown.
This gives the following expression for u ¯ λ , f ( s ) :
u ¯ λ , f ( s ) = 1 1 λ a L K ¯ ( s ) f ¯ ( s ) + λ a NL K ¯ ( s ) 1 λ a L K ¯ ( s ) NL ¯ [ u λ , f ] ( s ) .
Let us write H L ( s , λ ) = 1 / ( 1 λ a L K ¯ ( s ) ) and H NL ( s , λ ) = λ a NL K ¯ ( s ) / ( 1 λ a L K ¯ ( s ) ) . Three of the following situations are possible:
(a)
a NL = 0 and a L = 1 : this is the linear situation, for which the Volterra problem of the second kind can be described as a linear dynamics system with Laplace transfer function H L ( s , λ ) = u ¯ λ , f ( s ) / f ¯ ( s ) = 1 / ( 1 λ K ¯ ( s ) ) .
With the exception of the special case of the unit kernel K ( t , t ) = 1 , note that the transfer function H L ( s , λ ) does not have the usual form of a quotient of polynomials in s, as it is usually the case for LTIS. Several examples of linear situations are presented in Section 3.
(b)
a NL = 1 and a L = 0 : the result is a purely non-linear situation. To solve the problem numerically, the expression of the Laplace transform NL ¯ [ u λ , f ] ( s ) must be specified. Various approaches can be considered for this purpose, such as the Adomian polynomials method [15,16].
(c)
a NL 0 and a L 0 : this is a hybrid situation, in which the IIRFM numerical solution uses both the Laplace transfer function and the Adomian polynomials approach.
We now explain the next steps of the IIRFM, but only in the case of the linear situation ( a NL = 0 and a L = 1 ).
Step 2.
To solve the linear Volterra problem using the transfer function H L ( s , λ ) , it is first rewritten as a linear combination of M + 1 first-order transfer functions (first-order partial fractions decomposition, noted FOD):
H L M ( s , λ ) = 1 + i = 1 M a M , i ( λ ) b M , i ( λ ) + s ,
with a M , i , b M , i Z . This particular expression (other expressions are possible) of the transfer function H L ( s , λ ) is the founding point of the IIRFM applied to second-kind CVIEs. By introducing M + 1 partial transfer functions H L , i M ( s , λ ) , we can rewrite Equation (18) as H L M ( s , λ ) = i = 0 M H L , i M ( s , λ ) with H L , 0 M ( s , λ ) = 1 . With these notations, u ¯ λ , f ( s ) can still be written as:
u ¯ λ , f ( s ) = i = 0 M H L , i M ( s , λ ) f ¯ ( s ) = i = 0 M u ¯ λ , f , i ( s ) .
Depending on the convolution kernel considered, Equation (18) can be an exact form (see Section 3.1) or an approximation of H L ( s , λ ) (see Section 3.2, Section 3.3 and Section 4). In the latter case, coefficients a M , i and b M , i are determined by rational interpolation of the function G ( s , λ ) = H L ( s , λ ) 1 . We will come back to this point later.
The choice of (18) ensures that the approximate form H L M ( s , λ ) tends to be 1 for s , as does the exact form H L ( s , λ ) for the convolution kernels considered in this study, unit, generalized Abel, and logarithmic. This also reduces the interpolation interval to a set of s values close to 0, as will be shown later.
The solution obtained by the IIRFM is now denoted by u λ , f IIR , M ( t ) .
Step 3.
A bilinear transformation s 2 T e z 1 z + 1 ([12]) is applied to the transfer function H L M ( s , λ ) , where T e R + is a uniform time step (or sampling period) and z is the complex number introduced in the definition of the z-transform (noted here Z ) of a discrete-time signal x ( k T e ) : Z x ( k T e ) = x ¯ ( z ) = k = 0 x ( k T e ) z k , with  k N . After a few algebraic manipulations, we obtain the following z transfer function:
H L M ( z , λ ) = 1 + i = 1 M A M , i + A M , i z 1 1 B M , i z 1 = i = 0 M H L , i M ( z , λ ) ,
with H L , 0 M ( z , λ ) = 1 , H L , i > 0 M ( z , λ ) = ( A M , i + A M , i z 1 ) / ( 1 B M , i z 1 ) , A M , i = a M , i T e / ( 2 + b M , i T e ) , B M , i = ( 2 b M , i T e ) / ( 2 + b M , i T e ) and b M , i T e 2 .
The z-transform u ¯ λ , f IIR , M ( z ) obtained by the IIRFM satisfies the following relations:
u ¯ λ , f IIR , M ( z ) = H L M ( z , λ ) f ¯ ( z ) = i = 0 M H L , i M ( z , λ ) f ¯ ( z ) ,
where f ¯ ( z ) = Z f ( k T e ) is the z-transform of the source function f ( t ) . The solution u ¯ λ , f IIR , M ( z ) can be further decomposed into a sum of M + 1 partial solutions, as follows: u ¯ λ , f IIR , M ( z ) = i = 0 M u ¯ λ , f , i IIR , M ( z ) with u ¯ λ , f , i IIR , M ( z ) = H L , i M ( z , λ ) f ¯ ( z ) .
Step 4.
The solution is calculated in the discrete time space using recurrence relations deduced from transfer function (20). In order not to unnecessarily complicate the theoretical developments in this paragraph, it is assumed here that the function f ( t ) is zero at the initial instant, which also imposes a zero initial condition for the solution: u λ , f ( 0 ) = 0 . The case of a non-zero initial condition will be considered later in this work (see Section “Exponential Source Term f ( t ) = e t ”). The discrete-time solution at time t k = k T e is denoted by u λ , f IIR , M ( k T e ) = i = 0 M u λ , f , i IIR , M ( k T e ) .
Using time-shifting properties Z x ( k m ) T e = z m Z x k T e and linearity Z a m x k T e = a m Z x k T e of the z-transform, the following M + 1 partial recurrence relations are obtained:
u λ , f , 0 IIR , M ( k T e ) = f ( k T e ) k > 0 and u λ , f , 0 IIR , M ( 0 ) = 0 , u λ , f , i 0 IIR , M ( k T e ) = B M , i u λ , f , i 0 IIR , M ( ( k 1 ) T e ) + A M , i f ( k T e ) + f ( ( k 1 ) T e ) , k 1 , with u λ , f , i 0 IIR , M ( 0 ) = 0 .
Note that the second recurrence relation of (22) leads to the following explicit expression for the partial solution u λ , f , i 0 IIR , M ( k T e ) :
u λ , f , i 0 IIR , M ( k T e ) = A M , i j = 0 k 1 B M , i k ( j + 1 ) f ( j T e ) + f ( ( j + 1 ) T e ) .
This is a particularly important result of this study. Depending on the expression of the source term f ( t ) , the expression of Equation (23) can be further simplified, as will be seen with some of the examples to be studied in Section 3.
The complete solution u λ , f IIR , M ( k T e ) provided at time t k = k T e by the IIRFM is therefore of the following general form:
u λ , f IIR , M k T e = f k T e + i = 1 M A M , i j = 0 k 1 B M , i k ( j + 1 ) f ( j T e ) + f ( ( j + 1 ) T e ) .
In the most general case, the calculation of u λ , f IIR , M k T e involves two sums, with a total of M × k terms. As a result, calculation times can become significant as k increases. However, these calculation times can be significantly reduced when an explicit expression Ξ M , i k , T e , B M , i , f of j = 0 k 1 B M , i k ( j + 1 ) f ( j T e ) + f ( ( j + 1 ) T e ) exists. In this case, the numerical solution u λ , f IIR , M k T e can be written as
u λ , f IIR , M k T e = f k T e + i = 1 M A M , i Ξ M , i k , T e , B M , i , f ,
and then the calculation only involves a single sum of M terms, which significantly reduces computation times.
Equations (24) and (25) are the key results of the IIRFM applied to the numerical solution of linear CVIEs of the second kind. In contrast to the method of homotopic perturbations with Laplace transformation (HPM-L), which will be presented in the following paragraph, it is not necessary with the IIRFM to calculate the Laplace transform f ¯ ( s ) of the source function f ( t ) in order to numerically solve the problem under consideration. This is a considerable advantage of the IIRFM, as there are situations where the calculation of f ¯ ( s ) is not possible. This is the case, for example, if  f ( t ) includes a random component (like noise) or is the result of experimental measurements.
The block diagram shown in Figure 2 summarizes the main steps of the IIRFM as applied to the numerical solution of linear CVIEs of the second kind.

2.2. Homotopic Peturbation Method with Laplace Transform (HPM-L)

In what follows, some of the results obtained by the IIRFM will be compared, among others, with the results provided by the homotopic perturbation method with Laplace transformation (HPM-L), an approach that also relies on an initial application of the Laplace transformation to Equation (2). We therefore recall the main results that define this method in the linear case only ( a NL = 0 and a L = 1 ).
The homotopic perturbation method (HPM) was introduced by He [17,18,19] to solve ordinary differential equations (ODEs) or partial differential equations (PDEs), whether linear or non-linear. This method can also be applied to the case of linear CVIEs of type (2) and to ordinary integro-differential equations. In the present case, the HPM approach consists first of constructing a homotopy H ( v , p ) associated with Equation (2), where p [ 0 , 1 ] is an embedded parameter and v is an approximate solution to Equation (2). The homotopy is constructed in such a way that at the limit where p 1 , the solution v of the homotopic equation H ( v , p ) = 0 converges to the solution u of the problem under consideration. The method of homotopic perturbations with Laplace transform (HPM-L) used in this work consists of using the Laplace integral transformation L in the homotopy construction process. It is possible to construct H ( v , p ) directly from the Laplace transform of a differential Equation [20], which is generally well suited to the case of non-steady PDEs. In the case of EDOs and CVIEs, it is generally preferable to rewrite the equations using the direct L and inverse L 1 Laplace transformations [21,22]:
L L u λ , f ( t ) = L f ( t ) + λ L K ( t , 0 ) × L u λ , f ( t ) .
By noting L K ( t , 0 ) = K ¯ ( s ) , the Laplace transform of (2) can be written in the following form:
L u λ , f ( t ) = L f ( t ) + λ K ¯ ( s ) × L u λ , f ( t ) .
It remains to apply the inverse Laplace transformation L 1 to algebraic Equation (27), which leads to the following statement for u λ , f ( t ) :
u λ , f ( t ) = f ( t ) + λ L 1 K ¯ ( s ) × L u λ , f ( t ) .
Finally, from (28) the following H ( v , p ) homotopy can be constructed:
H ( v , p ) = v ( t ) f ( t ) p λ L 1 K ¯ ( s ) × L v ( t ) ,
which verifies the following relations:
H ( v , 0 ) = 0 v ( t ) = f ( t ) , H ( v , 1 ) = 0 v ( t ) = f ( t ) + λ L 1 K ¯ ( s ) × L v ( t ) .
The second equation in (30) shows that the solution to the homotopic equation H ( v , 1 ) = 0 is indeed the solution to the initial problem (2). To solve this homotopic equation, the perturbation method is generally used, which consists of writing the solution v H ( t , p ) of (29) as a series expansion of the embedded parameter p, as follows:
v H ( t , p ) = n = 0 p n v n ( t ) .
Assuming that series (31) is convergent, we then have the following result:
v ( t ) = lim p 1 v H ( t , p ) = v 0 ( t ) + v 1 ( t ) + + v n ( t ) + .
Two situations can be encountered:
  • lim p 1 v H ( t , p ) converges to a closed form, in which case the homotopic perturbation method leads to an analytic solution of Equation (2);
  • There is no known closed form of lim p 1 v H ( t , p ) , in which case the series is generally truncated at the N H + 1 first terms and the solution of the problem is written in the following approximate form:
    v N H ( t ) = n = 0 N H v n ( t ) .
It is the last situation that is of particular interest in the present work, because even if series (32) is convergent, the value of the parameter λ can considerably influence the speed of convergence of v N H ( t ) . Situations may therefore arise where the number N H + 1 of terms to be considered may become numerically prohibitive, making the HPM-L method imprecise, inoperative, or very difficult to implement.
The v n ( t ) functions are obtained by incorporating Equation (31) into the homotopic equation H ( v , p ) = 0 . Using the linearity of the Laplace transformation, the following recurrence relation of the given problem is obtained as follows:
n = 0 p n v n ( t ) = f ( t ) + λ n = 0 p n + 1 L 1 K ¯ ( s ) × L v n ( t ) .
By identifying the different powers of p, we deduce the equations verified by the following v n ( t ) functions:
n = 0 : v 0 ( t ) = f ( t ) , n = 1 : v 1 ( t ) = λ L 1 K ¯ ( s ) × L v 0 ( t ) , n = 2 : v 2 ( t ) = λ L 1 K ¯ ( s ) × L v 1 ( t ) , n = N : v N ( t ) = λ L 1 K ¯ ( s ) × L v N 1 ( t ) ,
In the following section, the IIRFM and HPM-L methods are used to numerically solve various linear CVIEs of the second kind.

3. Linear Convolutive Volterra Integral Equations of the Second Kind

Several linear examples ( a L = 1 and a NL = 0 ) that allow one to compare the relative accuracies and efficiencies of the HPM-L and IIRFM numerical methods (see Section 3.1) are presented in this paragraph. The IIRFM is also compared with other recent numerical approaches (see Section 3.2.3 and Section 3.2.4).

3.1. Unit Kernel

To obtain an analytical solution of the linear Volterra integral equation with unit kernel K ( t , t ) = 1 , we derive Equation (3) with respect to t, and then integrate the resulting ordinary equation using the variation of constants method. After integration by parts, we obtain the following expression for the solution u λ , f ( t ) :
u λ , f ( t ) = f ( t ) + λ 0 t e λ t t f ( t ) d t .
This equation is of the form of Equation (4), with  R ( t ) = λ e λ t .

3.1.1. Monomial Source Term f ( t ) = t n

We first consider a monomial source term f ( t ) = t n , which leads to a relatively simple analytical expression for the solution, denoted here by u λ , t n ( t ) :
u λ , t n ( t ) = t n 1 + e λ t ( λ t ) n Γ ( n + 1 ) Γ ( n + 1 , λ t ) ,
where Γ ( n , z ) = z x n 1 e x d x is the incomplete gamma function and Γ ( n + 1 ) = Γ ( n + 1 , 0 ) . For  t R + and n N , Equation (37) admits the following Maclaurin series expansion:
u λ , t n , m ( t ) = t n + λ t n + 1 n + 1 + λ 2 t n + 2 n + 1 n + 2 + + λ m t m + n n + 1 n + m + O ( t m + n + 1 ) .
Note that in the special case where n = 1 , the analytical solution u λ , t ( t ) is given by the following expression, deduced from (37):
u λ , t ( t ) = e λ t 1 λ ,
whose limit for t tends to 1 / | λ | when λ < 0 .

Numerical Solution by the IIRFM

The application of the IIRFM presented in Section 2.1, leads in the present case to the following relations, with  K ¯ ( s ) = 1 / s and H L ( s , λ ) = s / ( s λ ) :
L u λ , f ( t ) = H L ( s , λ ) L f ( t ) L u λ , f ( t ) = s s λ L f ( t ) .
For the case f ( t ) = t n , one finds L f ( t ) = Γ n + 1 s ( n + 1 ) , in which case it is possible to analytically invert the Laplace transform (40), giving back the analytical result (37).
Rather than an analytical solution, we are looking here for a numerical solution, designated by u λ , f IIR , deduced from relation (40). To this end, the bilinear transformation s 2 T e z 1 z + 1 is applied to the Laplace transfer function H L ( s , λ ) = s / ( s λ ) , without having to calculate the Laplace transform of the source function f ( t ) . The application of the bilinear transformation leads to the following z-transfer function:
H L ( z , λ ) = u ¯ λ , f IIR ( z ) f ¯ ( z ) = 1 z 1 1 λ T e 2 1 + λ T e 2 z 1 .
Using time-shifting properties TZ x ( k m ) T e = z m TZ x k T e and linearity TZ a m x k T e = a m TZ x k T e of the z-transform, we obtain the recurrence relation satisfied by the solution u λ , f IIR ( t k ) at the following discrete time t k = k T e :
u λ , f IIR ( t k ) = 1 + λ T e 2 1 λ T e 2 u λ , f IIR ( t k 1 ) + f ( t k ) f ( t k 1 ) 1 λ T e 2 ,
where u λ , f IIR ( t 0 ) = 0 and T e R + is the sample period, satisfying T e 2 / | λ | .
In the specific case f ( t ) = t n considered in this paragraph, f ( t k ) f ( t k 1 ) can be assimilated to T e f ( t k ) = n T e n k n 1 . It is then possible to derive an approximate explicit expression of u λ , t n IIR ( t k ) , which will be denoted by u ˜ λ , t n IIR ( t k )
u ˜ λ , t n IIR ( t k ) = n T e n w k Li 1 n w 1 λ T e 2 Φ w , 1 n , k + 1 1 + λ T e 2 ,
where w = ( 2 λ T e ) / ( 2 + λ T e ) ; Li m w = p = 1 w p / p m is the polylogarithm function and Φ w , a , b = p = 0 w p / ( p + b ) a is the Lerch transcendent function. These two functions are analytically extended for | w | > 1 .

Exact Solution by Homotopic Perturbation Method with Laplace Transformation (HPM-L)

The application of the HPM-L method presented in Section 2.2 leads in the present case to the following relationships, deduced from results (35) and from the Laplace transform L t n = Γ ( n + 1 ) / s n + 1 :   
0 : v 0 ( t ) = t n , 1 : v 1 ( t ) = λ L 1 L t n s 1 = λ Γ ( n + 1 ) L 1 s n 2 = λ Γ ( n + 1 ) Γ ( n + 2 ) t n + 1 , v 1 ( t ) = λ t n + 1 n + 1 , 2 : v 2 ( t ) = λ L 1 L λ t n + 1 n + 1 s 1 = λ 2 n + 1 Γ ( n + 2 ) L 1 s n 3 = λ 2 n + 1 Γ ( n + 2 ) Γ ( n + 3 ) t n + 2 , v 2 ( t ) = λ 2 t n + 2 n + 1 n + 2 , m : v m ( t ) = λ L 1 L λ m 1 t n + m 1 n + 1 n + 2 n + m 1 s 1 = λ m n + 1 n + 2 n + m 1 Γ ( n + m ) Γ ( n + m + 1 ) t n + m , v m ( t ) = λ m t n + m n + 1 n + m ,
The v ( t ) solution obtained by the HPM-L method is therefore written here, according to relation (32)
v ( t ) = t n + λ t n + 1 n + 1 + λ 2 t n + 2 n + 1 n + 2 + + λ m t n + m n + 1 n + m + , = n ! t n m = 0 λ t m n + m ! .
Since m = 1 λ t m n + m ! = e λ t ( λ t ) n Γ ( n + 1 ) Γ ( n + 1 , λ t ) / Γ ( n + 1 ) , we find analytical result (37).

Approximate Numerical Solutions for f ( t ) = t

The approximate numerical solutions obtained by the HPM-L and IIRFM methods, respectively, are now compared in the case where f ( t ) = t .
The HPM-L approximate numerical solution is obtained by truncating (44) at rank N H . Table 1 shows the absolute errors ε N H ( 5 | λ | ) = | u λ , t ( 5 | λ | ) v N H ( 5 | λ | ) | obtained when t = 5 | λ | , between the analytical value u λ , t ( 5 | λ | ) and the approximation v N H ( 5 | λ | ) , for different values of N H , when the parameter λ takes the values ± 0.5 , ± 1.0 and ± 1.5 .
For λ = 0.5 , 1.0 , and 1.5 , the results gathered in Table 1 show that the absolute error ε N H ( 5 | λ | ) is less than 0.01 for the respective truncation orders N H 6 , 16 , and 32. These results (as well as all the other results gathered in Table 1) show the influence of the value of λ on the speed of convergence of the HPM-L method. In particular, the larger the value of | λ | , the slower the convergence of the HPM-L method here. Figure 3a,b shows the evolution of the approximate solution v N H ( t ) as a function of t, for  n = 1 and, respectively, for λ = 1 and λ = 1.5 , for different values of N H . The numerical results (dashed curves) provided by the HPM-L method are compared with the analytical solutions (solid curves) u 1 , t ( t ) and u 1.5 , t ( t ) . These figures clearly show that the number N H of terms, required to obtain an approximate solution v N H ( t ) converging to the analytical solution u λ , 1 ( t ) , increases with the value of λ . From a numerical point of view, this can be a definite handicap for the HPM-L approach, especially when there is no closed form of the solution.
We are now looking for the numerical solution u ˜ λ , f IIR ( t k ) at t k = k T e , using the IIRFM. In the case of a unit kernel and for f ( t ) = t , result (43) takes on the following expression, valid whatever the sign of λ :
u ˜ λ , t IIR ( t k ) = 1 λ 2 + λ T e 2 λ T e k 1 λ .
This is a particularly simple expression, whose accuracy depends, for a given value of λ , on the choice of the value of T e . Note that lim t k u ˜ λ , t IIR ( t k ) 1 / | λ | when λ < 0 , which is consistent with the exact result of (39). Table 2 gathers the absolute errors ε T e ( 5 | λ | ) = | u λ , t ( 5 | λ | ) u ˜ λ , t IIR ( 5 | λ | ) | obtained when t k = 5 | λ | , between the analytical value u λ , t ( 5 | λ | ) and the approximate result u ˜ λ , t IIR ( 5 | λ | ) , for different values of T e , when λ takes the values ± 0.5 , ± 1.0 , and ± 1.5 .
The analysis of the results gathered in Table 2 shows that the IIRFM leads to highly accurate results for negative values of λ , even for sampling periods as large as T e = 0.1 . This is a very interesting result. As Equation (45) is simpler than Equation (39), which includes the incomplete gamma special function, it is also easier to implement numerically, and poses no particular convergence problems.
Figure 4 shows the variations of relative errors ε r NM ( t ) = | u NM ( t ) / u λ , t ( t ) 1 | as a function of t, where (NM) stands for “numerical method”, with NM≡HPM-L or IIRFM. The value of λ is negative here and equals λ = 1.5 . The calculation step is the same for both numerical approaches and is equal to δ t = T e = 10 2 . For  N H = 15 , we can clearly see from Figure 4 that the IIRFM performs significantly better than the HPM-L method as t deviates from zero.
Since v N H ( t ) has the same form as the Maclaurin series expansion (38), the approximate solution provided by the HPM-L method is very accurate in the vicinity of t = 0 , as shown in Figure 4a, even for small values of N H . In contrast, as soon as t deviates from the origin, the number N H of terms required to obtain a precise numerical value must increase accordingly, which can prove to be a problem from the point of view of numerical calculation (particularly in terms of calculation time). The numerical behavior of the IIRFM is very different from that of the HPM-L method. Indeed, as shown in Figure 4b, the accuracy of (45) increases with t k = k T e when λ < 0 .
Figure 5 shows the variations in the relative errors ε r HPM - L ( t ) and ε r IIRFM ( t ) as a function of t, when the parameter λ is positive and equal to λ = 1.5 , with all other quantities having the same values as those used for Figure 4. Even though the accuracy of the IIRFM decreases this time as t k deviates from the origin, it remains more accurate than the HPM-L method.
The order of accuracy n of the IIRFM approximation can be evaluated by computing the absolute value of the difference u ˜ λ , t IIR ( t k ) u λ , t ( t k ) , as a function of the sampling period T e . Using the expressions (39) and (45), we obtain
ε k ( T e ) = u ˜ λ , t IIR ( t k ) u λ , t ( t k ) = 1 | λ | e k λ T e 2 + λ T e 2 λ T e k .
Figure 6 shows the evolution of ε k ( T e ) as a function of T e [ 10 5 , 10 1 ] , when t k = 1 and λ = 1 . It can be seen from Figure 6 that the IIRFM approach is a second-order accuracy method in the present case.
In conclusion, it should be noted that, whatever the sign of λ , the accuracy of the IIRFM can be improved by reducing T e . This is not the case with the results obtained by the HPM-L method, whose accuracy is independent of T e .

3.1.2. Exponential Source Term f ( t ) = e t

As a preamble to this paragraph, it should be noted that it is not necessary for the initial condition on the solution to be null in order to apply the IIRFM in the case of the unit kernel.
The developments of the previous paragraph are repeated, but this time in the case of a source term f ( t ) of exponential form. In this case, there is also an analytical solution u λ , e t ( t ) of Equation (3), which is
u λ , e t ( t ) = e t + λ e λ t 1 + λ when λ 1 , ( 1 t ) e t when λ = 1 .
The most general case, where λ 1 , is considered in the following. For  t R + , Equation (47) admits a Maclaurin series expansion u λ , e t , m ( t ) of the following form to order m:
u λ , e t , m ( t ) = n = 0 m ( 1 ) m + λ m + 1 ( λ + 1 ) m ! t m + O ( t m + 1 ) .

Numerical Solution by the IIRFM

Since the transfer function H ( z , λ ) associated with Equation (3) is independent of the source function f ( t ) , the numerical solution u λ , e t IIR of the present problem still verifies recurrence relation (42). Equating the difference f ( t k ) f ( t k 1 ) with the function T e f ( t k ) = T e e k T e , it is also possible to find here an approximate explicit expression for the series u λ , e t IIR , T e ( t k ) , an expression which is denoted by u ˜ λ , e t IIR ( t k ) , as follows:
u ˜ λ , e t IIR , T e ( t k ) = e t k ω 1 + ω ω ω 1 t k / T e + e T e ω t k / T e e T e ω ,
where ω = ( 2 λ T e ) / ( 2 + λ T e ) , and ω = 2 T e / ( 2 + λ T e ) , with  T e 2 / λ . Note that here u ˜ λ , e t IIR , T e ( 0 ) = f ( 0 ) = 1 .

Solution by the HPM-L Method

The application of the HPM-L method leads to the following relationships:
0 : v 0 ( t ) = e t , 1 : v 1 ( t ) = λ L 1 L e t s 1 = λ L 1 1 s + s 2 , v 1 ( t ) = λ 1 e t , 2 : v 2 ( t ) = λ L 1 L λ 1 e t s 1 = λ 2 L 1 1 s 2 1 s + s 2 , v 2 ( t ) = λ 2 t 1 + e t , 3 : v 3 ( t ) = λ L 1 L λ 2 t 1 + e t s 1 = λ 3 L 1 1 s 3 1 s 2 + 1 s + s 2 , v 3 ( t ) = λ 3 1 t + 1 2 t 2 e t , m : v m ( t ) = λ L 1 L λ m 1 ( 1 ) m 2 k = 0 m 2 ( 1 ) k k ! t k + ( 1 ) m 1 e t s 1 , = λ m L 1 ( 1 ) m 2 k = 0 m 2 ( 1 ) k s k + 2 + ( 1 ) m 1 s + s 2 , v m ( t ) = λ m ( 1 ) m 1 k = 0 m 1 ( 1 ) k k ! t k + ( 1 ) m e t ,
The solution v ( t ) obtained by the HPM-L method is therefore written here, according to equations (32) and (50):
v ( t ) = m = 0 v m ( t ) = m = 0 λ m ( 1 ) m 1 k = 0 m 1 ( 1 ) k k ! t k + ( 1 ) m e t .
An approximate form of the solution, denoted by v N H ( t ) , is obtained by truncating series (51) at a rank N H < .

Numerical Results and Discussion

Figure 7a represents the numerical results obtained by the HPM-L (result (51)) and IIRFM (result (49)) methods, for  f ( t ) = e t , λ = 5 , T e = 10 2 , and t I t = [ 0 , 5 ] , using the semi-logarithmic representation. It can be seen that relatively high N H truncation orders have to be used for the numerical solution v N H ( t ) provided by the HPM-L method to converge to the exact solution (47) with acceptable accuracy over the entire interval I t . In contrast, it is observed that the numerical solution u ˜ λ , e t IIR , T e provides a very accurate approximation of the exact solution u λ ( t ) over the entire interval I t = [ 0 , 5 ] , even for a time step as large as T e = 10 2 .
The observations made in the previous discussion for λ = 5 are even more evident when the value λ = 5 is chosen. As shown in Figure 7b, convergence of the HPM-L method is slow in the case where λ = 5 . In the present case, no finite value of N H could be found that allows the series v N H ( t ) to converge to the exact solution u λ , e t ( t ) on the interval I t = [ 0 , 5 ] without using a convergence acceleration procedure. In contrast, we can observe that the IIRFM approach again gives a very accurate numerical approximation on I t , whatever the value of λ . It was also observed that the HPM-L method became unusable here from a numerical point when λ 10 .
In the present case, the IIRFM approach was found to be a first-order accuracy method.

3.2. Abel’s Kernel (with α = 1 / 2 )

In the case of Abel’s kernel, one obtains K ¯ ( s ) = Γ ( 1 α ) / s 1 α , and the transfer function H L ( s , λ ) is then given by
H L ( s , λ ) = 1 1 λ Γ ( 1 α ) s α 1 .
The results of Section 2 are applied here to the case α = 1 / 2 , for different source terms f ( t ) and different values of the parameter λ . Only the IIRFM approach is developed in detail in this section. To differentiate the results of this paragraph from those obtained in the previous paragraph for K ( t , t ) = 1 , we now agree to index all useful functions by the value of α .
In the present case, the Laplace transfer function H L , 1 2 ( s , λ ) of Abel’s problem is given, from Equation (17), by
H L , 1 2 ( s , λ ) = s s λ π .
In contrast to the case where K ( t , t ) = 1 studied in Section 3.1, there is no exact expression here for H L , 1 2 ( s , λ ) in the form of a quotient of polynomials in s. To apply the IIRFM in the present case, the function G 1 2 ( s , λ ) = H L , 1 2 ( s , λ ) 1 needs to be approximated by a rational fraction as accurate as possible. Several approaches make obtaining this approximation possible. In [23], an approximation of the transfer function was obtained using the least-squares method, while in [11,24] the rational fraction was constructed using Padé approximants. The second approach generally allows the construction of a very precise rational fraction, but mainly in the vicinity of a given point of the Laplace space. In order to obtain an accurate description over an extended domain of the Laplace space, and thus over an extended domain of the direct time space, a Chebyshev rational interpolation of the function G 1 2 ( s , λ ) , has been carried out here over a suitably chosen interpolation interval I s = [ s i , s f ] . The  I s interval is adjusted to obtain the best possible numerical resolution of Equation (3) over the entire time domain. This makes it possible to predict, with comparable accuracies, the behaviors of the solution of Equation (3) at both long and short times.

3.2.1. IIRFM Approach When α = 1 / 2 and λ 0

Situations where λ 0 are first considered. In this case, the denominator of the exact transfer function H L , 1 2 ( s , λ ) has no root with a positive real part, and therefore, H L , 1 2 ( s , λ ) has no discontinuity for s R + . The function G 1 2 ( s , λ ) = H L , 1 2 ( s , λ ) 1 is approximated in this case by a rational fraction noted as RF c M ( s , λ ) in the following form:
RF c M ( s , λ ) = i = 0 M 1 a M , i 1 2 , λ s i 1 + i = 1 M b M , i 1 2 , λ s i .
The 2 M coefficients a M , i , b M , i of the rational approximation RF c M ( s , λ ) are determined using a Chebyshev rational interpolation [12] of the function G 1 2 ( s , λ ) = H L , 1 2 ( s , λ ) 1 , on the interpolation interval I s = [ s i , s f ] , using a s k interpolation grid of the Gauss–Chebyshev type (composed of the zeros of first kind Chebyshev polynomials), as follows:
s k = 1 2 s i + s f + 1 2 s f s i cos 2 k 1 π 4 M , for k = 1 , 2 , , 2 M .
First-order partial fraction decomposition of (54) leads to writing the approximate transfer function H L , 1 2 M ( s , λ ) as a linear combination of first-order transfer functions of type (18), as follows:
H L , 1 2 M ( s , λ ) = 1 + i = 1 M a M , i 1 2 , λ b M , i 1 2 , λ + s ,
where the 2 M coefficients { a M , i , b M , i } are expressed as functions of the 2 M coefficients { a M , i , b M , i } .
As an example, an approximate expression of H 1 2 ( s , λ ) has been calculated as (56) for λ = 3 / 2 , with  s [ 0 , 5 ] and M = 8 . The coefficient values { a 8 , i , b 8 , i } obtained in this case are gathered in Table 3.
The relative error ε H r = | H L , 1 2 s , 3 2 / H L , 1 2 8 s , 3 2 1 | obtained when I s = [ 0 , 5 ] is shown as an insert to Figure 8. It can be seen that the rational approximation proposed here is very accurate over the I s domain and remains highly accurate for s [ 5 , 100 ] . This accuracy could be further improved by increasing the value of M, within the limits of stability of the transfer function H L , 1 2 M ( s , λ ) .

Monomial Source Term f ( t ) = t n

In the case of a monomial source term, f ( 0 ) = 0 induces u ( 0 ) = 0 . Unlike the case of the unit kernel studied in Section 3.1, there is no simple expression here for the analytical solution when α = 1 / 2 . Only the case n = 1 admits a relatively simple analytical solution, noted u λ , 1 2 , t , which will be used here to test the operation of the IIRFM approach as follows:
u λ , 1 2 , t ( t ) = e π λ 2 t 1 + erf λ π t 2 λ t 1 π λ 2 ,
where erf ( x ) is the error function. The numerical solution provided by the IIRFM is given by Equation (25), for which we can propose the following explicit expression when f ( t ) = t :
u λ , 1 2 , t IIR , M t k = t k + i = 1 M A M , i 2 t k T e + T e B M , i t k / T e 1 + B M , i B M , i T e + 2 t k B M , i 1 2 ,
where it is recalled that A M , i = a M , i T e / ( 2 + b M , i T e ) and B M , i = ( 2 b M , i T e ) / ( 2 + b M , i T e ) . The values of coefficients { a M , i , b M , i } are given in Table 3 when M = 8 .
Figure 9a shows the variations of u 3 2 , 1 2 , t ( t ) given by (57) when λ = 3 / 2 and u 3 2 , 1 2 , t IIR , 8 ( t ) given by (58), as functions of t, when f ( t ) = t , M = 8 , t [ 0 , 50 ] and T e = 0.05 .
The curves shown in Figure 9 were plotted using a Python code (version 3.8.5) calling on the numpy and scipy.special libraries, using single-precision floating-point numbers. In Figure 9a, expansion of the domain corresponding to the interval t [ 4.4 , 5 ] reveals anomalous discontinuities in the graphical representation of the analytical solution (57), whereas the numerical solution (Equation (58), dashed curve) provided by the IIRFM is continuous, as expected. Note that these discontinuities are also found in the graphical representation (see Figure 9b) of the function g λ ( t ) , defined by:
g λ ( t ) = 1 + erf λ π t .
The same numerical problem is observed with the free codes octave-8.4.0, scilab-2023.0.0 and julia-1.10.4 and other proprietary softwares used with single-precision floating-point numbers. The reason for this anomalous behavior lies in the internal evaluation of the function g λ ( t ) when λ π t 5 , as shown by Figure 9b. The Julia-1.10.4 code of Listing 1 graphically illustrates the pathological behavior of the function g λ ( t ) .
Listing 1. Julia-1.10.4 code illustrating the pathological behavior of g λ ( t ) .
using SpecialFunctions, Plots
function gL(t; La = -1.5)
     return 1 + erf(La*(pi*t)^0.5)
end
plot(t->gL(t, La = -1.5), 4.4, 5.5)
The IIRFM approach introduced in this work can be exploited to propose an approximation to the function g λ ( t ) that exhibits smooth numerical behavior in single precision, for any value of t 0 and for any value of λ , whether positive or negative. To this end, the function g λ ( t ) is first rewritten using expression (57):
g λ ( t ) = e π λ 2 t π λ 2 u λ , 1 2 , t ( t ) + 1 + 2 λ t .
By replacing the function u λ , 1 2 , t ( t ) by its IIRFM approximation (58), we finally obtain an approximate form of the function g λ ( t ) , noted here as g λ IIR , M :
g λ IIR , M ( t k ) = e π λ 2 t k π λ 2 u λ , 1 2 , t IIR , M ( t k ) + 1 + 2 λ t k ,
where t k = k T e 0 . Figure 9b shows the quality of the IIRFM approximation (60) for the g λ ( t ) function.
The present approach therefore provides an efficient numerical approximation of the g λ ( t ) function, which does not suffer in single precision from the anomalous discontinuities present in various free or proprietary codes.

Sinusoidal Source Term f ( t ) = sin ( t )

This is another source term that ensures a zero initial condition. An analytical expression u λ , 1 2 , f ( t ) of the solution also exists here:
u λ , 1 2 , f ( t ) = π λ 2 e π λ 2 t g λ ( t ) π 2 λ 4 + 1 + 2 π λ FS 2 π t π λ 2 sin ( t ) + cos ( t ) π 2 λ 4 + 1 2 π λ FC 2 π t + 1 π λ 2 cos ( t ) sin ( t ) π 2 λ 4 + 1 ,
where FC and FS are the Fresnel integrals, defined by FS ( x ) = 0 x sin π t 2 / 2 d t and FC ( x ) = 0 x cos π t 2 / 2 d t . The g λ ( t ) function with pathological numerical behavior in single precision is again found here. It has also been possible to determine here an explicit expression for the IIRFM numerical solution, as follows:
u λ , 1 2 , f IIR , M t k = sin t k + 2 cos T e 2 F M , T e ( t k ) ,
with F M , T e ( t k ) given by
F M , T e t k = i = 1 M A M , i B M , i + 1 B M , i t k / T e sin T e 2 B M , i sin t k + T e 2 + sin t k T e 2 B M , i cos T e 2 + sin 2 T e ,
where A M , i and B M , i have the same expressions as in previous examples.
Figure 10a compares the plots obtained using analytical (61) and numerical IIRFM (62) solutions. The problem of artificial discontinuities in the graphical representation of the analytical solution, due to the presence of the function g λ ( t ) in (61), is very pronounced here in the vicinity of t 5 . These defects are not present in the plot obtained by the IIRFM (dashed line), which here gives better numerical results in single precision than using the analytical expression.

Exponential Source Term f ( t ) = e t

For this last example, the case of a non-zero initial condition is considered, since here u λ , 1 2 , f ( 0 ) = f ( 0 ) = 1 . Using Equation (4), it is possible to find an analytical solution to the problem (3) when f ( t ) = e t and α = 1 / 2 , whose expression is as follows:
u λ , 1 2 , e t ( t ) = 2 λ DF t + π λ 2 e π λ 2 t g λ ( t ) + e t 1 π 3 / 2 λ 3 erfi t 1 + π λ 2 ,
where DF ( x ) = e x 2 0 x e y 2 d y is the Dawson function and erfi ( x ) = erf ( i x ) / i is the imaginary error function, with i 2 = 1 . The presence of the function g λ ( t ) in the analytical expression of the solution leads, as shown in Figure 10b, to considerable difficulties in graphically representing the function u λ , 1 2 , e t ( t ) given by Equation (64) in single precision.
When α 0 , applying the IIRFM to numerically solve Equation (3) assumes the initial condition on the solution to be null. If α 0 and the initial condition is non-zero but finite, then we introduce the auxiliary functions v λ , α , f 0 ( t ) and w λ , α , f ( 0 ) ( t ) , which satisfy the following auxiliary integral equations:
v λ , α , f 0 ( t ) = f 0 ( t ) + λ 0 t v λ , α , f 0 ( t ) ( t t ) α d t , w λ , α , f ( 0 ) ( t ) = f ( 0 ) + λ 0 t w λ , α , f ( 0 ) ( t ) ( t t ) α d t ,
with f 0 ( t ) = f ( t ) f ( 0 ) . The u λ , α , f ( t ) solution sought is written as u λ , α , f ( t ) = v λ , α , f 0 ( t ) + w λ , α , f ( 0 ) ( t ) , which satisfies Equation (3), with the initial condition u λ , α , f ( 0 ) = f ( 0 ) . Two important special cases should be noted: firstly, in the case where f ( 0 ) = 0 , the result is w λ , α , f ( 0 ) ( t ) = 0 , and the situation with zero initial condition studied in the previous paragraphs is retrieved. On the other hand, when f ( t ) = C = cste , in this case f 0 ( t ) = 0 results in v λ , α , f 0 ( t ) = 0 .
The calculation of v λ , α , f 0 ( t ) by the IIRFM has already been described in previous paragraphs, with f 0 ( 0 ) = 0 v λ , α , f 0 ( 0 ) = 0 . General result (24) is still applicable and it is also possible to find here an explicit expression of the IIRFM numerical solution for f ( t ) = e t f 0 ( t ) = e t 1 , which holds for α = 1 / 2 , as follows:
v λ , 1 2 , f 0 IIR , M t k = e t k 1 + i = 1 M A M , i 1 B M , i 1 + e T e e t k 2 1 e T e B M , i + B M , i t k / T e 1 + B M , i 1 e T e 1 B M , i 1 e T e B M , i .
There remains the need to explain the calculation of w λ , α , f ( 0 ) ( t ) . Since the source term f ( 0 ) is constant here, the auxiliary solution w λ , α , f ( 0 ) ( t ) can always be written as w λ , α , f ( 0 ) ( t ) = f ( 0 ) w λ , α 0 ( t ) , where the function w λ , α 0 ( t ) is independent of the expression of the source term f ( t ) considered in the integral Equation (3). This means that there is only one possible function w λ , α 0 ( t ) for a given set of parameters ( α , λ ) . More precisely, the function w λ , α 0 ( t ) is given by the following inverse Laplace transform:
w λ , α 0 ( t ) = L 1 1 s λ Γ 1 α s α .
Depending on the value taken by α , different approaches can be adopted to calculate w λ , α 0 ( t ) . In the case of the present paragraph where α = 1 / 2 , the inverse Laplace transform (67) can be computed explicitly and leads to the following analytical expression of w λ , 1 2 0 ( t ) :
w λ , 1 2 0 ( t ) = e λ 2 π t 1 + erf λ π t = e π λ 2 t g λ ( t ) .
Finally, the IIRFM solution of the initial problem is written here for f ( t ) = e t and α = 1 / 2 :
u λ , 1 2 , e t IIR , M t k = v λ , 1 2 , f 0 IIR , M t k + f ( 0 ) e π λ 2 t k g λ t k ,
where v λ , 1 2 , f 0 IIR , M is given by (66). This expression of u λ , 1 2 , e t IIR , M involves the function g λ ( t ) , whose anomalous numerical behavior in single precision has been highlighted in the previous paragraphs for λ π t 5 . Equation (69) is therefore acceptable as long as 0 t 1 π 5 λ 2 . When t > 1 π 5 λ 2 , we recommend not using the analytical expression of the function g λ ( t ) , but rather resorting to the following numerical form:
u λ , 1 2 , e t IIR , M t k = v λ , 1 2 , f 0 IIR , M t k + f ( 0 ) e π λ 2 t k g λ IIR , M t k .
Result (60) finally allows to write u λ , 1 2 , e t IIR , M t k in a fully numerical form:
u λ , 1 2 , e t IIR , M t k = v λ , 1 2 , f 0 IIR , M t k + f ( 0 ) π λ 2 u λ , 1 2 , t IIR , M t k + 1 + 2 λ t k ,
where v λ , 1 2 , f 0 IIR , M t k is given by (66) and u λ , 1 2 , t IIR , M t k is given by (58). This result is particularly interesting, since it holds (for the Abelian kernel α = 1 / 2 ) whatever the value of λ 0 and whatever the expression of f ( t ) , provided that f ( 0 ) is finite. Note that this result does not involve any numerical inversion of Laplace transforms, and that it generalizes the IIRFM for solving Abel integral equations of the second kind for α = 1 / 2 , to the case of non-zero but finite initial conditions.

Constant Term Source f ( t ) = C

In the present case, we consider f 0 ( t ) = f ( t ) f ( 0 ) = 0 and f ( 0 ) = C , leading to the following expressions for the auxiliary functions v λ , α , f 0 ( t ) and w λ , α , f ( 0 ) ( t ) :
v λ , α , f 0 ( t ) = λ 0 t v λ , α , f 0 ( t ) ( t t ) α d t = 0 , w λ , α , f ( 0 ) ( t ) = C + λ 0 t w λ , α , f ( 0 ) ( t ) ( t t ) α d t .
The analytical solution of the problem is deduced from the developments of Section “Exponential Source Term f ( t ) = e t ”:
u λ , 1 2 , C ( t ) = e π λ 2 t 1 + erf λ π t .
The expression for the IIRFM solution is deduced from Section “Exponential Source Term f ( t ) = e t ” for the present case, as follows:
u λ , 1 2 , C IIR , M t k = C π λ 2 u λ , 1 2 , t IIR , M t k + 1 + 2 λ t k .

3.2.2. IIRFM Approach When α = 1 / 2 and λ > 0

When λ > 0 , the denominator of the Laplace transfer function H α ( s , λ ) admits a positive root s λ , whose expression is s λ = λ Γ ( 1 α ) 1 / ( 1 α ) . The transfer function H L , α ( s , λ ) in this case presents a singularity at s = s λ . To allow a precise determination of the coefficients a M , i and b M , i of Equation (18), despite the singularity at s λ , we first calculate a Chebyshev rational interpolation (denoted by RF c M ) of the function ( s s λ ) × H L , α ( s , λ ) 1 . This interpolation is written as
RF c M ( s s λ ) × H L , α ( s , λ ) 1 = i = 0 M 1 a M , i α , λ s i 1 + i = 1 M 1 b M , i α , λ s i .
The following approximate form H L , α M ( s , λ ) of the transfer function H L , α ( s , λ ) is deduced:
H L , α M ( s , λ ) = 1 + i = 0 M 1 a M , i α , λ s i i = 0 M b M , i α , λ s i ,
with b M , 0 = s λ , b M , 1 = 1 b M , 1 s λ , b M , M = b M , M 1 and b M , m = b M , m 1 b M , m s λ for m [ 2 , M 1 ] . Further developments are strictly identical to those exposed in the previous paragraphs for λ < 0 and will therefore not be repeated here. Note also that all the results established in Section “Exponential Source Term f ( t ) = e t ” remain applicable here, but with λ > 0 .
As an example, let us consider the case where f ( t ) = e t , α = λ = 1 / 2 , for which it is found that s λ 0.7853982 . The rational approximation (75) has in this case been calculated over the interpolation interval I s = [ s λ , 5 s λ ] , for M = 4 . Figure 11a shows the analytical (64) and numerical IIRFM (71) solutions in the present case, for t [ 0 , 3 ] and T e = 5 × 10 3 . Figure 11b shows the relative error ε r = 100 × ( u 1 2 , 1 2 , e t IIR , 4 / u 1 2 , 1 2 , e t 1 ) in %, as a function of t. The correspondence between analytical and IIRFM resolutions is excellent, even for the low M value considered here.

3.2.3. Comparison of the IIRFM with the Mouley et al. [13] Approach

In this section, the IIRFM is applied to the numerical solution of two second-kind Abel linear integral equations, which Mouley et al. have studied using a recent numerical method based on an approximation involving Daubechies wavelets [13]. The accuracy of the results provided by the two methods is compared on two examples, for different orders of approximation (see Table 4 and Table 5).

First Example

Consider the following Abel integral equation of the second kind (there is a typographical error in the paper by Mouley et al. where 16 5 t 5 / 2 is found in place of 16 15 t 5 / 2 ):
u ( t ) = t 2 + 16 15 t 5 2 0 t u ( t ) t t d t ,
which admits as analytical solution the function u ( t ) = t 2 . This is a convolutive integral equation of the form (3), with α = 1 / 2 , λ = 1 , f ( t ) = t 2 + 16 15 t 5 / 2 , and u ( 0 ) = 0 . The developments of Section 3.2.1 are therefore directly applicable to the numerical solution of this integral equation using the IIRFM method.
Table 4 gathers the absolute errors ε Dau j = u ( t ) u Dau j ( t ) committed using the Mouley et al. method and absolute errors ε IIR M = u ( t ) u IIR M ( t ) committed using the IIRFM (with I s = [ 0 , 15 ] and T e = 10 3 ), for different values of t [ 0 , 1 ] , different resolutions j, and for different approximation orders M.
Table 4. Absolute errors in the numerical resolution of (77) using the Mouley et al. [13] and IIRFM methods (with I s = [ 0 , 15 ] and T e = 10 3 ).
Table 4. Absolute errors in the numerical resolution of (77) using the Mouley et al. [13] and IIRFM methods (with I s = [ 0 , 15 ] and T e = 10 3 ).
ε Dau j × 10 6 ε IIR M × 10 6
t Exact Values j = 4 j = 6 j = 8 M = 4 M = 6 M = 8
0.125 0.015625 11710100
0.25 0.062500 3710200
0.375 0.140625 22001200
0.5 0.250000 1600310
0.625 0.390625 12002310
0.75 0.562500 10003710
0.875 0.765625 8002110
It can be seen from these results that the accuracies obtained with the IIRFM are, in the case of integral Equation (77), quite comparable with those provided by the method of Mouley et al. [13].

Second Example

Let us now consider the following Abel integral equation of the second kind, also studied by Mouley et al. in [13]:
u ( t ) = 1 t + 1 + π 8 1 4 sin 1 1 t 1 + t 1 4 0 t u ( t ) t t d t ,
of which the analytical solution is u ( t ) = 1 / t + 1 . This is still a convolutive integral equation of the form (3), with α = 1 / 2 , λ = 1 / 4 , f ( t ) = 1 / t + 1 + π 8 1 4 sin 1 1 t 1 + t , and u ( 0 ) = 1 .
Table 5 gathers the absolute errors ε Dau j and ε IIR M committed by numerically calculating the solution of (78), for different values of t [ 0 , 1 ] , for different resolutions j, and different approximation orders M.
Table 5. Absolute errors ε in the numerical resolution of (78) using the Mouley et al. [13] and IIRFM methods (with I s = [ 0 , 15 ] and T e = 10 3 ).
Table 5. Absolute errors ε in the numerical resolution of (78) using the Mouley et al. [13] and IIRFM methods (with I s = [ 0 , 15 ] and T e = 10 3 ).
ε Dau j × 10 6 ε IIR M × 10 6
t 1 / t + 1 j = 4 j = 6 j = 8 M = 4 M = 6 M = 8
0.125 0.942809 21,732437010742933
0.25 0.894427 10,73927746832831
0.375 0.852803 856820915151511
0.5 0.816497 697116954173211
0.625 0.784465 58901433353921
0.75 0.755929 511312443072711
0.875 0.730297 452211012715410
This time, it was found that the accuracy of the numerical resolution of (78) by the IIRFM (with I s = [ 0 , 15 ] and T e = 10 3 ) is much better than the accuracy allowed by the method of Mouley et al., which uses Daubechy wavelets. Despite the greater complexity of the source function f ( t ) considered in this second example, the accuracy obtained here with the IIRFM remains comparable to that obtained with the same method for the previous example. This is a clear advantage of the IIRFM in the present case, as it is found that the accuracy depends relatively little on the complexity of the source function f ( t ) .

3.2.4. Comparison of the IIRFM with the Singha et al.  [14] Approach

In the present section, the results obtained with the IIRFM are compared with the numerical results obtained using Laguerre polynomials L i of order i [14], at approximation order M = 5 :
u S M ( t ) = i = 0 M i L i ( t ) ,
where u S M ( t ) is the approximation of Singha et al.
The accuracy of the results provided by the two methods is compared in two examples proposed by Singha et al. in [14].

First Example

Consider the following Abel integral equation of the second kind:
u ( t ) = π t 2 + t 0 t u t t t d t ,
of which analytical solution is u ( t ) = t . This is a convolutive integral equation of the form (3), with α = 1 / 2 , λ = 1 , f ( t ) = π t 2 + t , and u ( 0 ) = 0 . Table 6 gathers the absolute errors ε committed by numerically calculating the solution of (80), with the polynomial method proposed by Singha et al. [14] and with the IIRFM, for different values of t [ 0 , 2 ] , and for the same order of approximation M = 5 .
The results in Table 6 clearly demonstrate the benefits of the IIRFM for the numerical resolution of (80). For each value of t considered, the IIRFM is always much more accurate here than the approach proposed by Singha et al. [14].
The values of the Laguerre approximation u S M ( t ) = i = 0 M i L i ( t ) have been calculated here using the { i } coefficient values of Singha et al. [14]:
{ 0 , 1 , 2 , 3 , 4 , 5 } = π × 1 2 , 1 4 , 1 16 , 1 32 , 5 256 , 7 512 .

Second Example

Consider the following Abel integral equation of the second kind:
u ( t ) = π Γ 5 4 Γ 7 4 t 3 4 + t 1 4 0 t u t t t d t ,
which admits the analytical solution u ( t ) = t 1 / 4 . This is a convolutive integral equation of the form (3), with α = 1 / 2 , λ = 1 , f ( t ) = π Γ 5 4 Γ 7 4 t 3 4 + t 1 4 and u ( 0 ) = 0 . Table 7 gathers the absolute errors ε committed by numerically calculating the solution of (78) with the polynomial method of Singha et al. [14] and with the IIRFM, for different values of t [ 0 , 2 ] , and for the same order of approximation M = 5 .
It can also be seen here that the numerical resolution of Equation (82) by the IIRFM remains overall more accurate than the resolution made with the numerical method of Singha et al., and for which the absolute error suffers significant variations over the range of t values studied.
The values of the Laguerre approximation u S M ( t ) = i = 0 M i L i ( t ) have been calculated here using the { i } coefficient values of Singha et al. [14]:
{ 0 , 1 , 2 , 3 , 4 , 5 } = 0.906402 , 0.226601 , 0.084975 , 0.049569 , 0.0340787 , 0.025559

3.3. Logarithmic Kernel

We now consider the case of a logarithmic kernel K ( t , t ) = ln t t , for which K ¯ ( s ) = [ γ + ln ( s ) ] / s . The Laplace transfer function H L ( s , λ ) associated with Equation (2) becomes
H L ( s , λ ) = 1 1 + λ γ + ln ( s ) / s ,
where γ 0.577216 is the Euler constant. The transfer function H L ( s , λ ) in this case admits a pole s λ , which is the root of the non linear equation e γ = s e s / λ , given by the Lambert W function: s λ = λ W e γ / λ . The same process is used here as in Section 3.2.2, for which it is recalled that the first step consists of calculating a rational Chebyshev interpolation (noted RF c M ) of the function ( s s λ ) × H ( s , λ ) 1 .
In the following paragraphs, two examples of IIRFM solving Equation (2) with logarithmic kernel are presented when λ = 1 . Refer to Section 3.2.1 and Section 3.2.2 for details of the different steps involved in implementing the IIRFM in the present case.

3.3.1. First Example: f ( t ) = t + 1 4 t 2 [ 3 2 ln ( t ) ]

Equation (2) admits in the present case the following analytical solution u 1 , f ( t ) = t , with zero initial condition.
Figure 12 compares the numerical results obtained with the IIRFM (for an order of approximation M = 5 ) with the results obtained with the analytical solution. The relative error shown in Figure 12b here remains quite small, especially for increasing values of t, attesting to the stability of the IIRFM in this case, despite the presence of a pole in the transfer function H L ( s , λ ) .

3.3.2. Second Example: f ( t ) = 1 e t / τ + λ t + ( τ t ) ln ( t ) +   λ τ e t / τ γ ln ( τ ) + i π + t / τ e x x d x

Equation (2) admits in the present case the following analytical solution u 1 , f ( t ) = 1 e t / τ . Figure 12 compares the numerical results obtained with the IIRFM (for an order of approximation M = 5 ) with the results obtained with the analytical solution. Despite the relative complexity of the f ( t ) source function considered in this example, the IIRFM provides stable and accurate numerical results, as illustrated by the variations in relative error shown in Figure 13b.

4. Basic Applications in Thermics

To conclude this study, the IIRFM is applied to the numerical solution of two basic thermal problems: the study of the temperature time evolution of a source point placed in a homogeneous infinite solid medium (see Section 4.1), and the study of the time evolution of the surface temperature of a semi-infinite solid when radiative heat transfers can be neglected (see Section 4.2).

4.1. Heat Point Source

Consider a punctual heat source placed at a point M of Cartesian coordinates ( x , y , z ) , within a stationary, homogeneous, and opaque solid medium (m) of constant physical properties. At each time t, the temperature T ( x , y , z , t ) of the medium at a point M with coordinates ( x , y , z ) is written as T ( M , t ) = T 0 + v ( M , t ) , with v ( M , 0 ) = 0 K. Noting Q ˙ ( M , t ) the thermal power emitted by the source point placed at M , and the temperature variation v ( M , t ) is given by the following integral [25]:
v ( M , t ) = 1 8 ρ c π κ 3 / 2 0 t Q ˙ ( M , t ) e r 2 4 κ ( t t ) ( t t ) 3 / 2 d t ,
where r 2 = ( x x ) 2 + ( y y ) 2 + ( z z ) 2 ; κ = k / ρ c > 0 is the thermal diffusivity of (m); k, ρ , and c are the thermal conductivity, density, and specific heat of (m), respectively. It is also assumed that thermal power is generated by the Joule effect and is written as Q ˙ ( M , t ) = R ( M , t ) i 2 ( t ) , where R ( M , t ) is the electrical resistance of the source point and i ( t ) is the electrical intensity of the heating current. In the present example, the resistance of the source point is assumed to vary linearly with temperature, as is usually the case for examples with metallic materials:
R ( M , t ) = R 0 1 + α T v ( M , t ) ,
where α T 0 is the temperature coefficient of the source point. From Equation (85), the temperature v ( M , t ) is given in this case by the following relation:
v ( M , t ) = R 0 8 ρ c π κ 3 / 2 0 t i 2 ( t ) e r 2 4 κ ( t t ) ( t t ) 3 / 2 d t + α T R 0 8 ρ c π κ 3 / 2 0 t i 2 ( t ) v ( M , t ) e r 2 4 κ ( t t ) ( t t ) 3 / 2 d t .
Assuming that the point source is spherical, with a radius being a very small compared to the dimensions of the medium (m), located at the origin of the coordinates, the uniform temperature v ( t ) of the point source verifies the following integral equation deduced from (87):
v ( t ) = R 0 8 ρ c π κ 3 / 2 0 t i 2 ( t ) e τ a ( t t ) ( t t ) 3 / 2 d t + α T R 0 8 ρ c π κ 3 / 2 0 t i 2 ( t ) v ( t ) e τ a ( t t ) ( t t ) 3 / 2 d t ,
where τ a = a 2 / 4 κ . Assuming that the heating current of the point source is of constant intensity i ( t ) = I 0 , the temperature v ( t ) is now given by the following linear Volterra integral equation of the second kind:
v ( t ) = R 0 I 0 2 8 ρ c π κ 3 / 2 0 t e τ a ( t t ) ( t t ) 3 / 2 d t + α T R 0 I 0 2 8 ρ c π κ 3 / 2 0 t v ( t ) e τ a ( t t ) ( t t ) 3 / 2 d t ,
which is of the form (2), with:
λ = α T R 0 I 0 2 8 ρ c π κ 3 / 2 > 0 ; f ( t ) = λ α T 0 t e τ a ( t t ) ( t t ) 3 / 2 d t ; K ( t , t ) = e τ a ( t t ) ( t t ) 3 / 2 .
To the best of our knowledge, integral Equation (89) has no analytical solution, so we propose to solve it numerically using the IIRFM, following the various steps shown in Figure 2.
Using the result L [ e τ a / t / t 3 / 2 ] = π / τ a e 2 s τ a , the Laplace transfer function H L ( s , λ ) of the temperature v ( t ) is written as follows, according to (17):
H L ( s , λ ) = 1 1 λ π / τ a e 2 s τ a .
Note that the denominator of the transfer function (91) can have a pole s λ , given by the following relationship:
s λ = ln λ π / τ a 2 4 τ a when 1 λ τ a π < 1 .
The transfer function (91) is then decomposed into the form (18), using a Chebyshev rational interpolation. Applying the bilinear transformation finally leads to the result (24), which is written here:
v λ , f IIR , M p T e = f p T e + i = 1 M A M , i j = 0 p 1 B M , i p ( j + 1 ) f ( j T e ) + f ( ( j + 1 ) T e ) ,
with:
f ( j T e ) = λ α T 0 j T e e τ a ( j T e t ) ( j T e t ) 3 / 2 d t = R 0 I 0 2 4 π k a erfc τ a j T e ,
where the coefficients A M , i and B M , i have the usual expressions introduced in the previous paragraphs. Note that in the limit where α T 0 (case of constant electrical resistance), the temperature v ( t ) is then identified with the function f ( t ) .
Figure 14 shows the time variations in temperature v ( t ) of a heat point source with radius a = 0.5 mm, immersed in ethylene glycol at rest, at room temperature ( ρ = 1113.5 kg / m 3 , k = 0.254 W / m / K , c = 2380 J / K / kg ). The electrical heating current considered had a constant intensity I 0 = 250 mA and the temperature coefficient of the point source was α T = 0.0039 K 1 , for an electrical resistance R 0 = 1 Ω . The Chebyshev interpolation H L M ( s , λ ) was calculated with M = 4 , over the interpolation interval I s = [ 0 , 5 ] Hz 1 . The maximal relative error ε r max committed on the approximation of H L ( s , λ ) was in this case ε r max = 0.058 % for s [ 10 2 , 10 2 ] Hz 1 . The sampling period was T e = 0.01 s.
From the results shown in Figure 14, it can be seen that the temperature coefficient α T has a non-negligible influence on the time evolution of the source point temperature. When the influence of α T is taken into account, the self-heating of the source point leads to an increase in its electrical resistance (because α T is positive), and so at constant heating intensity, this leads to an increase in the Joule heating power and consequently in its temperature v ( t ) , compared to the situation where α T = 0 . These observations could be applied to various methods of the thermal characterization of materials.

4.2. Surface Temperature of Semi-Infinite Solids

The IIRFM is now applied to the solution of the linear integral equation verified by the interface temperature of a semi-infinite isotropic solid (defined by x > 0 ), of uniform initial temperature T i , subjected to mixed (or Robin, see Figure 15a) boundary conditions at x = 0 .
The temperature field θ ( x , t ) = T ( x , t ) T i in this case satisfies the usual heat equation without a source term, as follows:
θ t = κ 2 θ x 2 with θ ( x , 0 ) = 0 and θ ( + , t ) < ,
where κ = k / ρ c > 0 is the thermal diffusivity of the solid, assumed constant. Applying the Laplace transformation to Equation (95) leads to the ordinary differential equation d x 2 θ ¯ ( x , s ) s κ θ ¯ ( x , s ) = 0 , of which the solution is written in the general form θ ¯ ( x , s ) = A e s / κ x + B e s / κ x , with θ ¯ ( x , s ) = L θ ( x , t ) . The finite boundary condition θ ( + , t ) < implies B = 0 here.
Let us note q e ( 0 , t ) = q e ( t ) n e 0 , the heat flux exchanged by the solid with the external fluid at x = 0 . The continuity of heat fluxes at the interface x = 0 is written as q e ( 0 , t ) = k n e 0 T n e 0 = k n e 0 θ n e 0 . Applying the Laplace transformation to the continuity relation leads to the following expression for θ ¯ ( x , s ) :
θ ¯ ( x , s ) = κ k q ¯ e ( s ) e s / κ x s ,
where q ¯ e ( s ) = L q e ( t ) . By recalling that L 1 e β s / s = e β 2 / 4 t / π t and applying the convolution theorem, it is possible to invert θ ¯ ( x , s ) and write θ ( x , t ) in the form of the following convolution integral:
θ ( x , t ) = 1 k κ π 0 t q e ( τ ) ( t τ ) 1 / 2 e x 2 4 κ ( t τ ) d τ .
Finally, using Newton’s law q e ( 0 , t ) = h θ 0 , t θ ( t ) n e 0 leads to the following expression of θ ( x , t ) :
θ ( x , t ) = h k κ π 0 t θ ( τ ) θ ( 0 , τ ) ( t τ ) 1 / 2 e x 2 4 κ ( t τ ) d τ ,
where the Newton coefficient h is assumed to be constant. According to result (98), the surface temperature θ ( 0 , t ) finally verifies the following Abel integral equation of the second kind:
θ ( 0 , t ) = f ( t ) + λ 0 t θ ( 0 , τ ) ( t τ ) 1 / 2 d τ ,
with λ = h k κ π < 0 and f ( t ) = h k κ π 0 t θ ( τ ) ( t τ ) 1 / 2 d τ .
In the following paragraphs, Equation (99) is solved by the IIRFM for two different expressions of the temperature T ( t ) of the fluid far from the solid wall.

4.2.1. First Application: T ( t ) = cste

First, consider the usual case where the temperature T of the fluid far from the solid wall x = 0 is constant. In this case, it is found that f ( t ) = 2 h k θ κ / π t = f t n with n = 1 / 2 and θ = T T i . Equation (99) then admits an analytical solution θ a ( t ) = T ( 0 , t ) T i involving the complementary error function:
θ a ( t ) θ = 1 exp h 2 κ t k 2 erfc h κ t k .
Equation (99) is solved in the following by the IIRFM, using the general result (24). First of all, it is possible here to determine an explicit expression Ξ N , i for the sum j = 0 k 1 B N , i k ( j + 1 ) f ( j T e ) + f ( ( j + 1 ) T e ) , in the general case where f ( t ) = f t n :
Ξ N , i k , T e , B N , i , f = f T e n B N , i k 1 1 + B N , i Li n B N , i 1 + k n ( 1 + B N , i 1 ) Φ B N , i 1 , n , k ,
where the polylogarithm Li and transcendent Lerch Φ functions were defined in Section 3.1.1. Using result (101) with n = 1 / 2 , the following explicit expression is obtained for the IIRFM solution to the problem (99), at time t k = k T e (with k 1 ):
θ t f IIR , N t k f T e = k + i = 1 N A N , i B N , i k 1 1 + B N , i Li 1 2 B N , i 1 + k ( 1 + B N , i 1 ) Φ B N , i 1 , 1 2 , k .
The thermal problem (99), of calculating the surface temperature of a semi-infinite solid surrounded by a fluid, is solved using result (102) for the following physical properties: k = 40 W · K 1 · m 1 ; ρ = 7850 kg · m 3 ; c = 490 J · K 1 · kg 1 ; κ = 1.04 · 10 5 m 2 · s 1 ; for Newton coefficient h = 100 W · K 1 · m 2 ; t [ 0 , 200 ] s ; and T e = 1 s. Numerical (continuous line) and analytical results ( + symbols) are shown in Figure 15b. From this figure it can be seen that there is a good correspondence between the IIRFM numerical solution and the analytical solution. This good agreement is confirmed by the low values of the maximal error ε , N t IIR = 1.41 · 10 4 K and root mean square error ε 2 , N t IIR = 1.58 · 10 5 K, with
ε 2 , N t IIR = 1 N t + 1 k = 0 N t θ a ( t k ) θ t f IIR , N t k 2 and ε , N t IIR = max S θ a θ t f IIR , N ,
where S is the set of points t k , θ t f IIR , N t k n = 0 N t obtained using the IIRFM method.

4.2.2. Second Application: T ( t ) = T ^ sin ( ω t )

Thermal problem (99) is still considered, but now for a sinusoidal temperature T ( t ) with pulsation ω : T ( t ) = T ^ sin ( ω t ) , all other things being equal to what was considered in the previous example. The f ( t ) function in Equation (99) becomes here
f ( t ) = λ θ ^ 0 t sin ( ω τ ) ( t τ ) 1 / 2 d τ , = λ θ ^ π ω FC 2 ω t π sin ( ω t ) FS 2 ω t π cos ( ω t ) ,
where λ = h k κ / π ; θ ^ = T ^ T i ; FC and FS are the Fresnel integrals. An analytical solution θ a ( t ) also exists here, and can be obtained by inverting the Laplace transform θ ¯ ( 0 , s ) = θ ^ ω s 2 + ω 2 / 1 + k s h κ :
θ a ( t ) θ ^ = h 4 κ 2 sin ( ω t ) + h 2 k 2 κ ω cos ( ω t ) e h 2 κ t k 2 erfc h κ t k h 4 κ 2 + k 4 ω 2 + h k 2 κ ω h 4 κ 2 + k 4 ω 2 FC 2 ω t π k 2 ω sin ( ω t ) h 2 κ cos ( ω t ) h k 2 κ ω h 4 κ 2 + k 4 ω 2 FS 2 ω t π h 2 κ sin ( ω t ) + k 2 ω cos ( ω t ) .
Figure 16 compares the numerical results (solid line) obtained by the IIRFM to the results provided (black dashes) by the analytical expression (105), when ω = 2 π / T h with T h = 12 h and t [ 0 , t f = 6 T h ] . The correspondence between the results provided by the IIRFM and the analytical solution is again very good, with maximum error ε , N t IIR = 3.72 · 10 3 K and root mean square error ε 2 , N t IIR = 1.25 · 10 3 K, calculated for a sampling period T e = 216 s .

5. Conclusions

In this work, the first-order infinite impulse response filters method (IIRFM), usually reserved for the study of linear ordinary differential equations with constant coefficients, has been extended to the numerical solution of various linear Volterra convolution integral equations (VCIE) of the second kind. A wide variety of examples have been studied, illustrating the great flexibility of the IIRFM. Three types of convolutive kernels were discussed in detail in this work: unit, Abel, and logarithmic kernels.
In the case of the unit kernel, the IIRFM was compared to the homotopic perturbation method with Laplace transform (HPM-L). The HPM-L method has the great advantage of leading, for certain special cases, to an analytic solution of the linear convolutive Volterra integral equation of the second kind. However, due to its iterative nature, the numerical version of the HPM-L method may exhibit convergence speeds that are far too low to allow, in the absence of an analytical solution, an efficient numerical resolution of certain VCIEs. In contrast, we have observed that the IIRFM leads, in the case of the unit kernel, to numerical solutions that do not exhibit convergence problems. In many cases, it has also been possible to provide a simple, highly efficient analytical expression for the IIRFM numerical solution.
In the case of the Abel kernel, the numerical results provided by the IIRFM approach have been compared with the results provided by various recent numerical methods (Daubechies wavelets and approximations based on Laguerre polynomials). For most of the comparisons examined, the IIRFM proved to be the best-performing and most consistent overall.
Finally, the IIRFM has been applied to the numerical solution of Volterra integral equations of the second kind encountered in two linear thermal problems. In addition, a new approximate expression for the pathological function g λ ( t ) = 1 + erf ( λ π t ) has been proposed, giving accurate results for λ < 0 , even using single-precision floating-point numbers.
In a forthcoming paper, the IIRFM will be extended to the numerical solution of the non-linear Volterra convolution integral equations of the second kind. It will also be very interesting to consider the application of the IIRFM to fractional differential or integral equations [26].

Funding

This research received no external funding.

Data Availability Statement

The original contributions presented in this study are included in the article, further inquiries can be directed to the corresponding author.

Conflicts of Interest

The author declares no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:
BLTBilinear Transformation
CAIEConvolutive Abel Integral Equation
CIEConvolutive Integral Equation
CVIEConvolutive Volterra Integral Equation
FODFirst-Order Partial Fractions Decompositon
HPMHomotopy Perturbation Method
HPM-LHomotopy Perturbation Method with Laplace Transform
IIRInfinite Impulse Response
IIRFMInfinite Impulse Response First-Order Filters Method
ODEOrdinary Differential Equation
OIDEOrdinary Integro-Differential Equation
LTLaplace Transformation
LTISLinear Time Invariant Systems
PDEPartial Differential Equation

References

  1. De, S.; Mandal, B.; Chakrabarti, A. Use of Abel integral equations in water wave scattering by two surface-piercing barriers. Wave Motion 2010, 47, 279–288. [Google Scholar] [CrossRef]
  2. Mirceski, V.; Tomovski, Z. Analytical solutions of integral equations for modelling of reversible electrode processes under voltammetric conditions. J. Electroanal. Chem. 2008, 619–620, 164–168. [Google Scholar] [CrossRef]
  3. Wazwaz, A.M. Linear and Nonlinear Integral Equations—Methods and Applications; Springer: Berlin/Heidelberg, Germany, 2011. [Google Scholar]
  4. Polyanin, A.D.; Manzhirov, A.V. Handbook of Integral Equations, 2nd ed.; Chapman and Hall/CRC: Boca Raton, FL, USA, 2008. [Google Scholar]
  5. Kumar, S.; Kumar, A.; Kumar, D.; Singh, J.; Singh, A. Analytical solution of Abel integral equation arising in astrophysics via Laplace transform. J. Egypt. Math. Soc. 2015, 23, 102–107. [Google Scholar] [CrossRef]
  6. Thota, S. Solution of Generalized Abel’s Integral Equations by Homotopy Perturbation Method with Adaptation in Laplace Transformation. Sohag J. Math. 2022, 9, 29–35. [Google Scholar]
  7. Chowdhury, M.S.H.; Mohamed, M.S.; Gepreel, K.A.; Al-Malki, F.A.; Al-Humyani, M. Approximate Solutions of the Generalized Abel’s Integral Equations Using the Extension Khan’s Homotopy Analysis Transformation Method. J. Appl. Math. 2015, 2015, 357861. [Google Scholar]
  8. Liao, S. On the homotopy analysis method for nonlinear problems. Appl. Math. Comput. 2004, 147, 499–513. [Google Scholar] [CrossRef]
  9. Liao, S. Comparison between the homotopy analysis method and homotopy perturbation method. Appl. Math. Comput. 2005, 169, 1186–1194. [Google Scholar] [CrossRef]
  10. Bairwa, R.K.; Kumar, A.; Kumar, D. An Efficient Computation Approach for Abel’s Integral Equations of the Second Kind. Sci. Technol. Asia 2020, 25, 85–94. [Google Scholar]
  11. Heyd, R. Real-time heat conduction in a self-heated composite slab by Padé filters. Int. J. Heat Mass Transf. 2014, 71, 606–614. [Google Scholar] [CrossRef]
  12. Lahboub, D.; Heyd, R.; Bakak, A.; Lotfi, M.; Koumina, A. Solution of Basset integro-differential equations by IIR digital filters. Alex. Eng. J. 2022, 61, 11899–11911. [Google Scholar] [CrossRef]
  13. Mouley, J.; Panja, M.M.; Mandal, B.N. Approximate solution of Abel integral equation in Daubechies wavelet basis. CUBO A Math. J. 2021, 23, 245–264. [Google Scholar] [CrossRef]
  14. Singha, N.; Nahak, C. Solutions of the Generalized Abel’s Integral Equation using Laguerre Orthogonal Approximation. Appl. Appl. Math. Int. J. (AAM) 2019, 14, 1051–1066. [Google Scholar]
  15. Adomian, G.; Rach, R. Noise terms in decomposition series solution. Comput. Math. Appl. 1992, 24, 61–64. [Google Scholar] [CrossRef]
  16. Adomian, G. Solving Frontier Problems of Physics: The Decomposition Method; Kluwer: Alphen aan den Rijn, The Netherlands, 1994. [Google Scholar]
  17. He, J.H. Homotopy perturbation technique. Comput. Methods Appl. Mech. Eng. 1999, 178, 257–262. [Google Scholar] [CrossRef]
  18. He, J.H. A coupling method of a homotopy technique and a perturbation technique for non-linear problems. Int. J. Non-Linear Mech. 2000, 35, 37–43. [Google Scholar] [CrossRef]
  19. He, J.H. Homotopy perturbation method: A new nonlinear analytical technique. Appl. Math. Comput. 2003, 135, 73–79. [Google Scholar] [CrossRef]
  20. Madani, M.; Fathizadeh, M. Homotopy Perturbation Algorithm using Laplace Transformation. Nonlinear Sci. Lett. A 2010, 1, 263–267. [Google Scholar]
  21. Khan, Y.; Wu, Q. Homotopy perturbation transform method for nonlinear equations using He’s polynomials. Comput. Math. Appl. 2011, 61, 1963–1967. [Google Scholar] [CrossRef]
  22. Singh, J.; Kumar, D.; Rathore, S. Application of Homotopy Perturbation Transform Method for Solving Linear and Nonlinear Klein-Gordon Equations. J. Inf. Comput. Sci. 2012, 7, 131–139. [Google Scholar]
  23. Heyd, R.; Hadaoui, A.; Saboungi, M. 1D analog behavioral SPICE model for hot wire sensors in the continuum regime. Sens. Actuators A Phys. 2012, 174, 9–15. [Google Scholar] [CrossRef]
  24. Lotfi, M.; Mezrigui, L.; Heyd, R. Study of heat conduction through a self-heated composite cylinder by Laplace transfer functions. Appl. Math. Model. 2016, 40, 10360–10376. [Google Scholar] [CrossRef]
  25. Carslaw, H.; Jaeger, J.C. Conduction of Heat in Solids; Oxford University Press: Oxford, UK, 1959. [Google Scholar]
  26. Zhao, K.; Liu, J.; Lv, X. A Unified Approach to Solvability and Stability of Multipoint BVPs for Langevin and Sturm–Liouville Equations with CH-Fractional Derivatives and Impulses via Coincidence Theory. Fractal Fract. 2024, 8, 111. [Google Scholar] [CrossRef]
Figure 1. Differential or integral, linear, and invariant problems. (a) Schematic representation in direct space. (b) Schematic representation in Laplace space.
Figure 1. Differential or integral, linear, and invariant problems. (a) Schematic representation in direct space. (b) Schematic representation in Laplace space.
Mathematics 12 02416 g001
Figure 2. Main steps of the IIRFM, as applied to the numerical solution of linear CVIEs ( a L = 1 and a NL = 0 ) of the second kind, with LT: Laplace transform; FOD: first-order partial fractions decomposition; BLT: bilinear transformation; IIRFM: final result of the IIRFM method.
Figure 2. Main steps of the IIRFM, as applied to the numerical solution of linear CVIEs ( a L = 1 and a NL = 0 ) of the second kind, with LT: Laplace transform; FOD: first-order partial fractions decomposition; BLT: bilinear transformation; IIRFM: final result of the IIRFM method.
Mathematics 12 02416 g002
Figure 3. Dashed lines show the evolution of v N H ( t ) as a function of t for n = 1 , different values of λ , and different values of N H : (a) λ = 1 and N H = 3 , 5 , 7 , 9 ; (b) λ = 1.5 and N H = 5 , 9 , 13 , 17 .
Figure 3. Dashed lines show the evolution of v N H ( t ) as a function of t for n = 1 , different values of λ , and different values of N H : (a) λ = 1 and N H = 3 , 5 , 7 , 9 ; (b) λ = 1.5 and N H = 5 , 9 , 13 , 17 .
Mathematics 12 02416 g003
Figure 4. Relative errors as a function of t, when α = 0 , f ( t ) = t , and λ = 1.5 : (a) ε r HPM - L ( t ) with N H = 15 and T e = 10 2 ; (b) ε r IIRFM ( t ) with T e = 10 2 .
Figure 4. Relative errors as a function of t, when α = 0 , f ( t ) = t , and λ = 1.5 : (a) ε r HPM - L ( t ) with N H = 15 and T e = 10 2 ; (b) ε r IIRFM ( t ) with T e = 10 2 .
Mathematics 12 02416 g004
Figure 5. Relative errors as a function of t, when K ( t , t ) = 1 , f ( t ) = t , and λ = 1.5 : (a) ε r HPM - L ( t ) with N H = 15 and T e = 10 2 ; (b) ε r IIRFM ( t ) with T e = 10 2 .
Figure 5. Relative errors as a function of t, when K ( t , t ) = 1 , f ( t ) = t , and λ = 1.5 : (a) ε r HPM - L ( t ) with N H = 15 and T e = 10 2 ; (b) ε r IIRFM ( t ) with T e = 10 2 .
Mathematics 12 02416 g005
Figure 6. Absolute error ε k ( T e ) as a function of T e , when K ( t , t ) = 1 , f ( t ) = t , λ = 1 , and t k = 1 .
Figure 6. Absolute error ε k ( T e ) as a function of T e , when K ( t , t ) = 1 , f ( t ) = t , λ = 1 , and t k = 1 .
Mathematics 12 02416 g006
Figure 7. IIRFM and HPM-L numerical results when K ( t , t ) = 1 , f ( t ) = e t , T e = 10 2 , and N H = 5 , 15 , and 20: (a) λ = 5 and (b) λ = 5 .
Figure 7. IIRFM and HPM-L numerical results when K ( t , t ) = 1 , f ( t ) = e t , T e = 10 2 , and N H = 5 , 15 , and 20: (a) λ = 5 and (b) λ = 5 .
Mathematics 12 02416 g007
Figure 8. Exact transfer function H L , 1 2 s , 3 2 and its Chebyshev rational approximation H L , 1 2 8 s , 3 2 calculated using the interpolation interval I s = [ 0 , 5 ] . The insert shows the relative error ε H r committed on the rational approximation for s [ 0 , 100 ] .
Figure 8. Exact transfer function H L , 1 2 s , 3 2 and its Chebyshev rational approximation H L , 1 2 8 s , 3 2 calculated using the interpolation interval I s = [ 0 , 5 ] . The insert shows the relative error ε H r committed on the rational approximation for s [ 0 , 100 ] .
Mathematics 12 02416 g008
Figure 9. (a) Analytical solution (57) and numerical IIRFM solution (58) when f ( t ) = t , α = 1 / 2 , λ = 3 / 2 , M = 8 , and T e = 0.05 . (b) Function g λ ( t ) = 1 + erf λ π t when λ = 3 / 2 and t [ 4.4 , 5 ] .
Figure 9. (a) Analytical solution (57) and numerical IIRFM solution (58) when f ( t ) = t , α = 1 / 2 , λ = 3 / 2 , M = 8 , and T e = 0.05 . (b) Function g λ ( t ) = 1 + erf λ π t when λ = 3 / 2 and t [ 4.4 , 5 ] .
Mathematics 12 02416 g009
Figure 10. (a) Analytical solution (61) and numerical IIRFM solution (62) when f ( t ) = sin ( t ) , α = 1 / 2 , λ = 3 / 2 , M = 8 , and T e = 0.05 . (b) Analytical solution (64) and numerical IIRFM solution (71) when f ( t ) = e t , α = 1 / 2 , λ = 3 / 2 , M = 8 , and T e = 0.05 .
Figure 10. (a) Analytical solution (61) and numerical IIRFM solution (62) when f ( t ) = sin ( t ) , α = 1 / 2 , λ = 3 / 2 , M = 8 , and T e = 0.05 . (b) Analytical solution (64) and numerical IIRFM solution (71) when f ( t ) = e t , α = 1 / 2 , λ = 3 / 2 , M = 8 , and T e = 0.05 .
Mathematics 12 02416 g010
Figure 11. IIRFM solution of Equation (3), for λ = α = 1 / 2 , f ( t ) = e t , T e = 5 × 10 3 , and M = 4 . (a) Changes in analytical and numerical solutions as a function of t. (b) Relative error ε r = 100 × ( u 1 2 , 1 2 , e t IIR , 4 / u 1 2 , 1 2 , e t 1 ) in %.
Figure 11. IIRFM solution of Equation (3), for λ = α = 1 / 2 , f ( t ) = e t , T e = 5 × 10 3 , and M = 4 . (a) Changes in analytical and numerical solutions as a function of t. (b) Relative error ε r = 100 × ( u 1 2 , 1 2 , e t IIR , 4 / u 1 2 , 1 2 , e t 1 ) in %.
Mathematics 12 02416 g011
Figure 12. IIRFM resolution of Equation (2) considering a logarithmic kernel, with λ = 1 , f ( t ) = t + 1 4 t 2 [ 3 2 ln ( t ) ] , T e = 10 2 , M = 5 , and an interpolation interval I s = [ s λ , 5 s λ ] . (a) Changes in the analytical and numerical solutions as a function of t. (b) Relative error ε r = 100 × ( u 1 , f IIR , 5 / u 1 , f 1 ) in %.
Figure 12. IIRFM resolution of Equation (2) considering a logarithmic kernel, with λ = 1 , f ( t ) = t + 1 4 t 2 [ 3 2 ln ( t ) ] , T e = 10 2 , M = 5 , and an interpolation interval I s = [ s λ , 5 s λ ] . (a) Changes in the analytical and numerical solutions as a function of t. (b) Relative error ε r = 100 × ( u 1 , f IIR , 5 / u 1 , f 1 ) in %.
Mathematics 12 02416 g012
Figure 13. IIRFM resolution of Equation (2), considering a logarithmic kernel, with λ = 1 , f ( t ) = 1 e t / τ + λ t + ( τ t ) ln ( t ) + λ τ e t / τ γ ln ( τ ) + i π + t / τ e x x d x , T e = 2.5 × 10 2 , M = 5 , τ = 1 / 2 , and an interpolation interval I s = [ s λ , 5 s λ ] . (a) Changes in the analytical and numerical solutions as a function of t. (b) Relative error ε r = 100 × ( u 1 , f IIR , 5 / u 1 , f 1 ) in %.
Figure 13. IIRFM resolution of Equation (2), considering a logarithmic kernel, with λ = 1 , f ( t ) = 1 e t / τ + λ t + ( τ t ) ln ( t ) + λ τ e t / τ γ ln ( τ ) + i π + t / τ e x x d x , T e = 2.5 × 10 2 , M = 5 , τ = 1 / 2 , and an interpolation interval I s = [ s λ , 5 s λ ] . (a) Changes in the analytical and numerical solutions as a function of t. (b) Relative error ε r = 100 × ( u 1 , f IIR , 5 / u 1 , f 1 ) in %.
Mathematics 12 02416 g013
Figure 14. Time changes in the temperature v ( t ) of a Joule-heated source point at constant intensity, calculated using the IIRFM (relations (93) and (94)).
Figure 14. Time changes in the temperature v ( t ) of a Joule-heated source point at constant intensity, calculated using the IIRFM (relations (93) and (94)).
Mathematics 12 02416 g014
Figure 15. Temperature of a semi-infinite solid ( x > 0 ) with constant Robin boundary conditions at x = 0 , due to a fluid flowing at constant velocity u and temperature T far from the solid, with constant Newton coefficient h. (a) Temperature profiles at different times, for T > T i . (b) Time evolution of the T ( 0 , t ) surface temperature when T = cste . The symbols correspond to the analytical solution (100) and the solid line to the numerical solution (102) obtained by the IIRFM method.
Figure 15. Temperature of a semi-infinite solid ( x > 0 ) with constant Robin boundary conditions at x = 0 , due to a fluid flowing at constant velocity u and temperature T far from the solid, with constant Newton coefficient h. (a) Temperature profiles at different times, for T > T i . (b) Time evolution of the T ( 0 , t ) surface temperature when T = cste . The symbols correspond to the analytical solution (100) and the solid line to the numerical solution (102) obtained by the IIRFM method.
Mathematics 12 02416 g015
Figure 16. Time changes in the surface temperature θ ( 0 , t ) of a semi-infinite solid, for sinusoidal variations θ ( t ) = θ ^ sin ( ω t ) of the fluid temperature far from the solid and for a constant flow velocity u .
Figure 16. Time changes in the surface temperature θ ( 0 , t ) of a semi-infinite solid, for sinusoidal variations θ ( t ) = θ ^ sin ( ω t ) of the fluid temperature far from the solid and for a constant flow velocity u .
Mathematics 12 02416 g016
Table 1. Influence of λ parameter on absolute error ε N H ( 5 | λ | ) = | u λ , t ( 5 | λ | ) v N H ( 5 | λ | ) | committed at t = 5 | λ | with the HPM-L method, when K ( t , t ) = 1 , f ( t ) = t and for different values of N H . Numbers in parentheses indicate multiplying powers of ten.
Table 1. Influence of λ parameter on absolute error ε N H ( 5 | λ | ) = | u λ , t ( 5 | λ | ) v N H ( 5 | λ | ) | committed at t = 5 | λ | with the HPM-L method, when K ( t , t ) = 1 , f ( t ) = t and for different values of N H . Numbers in parentheses indicate multiplying powers of ten.
N H =357152030
λ = 0.5 2.67144 ( 1 ) 1.28311 ( 2 ) 3.42505 ( 4 ) 3.66374 ( 12 ) 8.88178 ( 16 ) -
λ = 1.0 1.09080 ( + 1 ) 5.69965 ( + 1 ) 1.97941 ( + 1 ) 1.02417 ( 2 ) 1.20352 ( 5 ) 6.82121 ( 13 )
λ = 1.5 5.10447 ( + 4 ) 4.95956 ( + 4 ) 4.47047 ( + 4 ) 5.46097 ( + 3 ) 3.04954 ( + 2 ) 4.77714 ( 2 )
Table 2. Influence of λ on absolute error ε T e ( 5 | λ | ) = | u λ , t ( 5 | λ | ) u ˜ λ , t IIR ( 5 | λ | ) | at t k = 5 | λ | using IIRFM, when α = 0 , f ( t ) = t and different values of T e . Numbers in parentheses indicate multiplying powers of ten.
Table 2. Influence of λ on absolute error ε T e ( 5 | λ | ) = | u λ , t ( 5 | λ | ) u ˜ λ , t IIR ( 5 | λ | ) | at t k = 5 | λ | using IIRFM, when α = 0 , f ( t ) = t and different values of T e . Numbers in parentheses indicate multiplying powers of ten.
T e = 10 1 10 3 10 5 T e = 10 1 10 3 10 5
λ = 0.5 1.49258 ( 4 ) 1.49218 ( 8 ) 1.50691 ( 11 ) λ = 0.5 1.81881 ( 3 ) 1.8179 ( 7 ) 1.32737 ( 10 )
λ = 1.0 2.80584 ( 5 ) 2.80777 ( 9 ) 2.58013 ( 11 ) λ = 1.0 6.20611 ( 1 ) 6.18386 ( 5 ) 2.5262 ( 8 )
λ = 1.5 1.81606 ( 7 ) 1.81704 ( 11 ) 3.11284 ( 12 ) λ = 1.5 1.09635 ( + 3 ) 1.08113 ( 1 ) 8.04227 ( 6 )
Table 3. Coefficient values { a 8 , i , b 8 , i } obtained by Chebyshev rational interpolation on I s = [ 0 , 5 ] when α = 1 / 2 , λ = 3 / 2 , and M = 8 . Numbers in parentheses indicate multiplying powers of ten.
Table 3. Coefficient values { a 8 , i , b 8 , i } obtained by Chebyshev rational interpolation on I s = [ 0 , 5 ] when α = 1 / 2 , λ = 3 / 2 , and M = 8 . Numbers in parentheses indicate multiplying powers of ten.
i12345678
a 8 , i 1.1929 ( 3 ) 1.0196 ( 2 ) 4.3171 ( 2 ) 1.4132 ( 1 ) 4.2352 ( 1 ) 1.3087 5.0193 53.639
b 8 , i 2.5442 ( 2 ) 1.5704 ( 1 ) 5.1737 ( 1 ) 1.3474 3.2497 8.2163 26.753 259.97
Table 6. Absolute errors ε in the numerical resolution of (80) by the Laguerre approximation (Singha et al. [14]) and by the IIRFM (with I s = [ 0 , 15 ] and T e = 10 3 ), and for the same order of approximation M = 5 .
Table 6. Absolute errors ε in the numerical resolution of (80) by the Laguerre approximation (Singha et al. [14]) and by the IIRFM (with I s = [ 0 , 15 ] and T e = 10 3 ), and for the same order of approximation M = 5 .
Approximated SolutionsAbsolute Errors: ε × 10 6
t t Laguerre IIRFM Laguerre IIRFM
0.2 0.44721360 0.42222142 0.44717924 24,99234
0.4 0.63245553 0.60061467 0.63247236 31,84117
0.6 0.77459667 0.75640000 0.77454516 18,19752
0.8 0.89442719 0.89246899 0.89446429 195837
1.0 1.00000000 1.01148725 1.00009060 11,48791
2.0 1.41421356 1.42927171 1.41397446 15,058239
Table 7. Absolute errors ε in the numerical resolution of (82) by the Laguerre approximation (Singha et al. [14]) and by the IIRFM (with I s = [ 0 , 15 ] and T e = 10 3 ), and for the same order of approximation M = 5 .
Table 7. Absolute errors ε in the numerical resolution of (82) by the Laguerre approximation (Singha et al. [14]) and by the IIRFM (with I s = [ 0 , 15 ] and T e = 10 3 ), and for the same order of approximation M = 5 .
Approximated SolutionsAbsolute Errors: ε × 10 6
t t 1 / 4 Laguerre IIRFM Laguerre IIRFM
0.2 0.66874030 0.63419588 0.66876039 34,54420
0.4 0.79527073 0.75840545 0.79527721 38,6536
0.6 0.88011174 0.86140153 0.88004956 18,71062
0.8 0.94574161 0.94609498 0.94584299 353101
1.0 1.00000000 1.01516222 1.00013177 15,162132
2.0 1.18920712 1.20439817 1.18901343 15,191194
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Heyd, R. Numerical Solution of Linear Second-Kind Convolution Volterra Integral Equations Using the First-Order Recursive Filters Method. Mathematics 2024, 12, 2416. https://doi.org/10.3390/math12152416

AMA Style

Heyd R. Numerical Solution of Linear Second-Kind Convolution Volterra Integral Equations Using the First-Order Recursive Filters Method. Mathematics. 2024; 12(15):2416. https://doi.org/10.3390/math12152416

Chicago/Turabian Style

Heyd, Rodolphe. 2024. "Numerical Solution of Linear Second-Kind Convolution Volterra Integral Equations Using the First-Order Recursive Filters Method" Mathematics 12, no. 15: 2416. https://doi.org/10.3390/math12152416

APA Style

Heyd, R. (2024). Numerical Solution of Linear Second-Kind Convolution Volterra Integral Equations Using the First-Order Recursive Filters Method. Mathematics, 12(15), 2416. https://doi.org/10.3390/math12152416

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop