To avoid singular generalized Jacobian matrix and further accelerate the convergence of the generalized Newton (GN) iteration method for solving generalized absolute value equations Ax - B|x| = b, in this paper we propose a new relaxed generalized Newton (RGN) iteration method by introducing a relaxation iteration parameter. The new RGN iteration method involves the well-known GN iteration method and the Picard iteration method as special cases. Theoretical analyses show that the RGN iteration method is well defined and globally linearly convergent under suitable conditions. In addition, a specific sufficient condition is studied when the coefficient matrix A is symmetric positive definite. Finally, two numerical experiments arising from the linear complementarity problems are used to illustrate the effectiveness of the new RGN iteration method.
1.
Introduction
Delay differential equations model problems in several domains, including biosciences, material science, and medicine[1,2]. Differential-difference equations are differential equations in which the system's evolution depends on both the system's historical context and its present state. A differential equation that consists of at least one shift term and whose highest-order derivative is multiplied by a small perturbation parameter is known as a singularly perturbed differential-difference equation. Singularly perturbed differential difference equations (SPDDEs) generally lead to solutions exhibiting boundary layers, and as the perturbation parameter goes to zero, the smoothness of the solution deteriorates.
The initial developments on the asymptotic analysis of singularly perturbed differential-difference equations have emerged from the articles by Lange and Miura [3,4]. The numerical explorations in this field can be found in [5], where the authors Kadalbajoo and Sharma have presented a ϵ-uniform numerical scheme comprising of a standard upwind finite difference operator on a fitted piecewise uniform mesh for a class of boundary value problems of singularly perturbed differential-difference equations with small shifts. In [6], the authors Kadalbajoo et al. developed a fitted operator and a fitted mesh finite difference method for a class of singularly perturbed difference-difference equations, the fitted mesh method being ϵ-uniform convergent of second order. Kadalbajoo and Ramesh [7] proposed a hybrid numerical method on Shishkin mesh for solving a singularly perturbed delay differential equation, wherein the solutions were compared with those obtained by using a simple upwind scheme and a midpoint upwind scheme. Sirisha et al. [8] presented a mixed finite difference method for singularly perturbed differential-difference equations with mixed shifts via domain decomposition on a constant mesh. Woldaregay and Duressa [9] presented a numerical scheme for singularly perturbed differential-difference equations with mixed small shifts by using the exponentially fitted operator finite difference method, and in [10] they applied an exponentially finite difference method to a singularly perturbed boundary value problem. Ranjan and Prasad [11] used an exponentially fitted three-term finite difference technique to approximate the solution of a singularly perturbed differential equation with small shifts. Kumar and Kadalbajoo [12] constructed a piecewise uniform mesh for solving singularly perturbed differential-difference equations with small shifts.
Jain[13] introduced the spline function approximation and shown that they are the consistency relations for the fundamental equations in discrete mechanics. Kadalbajoo and Bawa [14] proposed a variable mesh difference scheme for singularly perturbed boundary value problems using splines, and the method is shown to be quadratically convergent. Kadalbajoo and Patidar [15,16] applied spline techniques such as spline in compression and spline in tension to singularly perturbed two-point boundary value problems. Aziz and Khan [17] studied a spline method for second-order singularly perturbed boundary value problems, and the convergence of the method is shown to be dependent on the choice of the parameters. Mohanty and Jha [18] applied the variable mesh method using spline in compression for singularly perturbed two-point singular boundary value problems. Mohanty and Arora [19] applied tension spline on a non-uniform mesh for singularly perturbed two-point singular boundary value problems with significant first derivatives. Chakravarthy et al.[20,21] presented the numerical solution using spline in compression and spline in tension on a uniform mesh for second-order singularly perturbed delay differential equations. Ravi Kanth and Murali [22] presented a numerical technique for solving singularly perturbed nonlinear delay differential equations by the method of Spline in compression. The quasilinearization technique is applied in converting the nonlinear equation into a sequence of linear equations.
The aforementioned publications illustrate the implementation of spline methods for various types of singularly perturbed boundary value problems and motivate us in exploring the feasibility of constructing non-polynomial spline methods for solving singularly perturbed boundary value problems with mixed shifts. In the present paper, we constructed numerical methods using spline in compression, spline in tension and adaptive spline on a fitted mesh, and a comparative study is performed on the results. This method ensures a consistent level of accuracy regardless of the perturbation parameter and achieves reliable convergence.
The content of this paper is organized as follows: In Section 2, we introduce the problem under consideration and some properties of the solution. In Section 3, we discuss the mesh construction strategy for the problem. Section 4 is devoted to the proposed methods for the problem. In Section 5, we discuss the convergence of the proposed methods. In Section 6, we present the numerical results for test problems, and finally the conclusions follow in Section 7.
2.
Statement of the problem
Consider the general boundary value problem (BVP) for SPDDE containing both types of shifts:
for x∈(0,1),0<ϵ≪1, subject to the interval and boundary conditions
where 0<ϵ≪1 is the perturbation parameter, a(x),α(x),ω(x),β(x),f(x), ϕ(x), and χ(x) are smooth functions, and δ is the delay (negative shift) parameter and η is the advance (positive shift) parameter. As δ,η<ϵ, for a(x)−δα(x)+ηβ(x)>0 and α(x)+ω(x)+β(x)<0 ∀ x∈[0,1], the solution exhibits a boundary layer near x=0, while for a(x)−δα(x)+ηβ(x)<0 and α(x)+ω(x)+β(x)<0, the solution exhibits a boundary layer near x=1. Here we assume that α(x)≤M1, β(x)≤M2, and α(x)+ω(x)+β(x)≤−M<0. The function y(x) being continuous in [0,1] and differentiable in (0,1), satsifying (2.1) and (2.2), provides a smooth solution for (2.1) and (2.2).
Since the solution of y(x) of (2.1) and (2.2) is sufficiently differentiable, we expand the terms y(x−δ) and y(x−η) using Taylor series to obtain:
On substituting Eq (2.3) in Eq (2.1), the modified form of the Eq (2.1) is
subject to the conditions
where Υ(x)≈y(x), μ(x)=ϵ+α(x)δ22+β(x)η22, p(x)=a(x)−δα(x)+ηβ(x), q(x)=α(x)+ω(x)+β(x) and r(x)=f(x).
2.1. Some properties of the solution
We show that the operator L follows the minimum principle for the continuous problem Eq (2.4):
Lemma 1. Let Υ(x) be a smooth function, with Υ(0)≥0,Υ(1)≥0, then for x∈[0,1], Υ(x)≥0, whenever L(Υ(x))≤0 for x∈(0,1).
Proof. For proof of this lemma, the reader can refer to [5]. □
The bound for the solution of the continuous problem (2.4) is given in the following lemma:
Lemma 2. Let Υ(x) be the solution of (2.4) and (2.5), then, ‖Υ‖≤M−1‖r‖+max(|ϕ0|,|χ1|), ‖.‖ being the l∞ norm ‖Υ‖=maxs∈[0,1]|Υ(s)|.
Proof. For proof of this lemma, the reader can refer to [5]. □
Lemma 1 guarantees the uniqueness of the solutions of (2.4) and (2.5), and the existence of the solution is guaranteed as the given problem is linear. Also, the boundedness to the solution of the problem is implied by Lemma 2. Also, the bounds for the solutions of (2.4) and (2.5) and their derivatives are given in the following lemma.
Lemma 3. Considering Υ(x) to be the solution of (2.4) and (2.5), we have ||Υ(k)||≤2kC(2ϵ+δ2M1+η2M2)−k for k=1,2,3.
Proof. For proof of this lemma, the reader can refer to [5]. □
Theorem 1. Let Υ(x) be the solution of (2.4) and (2.5), and let Υ(x)=Υr(x)+Υs(x), where the regular component Υr(x) satisfies
and the singular component Υs(x)satisfies
where 0≤k≤3.
Proof. For proof of this theorem, the reader can refer to [5]. □
3.
Mesh construction
In this section, we discuss the mesh generation for the numerical solution of the singularly perturbed BVP (2.4) and (2.5).
The case when the boundary layer occurs at the left end of the domain D=[0,1], it is divided into two subdomains D1 and D2 such that D=D1∪D2=[0,τ]∪[τ,1], where τ is the transition parameter that is closer to x=0, and is defined by
where N is the number of mesh points in the domain D=[0,1] and τ0≥1|M|. It is clear that τ=12, the mesh is uniform; otherwise, the mesh condenses near the left boundary. It is assumed that N=2m, where m≥2 is an integer, which guarantees that there is at least one point in the boundary layer region. So, we consider equal number of mesh points in each subdomain and uniform partition over each subdomain with mesh points xi, as defined by
where h1=2τN and h2=2(1−τ)N on the domains D1 and D2 respectively.
Similarly, in the case when the boundary layer occurs at the right end of the solution domain D, we divide into subdomains D∗1 and D∗2 such that D=D∗1∪D∗2=[0,1−τ]∪[1−τ,1], where τ is so-called the transition parameter and is located near the point x=1. We consider equal number of grid points in each subdomain and uniform partition over each subdomain with grid points xi, as defined by
Now we show that L satisfies the discrete minimum principle:
Lemma 4. If the mesh function Υ(xi) satisfying Υ(x0)≥0,Υ(xN)≥0, then Υ(xi)≥0, 0≤xi≤1, for L(Υ(xi))≤0, 0<xi<1.
Proof. Let 0≤ˉzk≤1 be such that Υ(ˉzk)=minx∈[0,1]Υ(xi), and assuming that Υ(ˉzk)<0, clearly ˉzk∉{0,1}. Hence Υ′(ˉzk)=0 and Υ′′(ˉzk)≥0.
Now we have L(Υ(ˉzk))=μ(ˉzk)Υ′′(ˉzk)+p(xi)Υ′(ˉzk))+q(xi)Υ(ˉzk)>0, which is a contradiction to our assumption that Υ(ˉzk)<0. Therefore, Υ(ˉzk)≥0 and hence Υ(xk)≥0∀xi∈[0,1]. □
Lemma 5. Let Υi be any mesh function such that Υ0=ΥN=0. Then, for all 0≤i≤N, ‖Υj‖≤M−1max1≤j≤N−1|L(Υj)|.
Proof. Let us introduce two mesh functions ^υi± defined by
and for 1≤i≤N−1
Therefore, by the discrete minimum principle, we have ˆυ±i≥0 ∀i,0≤i≤N, which gives the required estimate. □
Lemma 6.
for each i.
Proof.
The above inequality is true for each j. Now we multiply these inequalities for j=i+1,...,N, and we obtain
Hence the result. □
Lemma 7. For i=0,1,...,N, we set
then for i=0,1,...,N−1, we have
Proof. It is easy to verify that
Now
from which the result follows. □
Lemma 8. There exists a constant C such that
for N/2≤i≤N.
Proof. suppose N/2≤i≤N. By [23]
It is easy to verify that N16(i−1/N)N−1lnN/(1+4N−1lnN) is bounded for any N≥2 from which the result follows. □
4.
Numerical methods
In this section, we present non-polynomial spline methods for solving the boundary value problems (2.4) and (2.5).
4.1. Spline in compression
The spline in compression SΔ(x) satisfies in [xi−1,xi] the differential equation
where SΔ(xi)=Υi,ψ>0,hi=xi−xi−1.
Solving (4.1) as a second-order differential equation, we obtain
Applying the interpolating conditions at xi−1 and xi; SΔ(xi−1)=Υi−1, SΔ(xi)=Υi,S′′Δ(xi)=Mi, and setting λi=√ψhi, we obtain the interpolating constants A and B and hence
Differentiating (4.3) and taking x→xi, we obtain
Considering the interval (xi,xi+1) and similarly we obtain
Equating the left and right hand dertivatives at xi, we obtain
This leads to a tridiagonal system
where ¯λ=1λ2i[λisinλi−1] and ¯¯λ=1λ2i[1−λicotλi].
The consistency relation for (4.5) may be expressed as λi2=tan(λi2), whose smallest positive root is λi≈8.986818916, which leads to the equation ¯λ+¯¯λ=12.
To obtain an approximation for Υ′i and Υ′′i, we use the Taylor series approximation for Υ about xi as:
Solving (4.6) and (4.7) for Υ′i and Υ′′i, we will obtain
Using the above approximations (4.8) and (4.9) in Υ′i+1=Υ′i+hi+1Υ′′i and Υ′i−1=Υ′i−hiΥ′′i, we obtain
We write Eq (2.4) as
Now we rewrite Eq (4.5) as
Substituting (4.12) in (4.13) and using the approximations (4.8), (4.10), and (4.11), we obtain the following tridiagonal scheme:
where
Using the Thomas algorithm, we can solve the above tri-diagonal scheme (4.14) subject to the boundary conditions (2.5).
4.2. Spline in tension
The spline in tension SΔ(x) in [xi−1,xi] satisfies the differential equation
where ψ>0 is a tension factor, SΔ(xi)=Υi, S′Δ(xi)=mi, S′′Δ(xi)=Mi, hi=xi−xi−1.
Solving (4.15) as a second-order differential equation, we obtain
Applying the interpolating conditions at xi−1 and xi and setting Λi=√ψhi, we obtain the interpolating constants A and B, and hence
Differentiating (4.17) and taking x→xi, we obtain
Considering the interval (xi,xi+1) and similarly, we obtain
Equating the left and right hand dertivatives at xi, we have
This leads to a tridiagonal system
where Λ1=1Λ2i[1−ΛisinhΛi] and Λ2=1Λ2i[ΛicothΛi−1].
We rewrite the Eq (4.19) as
Substituting (4.12) in (4.20) and using the approximations (4.8), (4.10), and (4.11), we obtain the following tridiagonal linear system:
where
Using the Thomas algorithm, we can solve the above tri-diagonal scheme (4.21) subject to the boundary conditions (2.5).
4.3. Adaptive spline
The function SΔ(x), which we call adaptive spline, satisfies the following differential equation:
Solving (4.22) and using interpolatory constraints SΔ(xi−1)=Υi−1, SΔ(xi)=Υi, we obtain
where νi=ψhi2Θ, zi=x−xi−1hi and Θ,ψ are constants.
The function SΔ(x) on the interval [xi,xi+1] is obtained by replacing i with i+1 (4.23).
Applying the conditions of continuity to the first or second derivative of SΔ(x) at xi, we obtain the following relationship:
which simplifies to the following form of tridiagonal system:
Some additional relations for the adaptive spline are listed as follows:
(i) mi−1=−hi(A1Mi−1+A2Mi)+Υi−Υi−1hi
(ⅱ) mi=hi(A3Mi−1+A4Mi)+Υi−Υi−1hi
(ⅱ) Mi−1=2νiςihi[−(A4mi−1+A2mi)+B1(Υi−Υi−1hi)]
(ⅳ) Mi=2νiςihi[(A3mi−1+A4mi)+B2(Υi−Υi−1hi)]
where ςi=νicothνi−12νi,
A1=14(1+2ςi)+ςi2νi, A2=14(1−2ςi)−ςi2νi, A3=14(1+2ςi)−ςi2νi and
A4=14(1−2ςi)+ςi2νi, B1=12(1−2ςi), B2=−12(1+2ςi).
In the limiting case, when νi→0, we have
ςi=0, ςiνi=16, A1=13, A2=16, A3=16, A4=13, B1=12, B2=−12
and the spline function (4.23) reduces to the cubic spline.
By introducing the parameter μi, we rewrite Eq (4.27) as
Substituting (4.12) in (4.28) and using the approximations (4.8), (4.10), and (4.11), we obtain the tridiagonal linear system of the form:
where
Using the Thomas algorithm, we can solve the above tri-diagonal scheme (4.29) subject to the boundary conditions (2.5).
5.
Convergence analysis
Here we perform the convergence analysis for the scheme described in Section 4.1.
Writing the tri-diagonal system Eq (4.14) in matrix-vector form, we obtain
in which A=[mi,j],1≤i,j≤N−1, is a tri-diagonal matrix of order N−1, with
and C=(di) is a column vector with di=h2i(¯λri−1+2¯¯λri+¯λri+1) with i=1,2,⋯,N−1 with T(hi)=O(h3i).
We also have
where (ˉΥ)=(¯Υ0,¯Υ1,⋯,¯ΥN)T and T(hi)=(T0(hi),T1(hi),⋯,TN(hi))T denote the actual solution and the local truncation error, respectively.
From Eqs (5.1) and (5.2), we obtain
Thus the error equation is
where E=ˉΥ−Υ=(eo,e1,e2,⋯,eN)T. Let |p(x)|≤c1, |q(x)|≤c2 and [mi,j] is the (i,j)th element of the matrix A. Then we have
For sufficiently small hi, we have
Hence the matrix is irreducible [24].
Let the ith row elements' sum of matrix A be Si, then we have
Let c1∗=min|p(x)|, c∗1=max|p(x)|, c2∗=min|q(x)|, c∗2=max|q(x)| and h=N−1maxi=1{hi,hi+1} so that 0<c1∗≤c1≤c∗1,0<c2∗≤c2≤c∗2.
Then for a given h, the matrix A is irreducible and monotone ([24,25]).
From (5.3), we have
At the end points i=0 and N, the above inequality holds, and for 1≤i≤N−1, we have
where ˆζ∈(xi−1,xi).
Since the mesh is piecewise uniform with step difference h, from (5.5) and (5.6), we obtain
Also, we have from [26]
as 0<h<1.
Using (5.7) and (5.8) in (5.5), we obtain
For a left layer problem, let the fine mesh points for the inside layer region be x1,⋯,xN/2, and the coarse mesh points in the outer region be xN/2+1,⋯,xN−1. Further Υr(x) and Υs(x) are the regular and the singular components of the numerical solution. Then, using (5.2) along with (3.1) and h1, h2 into (5.9) and Theorem 1, we obtain
Similarily,
The above results (5.10) and (5.11) can be concluded as
Theorem 2. a(x), \alpha(x), \omega(x), \beta(x), f(x) , \phi(x) , and \chi(x) be sufficiently smooth functions so that \Upsilon(x)\in C^3 [0, 1] . Let \Upsilon_i, \ i = 0(1)N be the approximate solution of (2.4), obtained using the fitted mesh finite difference method (4.14) with the conditions (2.5). Then, there is a constant \mathcal{M} independent of \epsilon and the mesh size such that
6.
Numerical results
To check the efficiency of the methods described in Sections 4.1–4.3, four test problems of SPDDEs are solved, of which two problems are of left end boundary layer type and the other two are right layer problems.
The double mesh principle is used for finding the maximum absolute errors, which is given by the formula:
and the numerical rate of convergence for the considered problems is calculated by the following formula:
The numerical techniques outlined in Sections 4.1–4.3 are applied to the test problems, and the maximum absolute errors and the numerical rate of convergence are evaluated. The numerical results are tabulated for a spectrum of values of \delta and \eta , smaller than \epsilon for all the test problems. The findings are displayed in Tables 1–8. Also, the \epsilon -uniform maximum absolute errors E_N for various values of the mesh parameter N and for \epsilon \in \{2^0, 2^{-1}, \cdots 2^{-20}\} is compared for each method described in Section 4 in Table 9. The numerical work illustrates the efficiency of the methods and is also consistent with those in literature. Graphs illustrating the influence of the shift parameters on the solution of the problem are depicted in Figures 1–8. The relationship between the error E_N and the number of mesh points N for the considered examples is plotted in Figures 9–12. These plots illustrate the efficiency of the methods presented in Sections 4.1–4.3.
Example 1. \epsilon y^{\prime\prime}(x)+ y^{\prime}(x)-2y(x-\delta)+y(x)-y(x+\eta) = -1 , y(x) = 1;-\delta \leq x\leq 0, y(x) = 1, 1\leq x \leq 1+\eta .
Example 2. \epsilon y^{\prime\prime}(x)+2.5y^{\prime}(x)-2e^{x}y(x-\delta)-y(x)-xy(x+\eta) = 1 , y(x) = 1;-\delta \leq x\leq 0, y(x) = 1, 1\leq x \leq 1+\eta .
Example 3. \epsilon y^{\prime\prime}(x)- y^{\prime}(x)-2y(x-\delta)+y(x)-2y(x+\eta) = 0 , y(x) = 1;-\delta \leq x\leq 0, y(x) = -1, 1\leq x \leq 1+\eta .
Example 4. \epsilon y^{\prime\prime}(x)-(1+e^{x^2})y^{\prime}(x)-xy(x-\delta)+x^2y(x)-(1-e^{-x})y(x+\eta) = 1 , y(x) = 1;-\delta \leq x\leq 0, y(x) = -1, 1\leq x \leq 1+\eta .
7.
Discussion and conclusions
In this paper, we proposed fitted mesh numerical methods for solving singularly perturbed boundary value problems of second-order ordinary differential equations with mixed shifts. The shifts that are smaller than the perturbation parameter are approximated using Taylor series and non-polynomial splines, namely, the spline in compression, the spline in tension, and the adaptive spline are applied to the Shishkin mesh. The methods presented are analyzed for convergence and shown to be first-order convergent. Numerical computations are carried out on two test problems that exhibit layer behavior on the left of the underlying interval and two right-layer problems. The maximum absolute errors and rates of convergence are tabulated, which show the first-order uniform rate of convergence. Graphs are plotted for the test problems for different values of the perturbation and shift parameters. From the figures, the effect of the shifts on the boundary layer behavior of the solution of the problems can be observed. As the shifts increase in magnitude, the thickness of the layer decreases for the left-layer problems, while for the right-layer problems, it increases. The methods presented in this paper have been found to be almost equally efficient in achieving \epsilon -uniform convergence and the numerical rate of convergence. Hence, it can be concluded that the presented methods provide considerable advantage for solving singularly perturbed linear second-order boundary value problems with mixed shifts.
Author contributions
Both the authors contributed equally to this work.
Acknowledgments
The authors wish to thank the National Board for Higher Mathematics, Department of Atomic Energy, Government of India, for their financial support under the project No. 02011/8/2021 NBHM(R.P)/R&D Ⅱ/7224, dated 24.06.2021.
Conflict of interest
The authors declare that they have no conflict of interest, relevant to the content of this article.