Numerical Solution of Direct and Inverse Problems for Time-Dependent Volterra Integro-Differential Equation Using Finite Integration Method with Shifted Chebyshev Polynomials

: In this article, the direct and inverse problems for the one-dimensional time-dependent Volterra integro-differential equation involving two integration terms of the unknown function (i.e., with respect to time and space) are considered. In order to acquire accurate numerical results, we apply the ﬁnite integration method based on shifted Chebyshev polynomials (FIM-SCP) to handle the spatial variable. These shifted Chebyshev polynomials are symmetric (either with respect to the point x = L 2 or the vertical line x = L 2 depending on their degree) over [ 0, L ] , and their zeros in the interval are distributed symmetrically. We use these zeros to construct the main tool of FIM-SCP: the Chebyshev integration matrix. The forward difference quotient is used to deal with the temporal variable. Then, we obtain efﬁcient numerical algorithms for solving both the direct and inverse problems. However, the ill-posedness of the inverse problem causes instability in the solution and, so, the Tikhonov regularization method is utilized to stabilize the solution. Furthermore, several direct and inverse numerical experiments are illustrated. Evidently, our proposed algorithms for both the direct and inverse problems give a highly accurate result with low computational cost, due to the small number of iterations and discretization.


Introduction
An integro-differential equation (IDE) is an equation which contains both derivatives and integrals of an unknown function. Several situations in the branches of science and engineering can be demonstrated by developing mathematical models which are often in the form of IDEs, such as in RLC circuit analysis, the activity of interacting inhibitory and excitatory neurons, the Wilson-Cowan model, and so on; see Reference [1] for more applications. In fact, many of these problems cannot be directly solved, because we may not know all necessary information or an incomplete system may be provided. This has led to the study of both direct and inverse problems for a certain type of one-dimensional IDE involving time, which is called the one-dimensional time-dependent Volterra IDE (TVIDE). Hence, in this study, we investigate the TVIDE of the following form v t (x, t) + Lv(x, t) = t 0 κ 1 (x, η)v(x, η)dη + x 0 κ 2 (ξ, t)v(ξ, t)dξ + F(x, t), x 0 κ 2 (ξ, t)v(ξ, t)dξ, while several studies in the literature have considered similar problems containing only one of these two terms.
The Volterra IDE containing only an integration term with respect to time arises in many applications, including the compression of poro-viscoelastic media, blow-up problems, analysis of space-time-dependent nuclear reactor dynamics, and so on; see Reference [2]. The existence, uniqueness, and asymptotic behavior of its solution have been discussed in Reference [3]. There are many authors who have studied the numerical solution of this type of problem by using techniques such as the finite element method [2], finite difference method (FDM) [4], collocation methods in polynomial spline [5], the implicit Runge-Kutta-Nyström method [6], the Legendre collocation method [7], and so on.
On the other hand, the Volterra IDE containing only an integration term with respect to space has also been studied in various areas, such as for the one-dimensional viscoelastic problem and one-dimensional heat flow in materials with memory [8], modeling heat/mass diffusion processes, biological species coexisting together with increasing and decreasing rates of growth, electromagnetism, and ocean circulation, among others [9]. Moreover, the existence and uniqueness for this type of Volterra IDE were shown in Reference [8]. Consequently, abundant numerical methods have appeared for finding solutions to this type of Volterra IDE using, for example, spline collocation method [10], collocation method with implicit Runge-Kutta method [11], decomposition method [12], and so on.
However, our problem deals with a Volterra IDE involving both temporal and spatial integrations. There have been no results in the literature regarding the existence and uniqueness of solutions to this type of problem. In this paper, we concentrate on providing a decent numerical procedure to find approximate solutions for both the direct and inverse problems of the proposed TVIDE (1).
Generally, it is well-known that the classification of problems involving differential equations was defined by Hadamard [13] in 1902. Mathematical problems involving differential equations are well-posed if the following conditions hold: existence, uniqueness, and stability. Otherwise, the problem is called ill-posed; this normally occurs in the inverse problem. Even though the initial and boundary conditions are prescribed, it is not sufficient to guarantee that our inverse problem (1) has unique solutions β(t) and v(x, t). Hence, additional conditions (e.g., the observation or measurement of data) need to be involved. In practice, there are many kinds of additional conditions; for example, a fixed point of the system, an average time of the system, or an integral of the system. After the additional conditions has been added as an auxiliary condition in our inverse problem (1), we can obtain the existence and uniqueness of β(t) and v(x, t). However, the additional condition may contains measurement or observation errors, which may cause the instability in the solutions; namely, a small perturbation in the input data can produce a considerable error, especially for β(t). Thus, some regularization techniques are required to overcome the ill-posedness and stabilize the solution.
There exist many schemes which are generally used to solve both direct and inverse problems of Volterra IDEs, such as the above-mentioned methods. However, those methods utilize the process of approximating differentiation. It is well-known that numerical differentiation is very sensitive to rounding errors, as its manipulation task involves division by a small step-size. On the other hand, the process of numerical integration involves multiplication by a small step-size and, so, it is very insensitive to rounding errors. In recent years, the finite integration method (FIM) has been developed to find approximate solutions of linear boundary value problems for partial differential equations (PDEs). The concept of FIM is to transform a given PDE into an equivalent integral equation, following which a numerical integration method, such as the trapezoid, Simpson, or Newton-Cotes methods (see References [14][15][16]), are applied. In 2018, Boonklurb et al. [17] modified the traditional FIM by using Chebyshev polynomials to solve one-and two-dimensional linear PDEs and obtained a more accurate result compared to the traditional FIMs and FDM. However, their technique [17] has not yet been utilized to overcome the direct and inverse problems of TVIDE, which are the major focuses of this work.
In this paper, we formulate numerical algorithms for solving the direct and inverse problems of TVIDE (1). We manipulate the idea of FIM in Reference [17] by using shifted Chebyshev polynomials, which we call the FIM with shifted Chebyshev polynomials (FIM-SCP), to deal with the spatial variable and use the forward difference quotient to estimate the time derivative. We further apply the Tikhonov regularization method to stabilize our ill-posed problem (1). The rest of the paper is organized as follows. In Section 2, the definition and some basic properties concerning the shifted Chebyshev polynomial are given to construct the shifted Chebyshev integration matrices. The Tikhonov regularization method is also presented in Section 2. In Section 3, we use the FIM-SCP and the forward difference quotient to devise efficient numerical algorithms to find approximate solutions to the direct and inverse problems of (1). Then, we implement our proposed algorithms through several examples, in order to demonstrate their efficiency compared with their analytical solutions. Furthermore, we also display the time convergence rate and CPU time (s) in Section 4. Finally, the conclusion and some directions for future work are given in Section 5.

Preliminaries
In this section, we introduce some necessary tools for solving the direct and inverse problems of TVIDE (1): the FIM-SCP and the Tikhonov regularization method.

Shifted Chebyshev Integration Matrices
We first introduce the definition and some basic properties of shifted Chebyshev polynomials [18], which are used to establish the first-and higher-order shifted Chebyshev integration matrices based on the idea of constructing integration matrices in Reference [17]. However, we slightly modify this idea by instead using a shifted Chebyshev expansion suitable for solving our problem (1) without domain transformation. We give the definition and properties as follows. Definition 1. The shifted Chebyshev polynomial of degree n ≥ 0 is defined by S n (x) = cos n arccos 2x Note that this shifted Chebyshev polynomial is symmetric, either with respect to the point x = L 2 or the vertical line x = L 2 over [0, L], depending on its degree. Next, we provide some important properties of the shifted Chebyshev polynomial, which we use to constructing the shifted Chebyshev integration matrix, as follows. Lemma 1. (i) For n ∈ N, the zeros of S n (x) are symmetrically distributed over [0, L] and given by (ii) For r ∈ N, the r th -order derivatives of S n (x) at the endpoint b ∈ {0, L} are d r dx r S n (x) (iii) For x ∈ [0, L], the single-layer integrations of shifted Chebyshev polynomial S n (x) arē (iv) Let {x k } n k=1 be a set of zeros of S n (x), the shifted Chebyshev matrix S is defined by Then, it has the multiplicative inverse S −1 = 1 n diag(1, 2, 2, ..., 2)S .
Next, we use the above definition and properties of shifted Chebyshev polynomials to construct the shifted Chebyshev integration matrices. First, let N be a positive integer and L be a positive real number. Define an approximate solution u(x) of a certain differential equation by a linear combination of shifted Chebyshev polynomials S n (x); that is, Let x k for k ∈ {1, 2, 3, ..., N} be the interpolated points which are meshed by the zeros of S N (x) defined in (2). Substituting each x k into (4), it can be expressed (in matrix form) as which is denoted by u = Sc. The unknown coefficient vector can be performed by c = S −1 u. Let us consider the single-layer integration of u(x) from 0 to x k , which is denoted by U (1) (x k ); we obtain for k ∈ {1, 2, 3, ..., N} or, in matrix form, We denote the above matrix by U (1) =Sc =SS −1 u := Au, where A =SS −1 := [a ki ] N×N is called the first-order shifted Chebyshev integration matrix for the FIM-SCP; that is, a ki u(x i ).
Next, consider the double-layer integration of u(x) from 0 to x k , which denoted by U (2) (x k ). We have for k ∈ {1, 2, 3, ..., N}. It can be written, in matrix form, as U (2) = A 2 u. Similarly, we can calculate the n-layer integration of u(x) from 0 to x k , which is denoted by U (n) (x k ). Then, we have for k ∈ {1, 2, 3, ..., N}, which can be expressed, in matrix form, as U (n) = A n u.

Tikhonov Regularization Method
In this section, we briefly present the idea of the Tikhonov regularization method [19], which is usually applied to stabilize ill-posed problems, such as our inverse problem. Normally, the considered inverse problem can be represented by the system of m linear equations with n unknowns, as where b is the vector in the right-hand side, which is perturbed by some noise , and x is the solution of the system (5) after perturbation. Tikhonov regularization replaces the inverse problem (5) by a minimization problem to obtain an efficiently approximate solution, which can be described as arg min where λ > 0 is a regularization parameter balancing the weighting between the two terms of the function and · is the standard Euclidean norm. To reformulate the above minimization problem (6), we obtain arg min Clearly, this is a linear least-square problem in x. Then, the above problem turns out to be the normal equation of the form To simplify the above equation, the solution x under the regularization parameter λ (denoted by x λ ) can be computed by We can see that the accuracy of x λ in (7) depends on the regularization parameter λ, which plays an important role in the calculation: A large regularization parameter may over-smoothen the solution, while a small regularization parameter may lose the ability to stabilize the solution. Therefore, a suitable choice of the regularization parameter λ is very significant for finding a stable approximate solution. There are many approaches for choosing a value of the parameter λ, such as the discrepancy principle criterion, the generalized cross-validation, the L-curve method, and so on. Nevertheless, the regularization parameter λ in this paper is chosen according to Morozov's discrepancy principle combined with Newton's method, which has been proposed in Reference [20]. We provide the procedure for calculating the optimal regularization parameter λ below, which can be carried out by the following steps: Step 1: Set n = 0 and give an initial regularization parameter λ 0 > 0.
Therefore, we receive the optimal regularization parameter λ, which is the terminal value λ n obtained from the above procedure. When the regularization parameter λ is fixed as the mentioned optimal value, we can directly obtain the corresponding regularized solution by (7).

Numerical Algorithms for Direct and Inverse Problems of TVIDE
In this section, we apply the FIM-SCP described in Section 2.1 to devise the numerical algorithms for solving both the direct and inverse TVIDE problems (1), in order to obtain accurate approximate results. Let u be an approximate solution of v in (1). Then, we have the following linear TVIDE over the domain Ω = (0, L) × (0, T]: subject to the initial condition and the boundary conditions for b ∈ {0, L} and r ∈ {0, 1, 2, ..., n − 1}, where t and x represent time and space variables, respectively. Additionally, κ 1 , κ 2 , F, φ, and ψ r are given continuous functions and L is the spatial linear differential operator of order n defined by L : is given and sufficiently smooth.

Procedure for Solving the Direct TVIDE Problem
First, we linearize (8) by uniformly discretizing the temporal domain into M subintervals with time step τ. Then, we specify (8) at a time t m = mτ for m ∈ N and use the first-order forward difference quotient to estimate the time derivative term u t . Next, we replace each x by x k for k ∈ {1, 2, 3, ..., N} as generated by the zeros of the shifted Chebyshev polynomial S N (x) defined in (2). Thus, we have The above equation can be written, in matrix form, as where each parameter in (12) can be defined as follows: Then, we consider the second integral term with respect to space by letting it be J m 2 (x k ) and using the idea of FIM-SCP (as described in Section 2.1) to approximate it. Then, we obtain for each x k ∈ {x 1 , x 2 , x 3 , ..., x N }. The above equation can be written, in matrix form, as where A =SS −1 is the shifted Chebyshev integration matrix defined in Section 2.1, Then, we apply the FIM-SCP (described in Section 2.1) to eliminate all spatial derivatives from (11) by taking the n-layer integral on both sides of (11), to obtain the following equation at the shifted Chebyshev node x k , as defined in (2), as  To simplify the n-layer integration of the spatial derivative terms of Lu m , by letting it be Q m (x k ) and using the technique of integration by parts, we have where d 1 , d 2 , d 3 , ..., d n are the arbitrary constants which emerge from the process of integration by parts. Then, we substitute each x k ∈ {x 1 , x 2 , x 3 , ..., x N } into the above equation and utilize the idea of FIM-SCP. Thus, we can express it, in matrix form, by Finally, we vary all points x k ∈ {x 1 , x 2 , x 3 , ..., x N } in (14) and rearrange them into matrix form by using the FIM-SCP with the derived matrix equations (12), (13), and (15); thus, we obtain or, factorizing the unknown solution u m explicitly, as Next, consider the given boundary conditions (10) at the endpoints b ∈ {0, L}. We can convert them into matrix form by using the linear combination of shifted Chebyshev polynomial (4) in term of the r th -order derivative of u at the iteration time t m and using (3). Then, we have for all r ∈ {0, 1, 2, ..., n − 1}. We can express the above equation, in matrix form, as which can be denoted by Bc m = Ψ m or BS −1 u m = Ψ m . Finally, we can construct the system of m th iterative linear equations from (16) and (17), which has N + n unknowns containing u m and d, as follows: where H m is the coefficient matrix of u m in (16) and E m 1 is the right-hand side column vector of (16). Consequently, the solution u m can be approximated by solving the system (18) starting from the given initial condition (9); that is, Note that, when we would like to find a numerical solution u(x, t) at any point x ∈ [0, L] for the terminal time T, we can calculate it by the following formula:

Procedure for Solving Inverse Problem of TVIDE
For the inverse problem in this paper, we specifically define the forcing term F(x, t) := β(t) f (x, t), where β(t) is a missing source function to be retrieved and f (x, t) is the given function. Thus, our considered time-dependent inverse TVIDE problem (1) becomes where u is an approximate solution of v and the other parameters are defined as in (8). The initial and boundary conditions of (19) are (9) and (10), which satisfy the compatibility conditions. Now, we remove all spatial derivatives from (19) and use the shifted Chebyshev integration matrix (as explained in Section 2.1). Then, we obtain the following matrix equation, based on the same process as in (16), as where .., f (x N , t m )] and the other parameters in (20) are as defined in Section 3.1. However, the occurrence of missing data is caused by the given conditions being insufficient to ensure a unique solution to our inverse problem. Hence, an additional condition or observed data needs to be involved. Thus, we use an additional condition, regarding the aggregated solution of the system, in the following form: where g(t) is the measured data at time t, which probably contains measurement errors. In order to illustrate the realistic phenomena of this problem, we assume that the measurement data of the aggregated solution g(t) involves some noise , which is denoted by g (t) (where g (t) − g(t) ≤ ) and define the noisy value by a random variable generated by the Gaussian normal distribution with mean µ = 0 and standard deviation σ = p|g(t)|, where p is the percentage of the noise to be input. Then, the additional condition (21) becomes Using the concept of FIM-SCP, the additional condition (22) at time t m can be written, in vector form, as where z = [S 0 (L),S 1 (L),S 2 (L), ...,S N−1 (L)] and eachS n (L) is as defined in Lemma 1(iii). Finally, we can establish the following system of m th iterative linear equations for the inverse TVIDE problem (19) by utilizing (20) and (23), which has N + n + 1 unknown variables including u m , d, and β m , as where H m is the coefficient matrix of u m defined in (20) and E m 2 is the right-hand side column vector of (20). Before seeking an approximate solution u m and source term β m , as we have mentioned, we must address that our inverse problem is ill-posed. When a noisy value is input into the system, it may cause a significant error. Hence, we need to stabilize the solution of (24) by employing the Tikhonov regularization method. We denote the linear system (24) by the simplified matrix equation as Applying the Tikhonov regularization method (6) in order to filter out the noise in the corresponding perturbed data, we can stabilize the numerical solution (25) by using (7). Thus, we have Finally, we can receive the optimal regularization parameter λ by using Morozov's discrepancy principle combined with Newton's method, as described in Section 2.2. Thus, we can directly obtain the corresponding regularized solution by (26).

Algorithms for Solving the Direct and Inverse TVIDE Problems
For computational convenience, we summarize the aforementioned procedures for finding approximate solutions to the direct (8) and inverse (19) TVIDE problems in Sections 3.1 and 3.2, respectively, as the numerical Algorithms 1 and 2, which are in the form of pseudocode. and F(x, t). Output: An approximate solution u(x, T). Update m = m + 1.

Algorithm 2 Numerical algorithm for solving the inverse TVIDE problem via FIM-SCP
Output: An approximate solution u(x, T) and the source terms β(t m ) at all discretized times. Set the measurement data g (t m ) = g(t m ) + , where ∼ N (0, p 2 |g(t m )| 2 ). Compute y λ n = (R R + λ n I) −1 R b . 10: Compute ∇y λ n = −(R R + λ n I) −1 y λ n .

17:
Find u m and β m by explicitly solving y λ using the matrix equation (7). 18: Update m = m + 1.

Numerical Experiments
In this section, we implement our devised numerical algorithms for solving the direct and inverse TVIDE problems through several examples, in order to demonstrate the efficiency and accuracy of the solutions obtained by proposed methods. Examples 1 and 2 are used to examine Algorithm 1 for the direct TVIDE problems (8). Examples 3 and 4 are inverse TVIDE problems (19), as solved by Algorithm 2. Additionally, time convergence rates and CPU times(s) for each example are presented to indicate the computational cost and time. The time convergence rate is defined by , where T is the terminal time, t m is a partitioned time contained in [0, T], u * (t m ) is the exact solution at time t m , u(t m ) is the numerical solution at time t m , and · ∞ is the l ∞ norm. Graphical solutions of each example are also depicted. Our numerical algorithms were implemented using the MatLab R2016a software, run on a Intel(R) Core(TM) i7-6700 CPU @ 3.40 GHz computer system. Example 1. Consider the following direct TVIDE problem, which consists of a second-order derivative with constant coefficient for x ∈ (0, 1) and t ∈ (0, T]: where subject to the homogeneous initial condition u(x, 0) = 0 for x ∈ [0, 1] and the Dirichlet boundary conditions u(0, t) = t and u(1, t) = te for t ∈ [0, T]. The analytical solution of this problem is u * (x, t) = te x .
In the numerical testing based on Algorithm 1, we first took the double-layer integral of both sides of (27) and transformed it into matrix form (16). Then, we obtained the approximate solutions u(x, T) for this problem (27) by applying the numerical Algorithm 1. The accuracy of our obtained approximate results was measured by the mean absolute error, which compared it to the analytical solution at different values of x ∈ {0.1, 0.3, 0.5, 0.7, 0.9} and the terminal time T = 1, as shown in Table 1. From Table 1, we observe that, when the partitioning number of the temporal domain M was fixed and nodal numbers N were increasingly varied, then the accuracy was significantly improved. Similarly, for a fixed nodal number N but various time partitioning numbers M, the accuracy results were also significantly improved. Moreover, the convergence rates with respect to the time in Algorithm 1 were estimated for various numbers of the time partition (M ∈ {5, 10, 15, 20, 25}) for the spatial points N = 10, as shown in Table 2. We can notice, from Table 2, that these time convergence rates for the ∞ norm indeed approached linear convergence for T ∈ {5, 10, 15}. The computational cost, in terms of CPU times (s), is also displayed in Table 2. Finally, a graph of our approximate solutions u(x, t) for different times t and the surface plot of the solution under the parameters N = 20, M = 20, and T = 1 are depicted in Figure 1.
We test the efficiency and accuracy of the proposed Algorithm 1 via the problem (28). First, we took a triple-layer integral on both sides of (28) and utilized the shifted Chebyshev integration matrix to transform it into matrix form (16). Next, we implemented Algorithm 1 to obtain numerical solutions u(x, T) for this problem (28). Table 3 shows the precision of our obtained approximate results at different values of x ∈ {0.1, 0.3, 0.5, 0.7, 0.9} and at the terminal time T = 1, through the mean absolute error. We can see that the accuracy was significantly improved according to an increase in the number of both the partitioning space and time domains. However, we observe that, in the case of fixed N, when M was increased, the mean absolute errors provide accurate results with a lower computational number M. Furthermore, the time convergence rates concerning the ∞ norm and CPU times (s) are demonstrated in Table 4  Example 3. Consider the following inverse TVIDE problem, which consists of a second-order derivative with constant coefficient and a continuous forcing function f (x, t) for x ∈ (0, 1) and t ∈ (0, T]: where subject to the initial condition u(x, 0) = e x for x ∈ [0, 1] and the boundary conditions u(0, t) = t + 1 and u(1, t) = t + e for t ∈ [0, T]. The additional condition, in terms of the aggregated solution of the system, is g(t) = t + e − 1. The analytical solutions of this problem are u * (x, t) = t + e x and β * (t) = e −2t .
Implementing the numerical Algorithm 2 by taking the double-layer integral of both sides of (29) and transforming it into matrix form (24), we obtained the approximate solutions u(x, 1) and β(t) for this problem (29). As the additional condition was measurement data, there may be an error in the measurement. Therefore, we perturbed the additional condition g(t) with a percentage p of the noise (p ∈ {0%, 1%, 3%, 5%}). In Table 5, we show the accuracy of the solutions u(x, 1) and β(t), in terms of the mean absolute error, respectively, denoted by E u = 1  Table 5, we can observe that the optimal regularization parameters λ were close to zero and the mean absolute errors for both E u and E β significantly increased with an increasing percentage p of the perturbation. Furthermore, we used the regularization parameter λ = 0 to explore the rates of convergence with respect to the ∞ norm and CPU times (s) for M = N ∈ {5, 10, 15, 20} with the final times T ∈ {1, 2, 3} as shown in Table 6. The graphical solutions of the perturbed functions u(x, 1) and β(t) for p ∈ {1%, 3%, 5%} are depicted in Figure 3. Table 5. Mean absolute errors of u(x, 1) and β(t) for optimal regularization parameter λ of Example 3.   Example 4. Consider the following inverse TVIDE problem, which consists of a second-order derivative with variable coefficient and the piecewise forcing function f (x, t) for x ∈ (0, 1) and t ∈ (0, T]: u t + u xx + u x − cos(xt)u = t 0 2 sin(x)u(x, η)dη − x 0 3t cos(ξ)u(ξ, t)dξ + β(t) f (x, t), where f (x, t) =       

Conclusions and Discussion
In this paper, we utilized FIM-SCP combined with the forward difference quotient to create efficient and accurate numerical algorithms for solving the considered direct and inverse TVIDE problems. According to the numerical examples in Section 4, we have demonstrated the performance of our proposed Algorithm 1 for seeking the approximate solutions of direct TVIDE problems in Examples 1 and 2. We can see that, for Example 1-which involved a second-order derivative with constant coefficients-Algorithm 1 provided an accurate result. Furthermore, for a problem involving a higher-order derivative with variable coefficients, it still provided high accuracy, in terms of solutions, as demonstrated in Example 2. Moreover, we handled inverse TVIDE problems using Algorithm 2, the effectiveness of which was illustrated in Examples 3 and 4. We used the Tikhonov regularization method to deal with the instability of the inverse problem; it can be seen that, in the examples, the regularization parameter λ was close to zero. Algorithm 2 could handle both continuous and piecewise-defined forcing terms with high accuracy, as demonstrated in Examples 3 and 4. Furthermore, when we perturbed the problems by adding noisy values, our Algorithm 2 still overcame the noise and provided approximate results that approached the analytical solutions. We further notice that our presented methods provide high accuracy, even when using only a small number of nodal points. Evidently, when we decrease the time step, they will furnish more accurate results. The rates of convergence with respect to time (based on the ∞ norm) of our methods were observed to be linear.
Finally, we also depicted the computational times for each example. However, we realize that there exist no theoretical error analysis results for the proposed numerical algorithms. Thus, our future research will study the error analysis, in order to find theories for order of accuracy and rate of convergence for our method. Another interesting direction for our future work is to extend our techniques to solve other types of IDEs and non-linear IDEs.