A Fitted Mesh Cubic Spline in Tension Method for Singularly Perturbed Problems with Two Parameters

A numerical treatment via a diﬀerence scheme constructed by the Crank–Nicolson scheme for the time derivative and cubic spline in tension for the spatial derivatives on a layer resolving nonuniform Bakhvalov-type mesh for a singularly perturbed unsteady-state initial-boundary-value problem with two small parameters is presented. Error analysis of the constructed scheme is discussed and shown to be parameter-uniformly convergent with second-order convergence. Numerical experimentation is taken to conﬁrm the theoretical ﬁndings.

Singularly perturbed problems (SPPs) widely appear in modeling physical problems such as fluid dynamics, heat and mass transfer in chemical engineering, quantum mechanics [7], elasticity, the theory of plates and shells, oil and gas reservoir simulation, and magnetic-hydrodynamic flow in general [2] and two-parameter singularly perturbed problems (TPSPPs) occur in chemical flow reactor theory [5] as well as in the case of boundary layers controlled by suction (or blowing) of some fluid [8] in particular.
Several articles dealt with the solution of an SPP and showed that it is revealed to have a sudden change in the boundary layer region and behaves normally in the outer region.Due to this multiscale character, the classical numerical methods cease to give satisfactory numerical results and are not competitive computationally.erefore, to overcome the shortcomings, researchers are forced to construct numerical methods that are flexible enough to provide a solution whose accuracy is independent of the perturbation parameter(s), and its singular nature can easily be captured.e sounding numerical methods on hand are the specially designed (fitted) ones.Just to list some of these, the authors of [9][10][11][12][13][14][15] used the fitted operator method, the authors of [16][17][18] implemented the fitted mesh method, Gowrisankar and Natesan [19] constructed a layer-adapted scheme through the equidistribution of a positive integrable monitor function, and Khan et al. [20] applied variable mesh in which more mesh points are in the layer region.On the other part, the authors of [21,22] used unfitted methods and the results yielded are not ε− uniform.
In this paper, we construct a numerical scheme for the parabolic convection-diffusion problem (1-2) through the fitted mesh method in which Crank-Nicolson discretization is used for the time variable and cubic spline in tension on layer-adapted meshes of Bakhvalov-type for the spatial variable discretization.e scheme constructed from this composition is in its first use for the problem under consideration.

Properties of the Continuous Problem
In this section, we consider some properties of the continuous problem through the maximum principle and bounds on the solution and on its derivatives that enable us to see a priori estimates for the solution u(x, t) and its derivatives.
is contradicts the given statement.erefore, we can conclude that the minimum of Π(x, t) is nonnegative.

□ Lemma 2. A solution u(x, t) of equation (1) satisfies the estimate
where ‖.‖ is the maximum norm.
Proof.For the proof, one can refer [23].

□
Lemma 3. e derivatives of the solution u(x, t) satisfy the following bound for all nonnegative integers i, j such that i + 3j ≤ 4.
where the constant C is independent of ε and μ and depends only on the bounded derivatives of the coefficients and the source term.

Discretization and Mesh Generation
Consider a uniform mesh on the time domain [0, T] with M mesh intervals spaced by k � T/M and given by Applying the Crank-Nicolson method for the time variable in equations ( 1) and (2), we obtain the following two-level semidiscrete scheme: where q(x) A characteristic equation corresponding to equation ( 7) that is helpful to describe the boundary layers occurring at the two end points of Ω is is defines two continuous functions Let Lemma 4. For a fixed number 0 < ρ < 1, a certain order δ, and the solution  u(x) of equation ( 3), the following holds: 2 International Journal of Mathematics and Mathematical Sciences where C is a constant independent of ε, μ.
Proof.For the proof, readers can refer [24].

□
Lemma 5. e error estimate in the temporal direction satisfies Proof. e proof is given in [11].
, and [1 − τ 2 , 1], such that the second subinterval is with equidistant N/2 mesh subintervals and the rest two are with gradually graded N/4 mesh subintervals.e two transition points τ 1 and τ 2 are chosen such that exp(− Since we are using nonuniform meshes in the spatial direction, we do not assume τ 1 and τ 2 to be 1/4.Furthermore, the relation η − 1 j log(η j ) ≤ CN − 1 , j � 1, 2 is given in [24].
For the two boundary layers, following the technique in [3], we have two generating functions , and ψ 2 (1) � 0. Analogous to [24], the mesh points in the spatial direction can be defined by is gives condensed and refined meshes in the layer regions and points where the mesh switches from fine to coarse and vice versa are chosen to be x (N/4)+1 � τ 1 and x (3N/4)+1 � 1 − τ 2 .e step sizes are related by However, from the mean value theorem, there exists

International Journal of Mathematics and Mathematical Sciences
For n � N/4 + 1, . . ., 3N/4, from equation ( 15), we have where m n and θ is a free parameter used to maintain consistency.e first derivatives of Q(x, t m ) from the right and left of x � x n are given, respectively, by en, from the continuity condition, we have where Application of the L'Hopital's rule gives n , from equation ( 7), we can have where, from [25], 4 International Journal of Mathematics and Mathematical Sciences Now, adopting the cubic spline in tension used in [25] yields Substituting equations ( 24) and ( 25) into equation ( 26) and writing it in compact form, we arrive at where n � 2, 3, . . ., N, m � 2, 3, . . ., M, International Journal of Mathematics and Mathematical Sciences (28) e system (27) together with the given conditions is uniquely solvable.Furthermore, the matrix associated with this system of equations satisfies the property of an M-matrix that is given in [26].

Convergence and Stability Analysis
Obtaining stability conditions for the Crank-Nicolson by using either the Fourier analysis or the matrix method is subjected to severe restrictions.is is worse for problems with variable coefficients like in our case.us, the maximum principle analysis can be taken as an alternative means [27].So, to establish the parameter uniform convergence of the proposed scheme, let us see the following lemmas.

Lemma 7. Assume that μh
en, under these assumptions, the matrix associated with equation ( 27) at each time level is an M-matrix.Hence, the operator L M,N ε,μ satisfies the following discrete maximum principle.

Lemma 8. (discrete maximum principle). Let the discrete function
Proof.To follow the proof by contradiction, let there exist a mesh point (l, m) for l ∈ 1, 2, . . .
and suppose Ψ m l ≤ 0. en, this indicates l ≠ 1, N + 1, and for l ∈ 2, 3, . . ., N { }, we obtain □ Lemma 9.For any discrete function Φ defined on the tensor product Ω N x × Ω M t such that Φ(x n , t m ) � 0 for all (x n , t m ) on the boundary of Ω N x × Ω M t , then for q m n ≥ q * > 0.
Proof. e proof follows from the application of the discrete maximum principle to a mesh function given by □ Lemma 10. e truncation error in the spatial direction, denoted by TE x , satisfies Proof.Let the truncation error in the spatial variable with fixed time be given by e Taylor series expansion of each term in space up to the third order derivative and collecting where upon simplification and restricting the expansion of the coefficients to their first term yield International Journal of Mathematics and Mathematical Sciences  17) and ( 23), it can be further simplified to As the operator L M,N ε,μ satisfies the discrete maximum principle on the tensor product D M t × D N x and the error is estimated, the current method is uniformly stable in the maximum norm.Applying the triangular inequality from equations ( 12) and (38), we can have the following theorem.

□
Theorem 1.Let u(x, t), (x, t) ∈ D be the solution of equations (1) and (2) and U(x n , t m ), (x n , t m ) ∈ D N x × D M t be the solution of equation (13).en, the error estimate for the fully discrete scheme is given by erefore, the presented method is convergent independent of the perturbation parameters and its rate of convergence is two.

Computational Results and Discussion
To verify the theoretical results experimentally, we solve the following examples.
Example 1.For the problem in [28], Example 2. Consider the following singular perturbation initial-boundary-value problem in [28], Example 3. Consider the following singular perturbation initial-boundary-value problem in [2], As the considered examples have no exact solution, we calculate the absolute maximum errors using the double mesh principle as follows: where U k h is an approximate solution obtained using k and h step sizes in the t and x directions, respectively, and U k/2 h/2 is an approximate solution obtained by bisecting the step sizes.As well, the corresponding rate of convergence is defined by e maximum absolute error and rate of convergence for Examples 1-3 are given in Tables 1-3.From these results, one can observe that the current method converges independently of the perturbation parameters and gives more accurate numerical results than the existing results available

Conclusions
e present study provides the numerical solutions for a singularly perturbed unsteady-state initial-boundary-value problem with two small parameters.e method comprises of developing a scheme through discretization of the time variable by the Crank-Nicolson method on a uniform mesh and the space variable by cubic spline in tension with Bakhvalov-type mesh.e method is flexible and successful
e proof in the third subinterval can be shown in a similar way as in the first subinterval.□ 3.2.Total Difference Scheme.Let Q(x, t m ) interpolate the solution U(x n , t m ) � U m n and pass through the points (x n , t m , U m n ) and (x n+1 , t m , U m n+1 ) such that

Table 1 :
Comparison of the maximum pointwise error and rate of convergence for Example 1 for μ � 10 − 6 and different values of ε.

Table 2 :
Comparison of the maximum pointwise error and rate of convergence for Example 2 for μ � 10 − 6 and different values of ε.

Table 3 :
Comparison of the maximum pointwise error and rate of convergence for Example 3 for μ � 10 − 7 and different values of ε.International Journal of Mathematics and Mathematical Sciences in the literature.eseresultsarealso in excellent agreement with the theoretical findings.esolutionprofilesillustrated in Figures1-3indicate the meshes in the layer regions are more condensed and refined than the outer layer.From the log-log plots in Figures4-6, we can conclude that our method is parameter uniform.