Singular perturbation boundary and interior layers problems with multiple turning points

In the study of singularly perturbed boundary problems with turning points, the solution undergoes sharp changes near these points and exhibits various interior phenomena. We employ the matching asymptotic expansion method to analyze and solve a singularly perturbed boundary and interior layers problem with multiple turning points, resulting in a composite expansion that ﬁts well with the numerical solution. The solution demonstrates a strong association with special functions, which is veriﬁed by the theory of diﬀerential inequalities


Introduction
Many phenomena in biology, chemistry, engineering, physics, etc. can be described by boundary value problems associated with various types of differential equations or systems.Boundary value problems often apply in mathematical physics, such as wave equations, Laplace equations [1], etc., where reaction-diffusion systems [2][3][4] can be applied in biology, geology, physics (neutron diffusion theory), and ecology.In biology the moving boundary model for oxygen diffusion in a sick cell explores the possibility of numerical approximations to different problems, which not only provide many numerical experiments, but also comparisons with analytical solutions [5,6].
For some boundary value problems, the singular perturbation method is also an effective approximation tool to solve its asymptotic solution, which is a seminumerical and semianalytical method, such as the boundary layer function method [7], matched asymptotic expansion method, multiscale method, and WKB method [8].These methods can not only serve as a tools, but also provide valuable insights into studying general singular perturbation boundary value problems.
The singular perturbation boundary value problem is a class of boundary value problems with small parameters before the highest derivative.The most significant property of the singularly perturbed boundary value problems is the reduction of the order of the equation as the perturbation parameter ε approaches 0. Singular perturbation problems with turning points arise as mathematical models [9] for various physical phenomena.Many problems in quantum physics belong to them, such as the famous Schrödinger equation.
The turning point is the point that changes the stability of the problem.As the original problem passes through the turning point, the properties of the function will change.In classical mechanics, when the particle incident on these points, all the kinetic energy is converted into potential energy, and the particle begins to turn and reflect back [10].
The turning point theory is a branch of the asymptotic theory of singular perturbation problems.The "turning point" is an exception of this theory because its precise definition is not an ordinary thing in general.Because of the particularity of the turning point, its analysis requires a complete understanding of the asymptotic properties of the solutions of differential equations.Turning points in singular perturbation have been a difficult problem, so that discussion of these problems was avoided in the early days of the study [11].A word of explanation is perhaps appropriate.For a singularly perturbed boundary value problem εy = f (x, y, y ), in the case of a linear function f , a turning point is defined as a point at which the coefficient of y vanishes.We also could extend this definition to nonlinear functions f by requiring that the partial derivative ∂f ∂y = f y vanishes at the turning point in a rather large domain of variation of y and y .
According to the location of turning point, it can be divided into the internal turning point problem, boundary turning point problem, and multiturning point problem.Since the stability of the problem changes around the turning point, the solutions often change dramatically near these turning points.In the study of singularly perturbed boundary value problems with turning points, the solution often changes sharply near the turning point and then presents some interior layer phenomena [12,13], such as shock layer [14,15], nonmonotone transition layer, spike layer [16], etc.
There is another kind of situation that is different from the above cases, because the rapid changes will be in the slope, and not in the value of the solution itself.This is usually called the corner layer phenomenon.One of the most important applications for corner layer problems is the singularly perturbed problem that combines combustion theory and reaction-diffusion equations.
There are two main methods for solving singular perturbation problems, asymptotic and numerical methods.Asymptotic methods can help us understand the qualitative behavior of the problem, whereas numerical methods perform quantitative analysis.The most suitable method for boundary layer and interior layer problems is the matching asymptotic expansion method [17], including the Prandtl matching principle, the Vandyke matching principle, and the intermediate variable matching principle.The basic idea of the matched asymptotic expansion method [18] is that although the asymptotic solution cannot be given by the expansion of a single scale, it can be given by the expansion of different scales.These solutions are valid in their regions, and the whole region can be represented by the union of these regions.Then they are matched by the overlapping parts.Finally, the composite asymptotic solution is obtained uniformly and efficiently in the whole region.
The method of upper and lower solutions has been developed for more than a century with the aim of studying boundary value problems associated with ordinary and partial differential equations of various types.Since 1893, Picard [19] introduced lower and upper solutions to prove the existence of solutions for separated boundary value problems associated with scalar second-order ordinary differential equations.The theory was later developed by Nagumo [20] in the 1930s into differential inequality theory and then developed by Jackson [21].Howes [22] simplified their results.By using the comparison principle and the monotone iterative technique [23,24] combined with the method of upper and lower solutions, Cherpion et al. [25] obtained the existence and iterative methods of extremal solution to the system.The method of upper and lower solutions and monotone iteration technique are important and interesting tools to study the existence of the solution of nonlinear problems.
The monotone iterative method is a constructive approach that yields both the maximum and minimum solutions to a differential equation problem.The primary objective of the monotone iterative technique is to convert the nonlinear problem into a linear problem by the following specific steps: firstly, constructing an approximate solution sequence, secondly, demonstrating that this sequence is monotonically bounded and possesses a convergent limit, and, finally, proving that this limit corresponds to the solution of the original problem.
The utilization of differential inequality techniques in the investigation of singular perturbation problems for ordinary and partial differential equations has a concise yet captivating history.Later, researchers have extended this theory with Neumann and Robin boundary data and to second-order systems with fairly general types of boundary conditions.
The requirements for the smoothness of the upper and lower solutions are extended by allowing the upper and lower solutions to have a limited number of corner points.As long as the upper and lower solutions satisfy a certain size relation at the corner points, the theoretical system of second-order differential inequalities is basically formed.The differential inequality theory has become an important means to deal with the singular perturbation boundary value problem.Its characteristic is that it can not only prove the existence of the solution of the perturbation problem, but also obtain an accurate estimation of the perturbation solution by constructing appropriate inequalities.By using differential inequality theory the results proved by other methods can be obtained concisely and effectively.It can also deal with more complex problems and reveal the nature of their progressive processes, and this research method has been found to be useful in many different applications.
In the theory of differential inequalities in singular perturbation [26], the combustion problem is briefly introduced.Scholars [27,28] combine the matched asymptotic expansion method with the corner layer problem to study this type of problem and then use differential inequality theory to prove it.
In the Introduction to Perturbation Methods, Holmes used the matched asymptotic expansion method to study the famous corner layer problem, including the Prandtl matching principle and the intermediate variable matching method.The intermediate variable matching method is also mentioned by Holmes in studying Kummer functions.In general, this method involves segmenting for matching and introducing a new scale variable between the inner and outer layer variables before matching.Subsequently, Tingting et al. [29] considered a broader singularly perturbed boundary and corner layer problem with two turning points.By using the matched asymptotic expansion method scholars [30] discussed a class of nonlinear singularly perturbed problems with changed position of boundary layer and showed that the solution phenomena of singularly perturbed problems have full richness.When studying the oxygen diffusion problem, Boureghda [31] not only analyzed the numerical scheme theoretically, but also provided some numerical experiments with comparisons with analytical solution.
In this paper, we consider a second-order singularly perturbed boundary and interior layers problem with multiple turning points: where -1 ≤ x ≤ 1, 0 < a < 1, k ≤ -1, 0 < ε 1, and A, B = 0.This is a singularly perturbed turning point problem with three parameters.The properties and behaviors of the solutions will change according to the different values of the parameters, and the boundary and interior layer phenomena will appear from different positions.For singularly perturbed problems with turning points, few studies analyze the problem as a whole based on asymptotic solution and numerical fitting and prove the consistent validity of the solution through differential inequality theory.
In this paper, by analyzing the phenomena of boundary and inner layers, we use a variety of matching methods to solve them and then construct the upper and lower solutions to verify the existence of solutions.Finally, we obtain a composite expansion with good numerical solution fitting effect, which enabled us to have a deep understanding of this kind of singularly perturbed boundary and interior layer problems with multiple parameters.
The paper is organized as follows.In Sect.2, we divide the original problem into two main cases for matching and solving.We construct the asymptotic solution of the original problem and analyze the phenomena of the solution at the boundary layer and turning points.We use the relevant matching method to get the asymptotic solutions that are uniformly valid.In particular, the solution has a strong correlation with two kinds of special functions.In Sect.3, we provide an estimate of the remainder and demonstrate the existence of solutions with the help of the generalized Nagumo theorem proposed by Howes.We took advantage of the comparison principle and monotone iteration technique combined with the upper and lower solution method to obtain the existence of the extreme value solution.The requirements on the smooth degree of the upper and lower solutions are extended.By constructing the upper and lower solutions in the well-ordered case the theorem provides a way to prove the singular perturbation boundary and interior layers problem with multiple turning points.In Sect.4, we utilize the bvp4c solver in Matlab to investigate the numerical solutions for various specific instances.To further emphasize the significance of our research, we employ both asymptotic and numerical solutions.Finally, in Sect.5, we summarize some conclusions and possible areas of future research.

Main result
In this section, we discuss the singularly perturbed boundary and interior layers problems with multiple turning points, which can be categorized into two-turning-point and threeturning-point cases.During the solving process, we observe a strong correlation between singular perturbation problems with turning points and special functions.To address this issue, we employ the matched asymptotic expansion method to obtain a composite expansion within the corresponding interval.
To illustrate and further solve the problem, we will first scrutinize the potential phenomena that may arise.This problem is characterized by three parameters n, a, and k.
It is obvious that any alterations to these values will correspondingly impact the original problem.
Firstly, since 0 < a < 1, we have -a < 0, 1a > 0, and -1 < -a < 0 < 1a < 1, and Eq.(1) has two turning points x = -a and x = 1-a.The original interval is divided into three parts by these two points.On these three intervals, the sign of the coefficient at y is > 0, < 0, and > 0 respectively, that is, the algebraic sign of the coefficient of y changes at these two points.So we guess that the corresponding interior layer phenomena will occur at these points.
Secondly, since x ∈ [-1, 1] and k ≤ -1, when k < -1, (xk) 2n+1 > 0 is always true, and when k = -1, x = -1 becomes a higher-order turning point of this problem.According to the value of k, the original problem is divided into two-and three-turning-point problems.The parameter k makes the discussion of the solution of the original problem different.
Finally, the parameter n determines the power of the small parameter ε and (xk), as well as the thickness of the boundary layer and interior layer.When n = 0 and k = -1, x = -1 is a first-order turning point of the original problem, whereas when n ≥ 1 and k = -1, x = -1 is a higher-order turning point of the original problem.
We seek the outer expansion in the form where n = 0, 1, 2, . . ., and then substituting Eq. ( 4) into Eq.( 1), we obtain By comparing the coefficients at of ε on the left and right sides, we obtain the reduced problem which is a first-order differential equation that can be determined with only one boundary condition.
We commence by ascertaining the presence of boundary layers and interior layers.According to the previous analysis, the coefficient of y on (-1, -a) is > 0, and the layer phenomenon appears near the left endpoint x = -1.The coefficient at y on (1a, 1) is also > 0, so there is no layer phenomenon near the right endpoint x = 1.Since the coefficient at y on (-a, 1a) is < 0, there is a layer phenomenon near the turning point x = 1a.Therefore we speculate that the original problem exhibits boundary layer phenomenon near x = -1 and interior layer phenomenon near x = 1a.The specific type of interior layer phenomenon near x = 1a needs further analysis.
We hypothesize the presence of an interior layer phenomenon at the turning point x = -a.However, our analysis reveals no corresponding layer phenomenon near x = -a.We can confirm that the problem will not exhibit boundary layer phenomena near x = 1.Therefore the reduced problem (5) satisfies condition (3).In addition, we can note that y 0 = 0 is also a solution of Eq. ( 5).Thus we obtain the equation The presence of a singularity at the point x = -a is evident in this solution.To eliminate this singularity and avoid subsequent complications, it is necessary to consider the previous analysis on layer phenomena classification.Consequently, we can derive the zerodegree outer solution The outer solution y 0 is continuous, but it lacks differentiability at x = 1a.For this reason, there is a corner layer at x = 1a.In general, the exact solution is not zero at the turning point x = -a, so we guess that the spike layer appears at x = -a, that is, the pulse state.At this time, we can find that y 0 (-1) = 0 = A, which also confirms that the boundary layer does appear at x = -1.In the following, we construct boundary, spike layer, and corner layer expansions at x = -1, x = -a, and x = 1a, respectively.The matched asymptotic expansion method is used to match them with the outer expansions so as to obtain the uniformly valid composite expansions in the interval [-1, 1].
The solution is first analyzed and discussed at the point x = -a.Given the suspicion of a spike layer phenomenon [32,33] at this particular point, the behavior of the solution within the interval (-1, 1a) can be described as follows: where s = 0.For a semilinear boundary value problem [32], where x = 0 is a turning point, that is, f (0) = 0. We assume that (a) f (x) and g(x, y) are continuously differentiable on [-1, 1] and [-1, 1] × R, respectively, where x = 0 is the turning point, f (0) = 0, and there exists d > 0 such that 0 < -f (0) ≤ d.
For Eq. ( 1) with f (x) = -(xk) 2n+1 (x + a)(x -(1a) and g(x, y) = (xk) 2n+1 y, we will verify whether the pulse state occurs at x = -a; we can prove that f (-a) = -ak > 0, where k ≤ -1, so that this problem does not satisfy the spike layer existence theorem at x = -a.Hence, when s = 0, indicating that f (-a) = 0, our focus lies in constructing boundary layer and corner layer expansions at x = -1 and x = 1a, respectively, with the intention of subsequently matching them.
Remark 1 If problem (1)-( 3) is divided into left and right problems according to x = -a for analysis, then we can also obtain that there is no any layer phenomenon at x = -a and f (-a) = 0.
Remark 2 It is worth noting that no matter k < -1 or k = -1, the boundary layer appears near the left endpoint, but the specific equations of the boundary layer are completely different in these two cases, and the thickness of the boundary layer is also different.Therefore it is necessary to discuss these cases separately.

Case 1. k < -1
We obtained the outer solution (7) for the original problem ( 1)-( 3), but it could not satisfy boundary condition ( 2)), and thus we need to construct the boundary layer solution to correct it near x = -1.When k < -1, the coefficient at y is constant > 0, and we introduce the stretched variable ξ = x+1 ε v > 0. Denoting the inner solution by Y and making the corresponding change of variables in Eq. ( 1), we get and the main part is From the principle of minimum degradation we obtain v = 2n + 1.We look for the asymptotic of the solution in the form of the inner expansion, and thus near the origin, the equation has the form Then substituting Eq. ( 10) into Eq.( 9) and comparing the coefficients of O( 1 ε 2n+1 ) on the left and right sides, we get Then Eq. ( 11) becomes Y 0 + mY 0 = 0, so thatwe obtain where c 1 , c 2 are arbitrary constants, and Y (ξ ) satisfies condition (2).When x = -1 and ξ = 0, we obtain Y 0 (0) = A and c 1 + c 2 = A.This must match with the outer solution.The matching condition by the matching asymptotic expansion method is Y 0 (+∞) = y 0 (-1 + ), and through matching we have that c 1 = 0. Therefore we determine the boundary layer solution where ξ = x+1 ε 2n+1 > 0. To determine the solution in the corner layer at x = 1a, we introduce the corner-layer stretched variable η = x- (1-a)  ε v .We denote the solution in this region by Ỹ .Substituting it into Eq.( 1), we obtain and the main part is The distinguished limit in this case occurs when v = 2n+1 2 .Unlike our previous assumption, we now consider the asymptotic expansion of corner solution For this problem, y 0 (1a) = 0, and the constant m will be determined from the matching.
Remark 3 Assuming the previous asymptotic expansion of corner solution, that is, we can also get the first term Ỹ00 = 0. We will verify that the unknown constant m = 2n+1 2 .Expansion (15) is the unique form of the asymptotic expansion of the corner layer solution, because the outer solution is continuous but nondifferentiable at the corner point.The solution from the corner-layer region should provide a smooth transition between these linear functions.Thus, substituting Eq. ( 15) into Eq.( 14) and comparing the coefficients of O(ε m ), we have It is possible to solve this equation using power series methods.However, a simpler way is to notice that Ỹ0 = η is a particular solution.Then using the method of order reduction, we find that the general solution is where C and D are arbitrary constants.
Remark 4 Since the error function [34] is defined as erf(x) = 2 √ π x 0 e -y 2 dy, Eq. ( 17) is equivalent to the expression without integration.Then the corner layer solution is Now we will use the matching asymptotic expansion method to match the corner solution with the outer solution, whereas in the previous part, the boundary layer analysis and matching were performed simultaneously.However, this approach is not applicable here as the matching process in this layer exhibits slight variations that warrant further consideration.To facilitate the matching procedure, we introduce the intermediate variable ε p , where 0 < p < 2n+1 2 .The outer solution ( 7) can be reformulated in this variable as Since x η ε p ∈ (0, a) ⊂ (0, 1), 1 + x η ε p ∈ (1, a + 1), and thus the variable a ∈ (0, 1) suggests that as a approaches 0, the approximate effect is enhanced.Subsequently, we will conduct matching verification by comparing the fitting effect of both asymptotic and numerical solutions.Since +∞ 0 exp(-t 2  2 ) dt = π 2 , we rewrite Eq. ( 18) as follows: 2 ), x η < 0, To match Eqs. ( 19) and ( 20), we must have m = 2n+1 2 .In this case, we obtain C = B(1+a) 2a and D = - by the intermediate variable matching principle.The final step involves consolidating these expansions into a composite expression.This is done in the usual way of adding the expansions together and then subtracting the common parts.The construction of a composite expansion for and for x ∈ (1a, 1], Remark 5 The error function was ultimately employed to represent the composite solution, not only for the convenience of numerical fitting, but also to enhance its clarity and intuitive presentation.Now analyzing the properties of the corner point x = 1a, we find that and as ε approaches 0, By the definition of the derivative we obtain Since the left and right derivatives exist and are equal, it is known that y(x, ε) is continuous and differentiable at x = 1a.Thus Eqs. ( 21) and ( 22) are the composite expansions over the entire interval [-1, 1].

Case 2. k = -1
In this case, the original problem becomes and (x + 1) 2n+1 = 0 when x = -1.Then x = -1 is not only a boundary layer point, but also a boundary turning point.Therefore there may be a multilayer phenomenon near this point; at least, the boundary layer solution will be different from Case 1.
In this case the outer expansion is still (7).Based on the preceding analysis, we investigate the possibility of the boundary layer.This is done by introducing the boundary-layer coordinate ξ = x+1 ε v , that is, x = ε v ξ -1.Also, as in the previous case, we denote by Y the solution in the layer.Now substitution of it into Eq.( 23) yields By sorting it out we get and the balance in this layer is between these three terms.We first match 2n + 1 -2v = 2nv and get v = 2n+1 2n+2 .Then Substituting Eq. ( 26) into Eq.( 25) and comparing the same power coefficients of ε, we get and then we obtain Y 0 + qξ 1 2n+1 Y 0 = 0 by noting q = (-1 + a)(-2 + a) > 0. The solution is Now we match 2n +1-2v = (2n +1)v and get v = 2n+1 2n+3 .Then and then we get The general solution of the latter is Ȳ0 = C, where C is an arbitrary constant.Because of 2n+1 2n+3 < 2n+1 2n+2 , ξ 1 is the boundary layer extension variable, and ξ 2 is the intermediate layer extension variable.We find that the boundary layer solution and the outer solution can be matched directly, so we guess that the intermediate layer solution C = 0.After the matching verification of the intermediate layer solution and the outer solution, we get C = 0, that is, there is no intermediate layer in this situation.So Eq. ( 28) must match with the outer solution.The matching condition is Y 0 (+∞) = y 0 (-1 + ), and Eq. ( 28) also satisfies the boundary condition (2).Then we have By the definition of gamma function, (x) = +∞ 0 t x-1 exp(-t) dt, we obtain a 0 = -q . Then we determine the following boundary layer solution: In this case, it is easy to obtain the corner layer solution as Similarly, the determination of the first-term composite expansion of the solution over the entire interval can be inferred from the aforementioned equations.
The inclusion of integrals in the form of solutions is not recommended, as mentioned in Remark 5.By utilizing incomplete gamma functions and their related properties [35] we can transform the solution from an integral form into a nonintegral form.We get ) when n = 0 and , the error function is one of the incomplete gamma functions, that is, the case n ≥ 1 includes the case n = 0. Therefore, for x ∈ [-1, 1a], the first-term composite expansion without integral form of the solution over the entire interval is and when x ∈ (1a, 1], Remark 7 The incomplete gamma function studied in this section is related to the confluent hypergeometric function [36].It has received a great deal of attention because the function has been found to have some rather interesting properties.Since γ (s, z) = z s s M(s, s + 1, -z), Eq. ( 34) can also be written as

The existence of solutions and estimation of the remainder
In this section, by using the generalized Nagumo theorem proposed by Howes and the monotone iteration technique combined with the upper and lower solution we obtain the existence of the solutions and estimate the remainder.In recent years, the theory of differential inequalities has become an important means to deal with (non)linear singular perturbation problems.It can not only prove the existence of perturbation problems, but also obtain accurate estimates of perturbation solutions by constructing appropriate upper and lower solutions.
Considering first the general Dirichlet problem (37) where f is a continuous function on [a, b] × R 2 , the method is based on the observation that there exist smooth functions α(t) and β(t) such that Theorem 1 [15] Assume that there exist upper and lower solution functions α(t) and β(t) belonging to C (2)   (2) on [a, b], then as long as α ≥ f (t, α, α ) and β ≤ f (t, β, β ) on every subinterval (t i-1 , t i ) and α (t -) ≤ α (t + ) and β (t -) ≥ β (t + ) for all t ∈ (a, b), then the corresponding result is valid.
The generalized Nagumo condition is that there exists a continuous function ϕ(s) > 0 on [0, +∞) such that for a ≤ t ≤ b, α(t) ≤ y ≤ β(t), and |y | < +∞, we have |f (t, y, y The above theorem has been generalized in several directions.Cherpion et al. [37] proved an existence result in the well-ordered case and nonwell-ordered case.For wellordered cases, the well-ordered lower and upper solutions are given by the maximum principle, whereas for nonwell-ordered cases, the lower and upper solutions in the reversed order are given by the antimaximum principle.They approximate the sequence by replacing global Lipschitz conditions with local Lipschitz conditions and adding Nagumo conditions.Then the existence of convergent sequences is obtained again by replacing two-sided Nagumo conditions with one-sided conditions.The monotone iteration technology has been paid attention to and developed continuously.De Coster et al. studied the monotone iteration technology in boundary value problems.
In addition, the method of proof requires a growth restriction.Howes employed the Nagumo-Jackson theorem as well as an extended version of the lower and upper solution.Two sets of monotone sequences are defined in the proof.According to the continuity of the original problem and the monotonic boundedness of the sequences, it is easy to show that they are uniformly convergent.The minimum and maximum solutions are found, and the contraction mapping principle shows that the solutions are unique.
Theorem 1 extends the Nagumo-Jackson theorem described above by allowing the bounding functions α and β to have finitely many "corners", provided that the correct inequalities are satisfied at the corner points t i .Its proof is an easy modification of the given one if we note that the maximum of a finite number of lower solutions (i.e., functions satisfying the α-differential inequalities) is also a lower solution, whereas the minimum of a finite number of upper solutions (i.e., functions satisfying the β-differential inequalities) is an upper solution.
We use the ideas introduced and developed in the above paper to prove other problems.Consider the corner layer behavior of a general second-order nonlinear boundary value problem Definition 1 Let constants δ 2 > 0 and δ 3 > 0 be such that f y (t, y, y , where Then the solution u = u(t) of the degenerate equation f (t, u, u ) = 0 of Eq. ( 39) is said to be locally weakly stable around t = t 0 ∈ (a, b).
Definition 2 Suppose there exist constants m > 0 and δ 1 > 0 such that the function h(t, y) = f (t, y, u ) has a partial derivative of order 2q + 1 in and Then the solution u = u(t) ∈ C 2 [a, b] of the degenerate equation f (t, u, u ) = 0 of Eq. ( 39) is said to be (I q ) stable on [a, b].
If [H1] and [H2] are satisfied, and the reduced solution u = u(t) is locally weakly stable near t = t 0 and (I q ) stable on [a, b], then we assume that u L (t 0 ) < u R (t 0 ).To estimate the solution using the theory of differential inequalities, it is necessary to construct the lower and upper solution functions and where the undetermined positive function V (t, ε) is the inner function operating near t = t 0 , and C is the undetermined constant.(The u L (t 0 ) > u R (t 0 ) also has a corresponding result.)Now we prove the problem studied in this paper.We only prove the case of n = 0, that is, For n = 1, 2, . . ., the same result can be proved.
Theorem 2 With Nagumo condition satisfied, there exists a sufficiently small positive number ε 0 such that for every 0 < ε ≤ ε 0 , problem (40)-( 41) has a solution y(x, ε) with the boundary layer property at x = -1 and the corner layer property at x = 1a, and as ε approaches 0, we have where × γ ( .
Proof A differentiable function is continuous and bounded on any closed interval.From There exists M > 0 such that |F(x, y, y )| ≤ M(1 + |y |), and under the generalized Nagumo condition, Since the coefficient at y in Eq. ( 40) changes sign at the points x = -a, 1a, the left boundary layer has little influence on the right solution by the boundary layer function method.Moreover, from the previous analysis we know that the point x = -a is a turning point that has little influence on the original problem, and the value at this point can only be 0, that is, y(-a) = 0, and a corner layer is generated near x = 1a.Then we divide the problem into left and right problems, as well as the interval into left and right intervals, that is, to [-1, -a] and [-a, 1].
The left problem is and the right problem is Then the left problem is a general singular perturbation problem with a left boundary layer, and the right problem is a singular perturbation problem with a corner layer.
For the left problem, we can construct lower and upper solutions as follows: so that, when x ∈ [-1, -a], as long as γ > 0, all verification conditions can be met.For ε small enough, the left problem has a solution y L (x, ε) on [-1, -a] that satisfies α L (x, ε) ≤ y L (x, ε) ≤ β L (x, ε), and |y L (x, ε) -A exp(- since the solution of the reduced problem at this time is locally weakly stable near x = 1-a, that is, there is an arbitrarily small constant δ 1 > 0 such that F y (x, y, y Then the reduced solution is (I 0 ) stable on the interval [-a, 1], that is, there exists constants m > 0 and δ 2 > 0 such that F(x, u, u ) = 0 and so that the lower and upper solutions to the right problem can be constructed based on these properties. When As long as γ ≥ 2B(1+a) (1-a-k)a , all conditions are established.Similarly, when u L (1-a) > u R (1-a), that is, when B < 0, we have and γ needs to satisfy the condition γ ≥ -2B(1+a) (1-a-k)a .Then when x ∈ [-a, 1], as long as γ ≥ 2|B|(1+a) (1-a-k)a , all the test conditions are satisfied.Thus as long as ε is small enough, there exists a solution y R (x, ε) on [-a, 1] satisfying α R (x, ε) ≤ y R (x, ε) ≤ β R (x, ε), and as ε approaches 0, ), x > 1a.
Remark 8 When employing the theory of differential inequalities to establish the existence of a solution, it is permissible for all functions to incorporate a boundary layer at the endpoint.However, their stability is compromised upon traversing the turning point.
In this paper, the presence of a boundary layer in the left problem does not influence the occurrence of the corner layer phenomenon in the right problem.Consequently, by combining these two problems we establish the existence of solutions and provide estimations for the remaining factors.As a result, Theorem 2 has been proven.

Numerical examples
In this section, to further investigate the approximate effect of the asymptotic solution obtained, we explore examples featuring a boundary layer at the left end of the interval and a corner layer inside it.We consider four specific cases, including the two cases studied in Sect.2, where multiple possibilities are discussed for each case.

Example 1
When k = -1 and n = 0, we provide an example with a corner layer at x = 1 2 .Consider the problem As shown in Fig. 1, we can roughly find the demarcation point p * between the boundary layer and the left outer layer and then represent the region of the boundary layer from Fig. 2(a).Similarly, the two demarcation points q l * and q r * between the outer and corner layers can also be found, and the general area of the corner layer can be marked from Fig. 1.
The fitting figures of the numerical and asymptotic solutions are shown in Figs.2(a) and 2(b) with ε = 10 -2 and ε = 10 -3 , respectively.From Fig. 2(a) we can see that x = 1 2 is the corner point, and the thickness of boundary and corner layers is O(ε 1 2 ).The smaller the ε, the better the fitting effect, and it is obvious that the thickness of boundary and corner layers also becomes thinner with the decrease of ε.
On the basis of Example 1, we choose a = 1 4 , A = -3, B = 2, and ε = 10 -3 to make Fig. 3, which further verifies that the asymptotic solution has a good approximation.

Example 2
Unlike Example 1, we take n = 1 to study the following problem with a = 1 2 : As shown in Fig. 4(a), the thickness of boundary and corner layers are changed in the example, which now are O(ε 3 4 ) and O(ε 3 2 ), respectively.We take ε = 10 -1 and 10 -2 , respectively, and observe the changes in Fig. 4. We find that changing the value of n does not affect an excellent approximation, but it does influence the thickness of each layer.

Example 3
Now we take k = -1.5 and n = 0 to investigate the problem We choose a = 1-a = 1 2 and ε = 10 -2 , 10 -3 , respectively.As shown in Fig. 6, the thickness of the boundary layer is O(ε), and the thickness of the corner layer is O(ε 1 2 ).We mark the layer phenomena in Fig. 5, and the numerical and asymptotic solutions fit well.As we calculated previously, the difference between k < -1 and k = -1 is in the thickness of the left boundary layer.
If we take a = 1 8 and a = 7 8 , then the corner points are x = 7 8 and 1 8 , respectively.As shown in Figs.7 and 8, when ε = 10 -2 , it is obvious that the fitting effect of a = 1  8 is better than that of a = 7  8 , but by adjusting ε = 10 -3 the fitting effect of a = 7 8 will also be better.This also confirms our previous conjecture that the closer the a to 0, the better its asymptotic effect.To verify the reliability of the conjecture, if a approaches 0, then the composite solution can be approximated as a solution without singularity, and the conjecture can also be confirmed by comparison of Figs. 7 and 8.We now consider two critical cases where a = 0 and a = 1 in problem (1): The solutions are shown in Fig. 9. Comparing Fig. 9 with Figs. 7 and 8, we clearly see that the trend of the figures is consistent, which further confirms the correctness of our previous conjecture.

Example 4
We take n = 1 to study the problem (45) As we discussed in Example 2, the thickness of boundary and corner layers changes; the thickness of the boundary layer is O(ε 3 ), and the thickness of the corner layer is O(ε 3 2 ).It is worth mentioning that since the value of n determines the power of ε and (xk) in problem (45), when n = 1, the order of ε is 3, as shown in Fig. 10, and then the fitting effect of ε = 10 -1 is already good.

Conclusions and future research plan
In this paper, we discuss the singularly perturbed boundary and interior layers problems with multiple turning points, specifically focusing on the cases involving two or three turning points.Subsequently, we derive composite expansions that effectively capture the numerical solutions within their respective intervals.Furthermore, based on our analysis process and solution method, we can extend our discussion to scenarios featuring a greater number of turning points.
In the process of solving, we observe a strong correlation between singular perturbation problems with turning points and special functions.This paper explores various special functions such as the error function, the (incomplete) gamma functions, and the confluent hypergeometric function.These special functions possess specific definitions and properties, which are often relevant to the equations we need to solve.Moreover, the inner and outer solutions can also be matched with their related properties.Given the diverse range of special functions available, it is crucial to consider any potential restrictive conditions when utilizing them.
Compared to the previous studies on this type of problem, in this paper, we present a relatively comprehensive theoretical framework for singular perturbation analysis.It is widely acknowledged that Eq. ( 1) exhibits a distinct spike phenomenon at x = -a.However, through rigorous verification, this paper demonstrates that although x = -a acts as a turning point altering the sign of the original problem, it does not give rise to an inner layer phenomenon.In fact, such turning point problems often exhibit such characteristics.Throughout the process of solving corner layer issues, this paper effectively showcases the adaptability and versatility of matching methods.
There are many matching methods in singular perturbation problems.However, when the general matching method proves inadequate, alternative approaches become necessary.The intermediate variable matching method serves as an effective tool for addressing complex singular perturbation problems by establishing a connection between the inner and outer expansions through careful selection of appropriate scale variables.So it can deal with some complicated singular perturbation problems more effectively.The combination of intermediate variable matching method and special functions can well work out this kind of interior problem.
After proving the existence of a solution, estimating the remainder and verifying the asymptotic effect by numerical examples, we find that the zero-order approximation obtained in this paper is good, and we try to continue to obtain higher-order approximations.When k < -1, there is no special function in the boundary layer solution, and we can solve the boundary layer first-order solutions at this time.For example, when n = 0, we take a = 1  2 , and as shown in Fig. 11, the asymptotic effect of the first asymptotic solution (the solid red line) on the interval [-1, -a] is slightly better than that of the zeroth asymptotic solution(the dotted green line).Therefore we guess that the higher-order asymptotic solution will improve the approximation of this problem.However, in the process of solving, the singularity will be stronger, and the appearance of various special functions will increase the difficulty of solving the solutions.Solving higher-order asymptotic solutions and eliminating singularity is our direction in the future.
Finally, if we can generalize the problem to a singularly perturbed problem with an increased number of turning points or higher nonlinearity, it is anticipated that the complexity of the problem will be significantly amplified.However, this plan remains a potential avenue for future exploration.
We can also consider problem (46)-(47): where k ≥ 1.Through analysis we know that this problem will produce a boundary layer phenomenon at the right endpoint and a corner layer phenomenon at x = -a, and the asymptotic effect is better as a approaches 1.As shown in Fig. 12, we let a = 0.25, 0.5, 0.75, respectively, and observe their properties; just as we have analyzed, this problem has a right boundary layer and an inner corner layer at x = -a.
) assuming the following conditions: (H1) f , f y , f y ∈ C(D(u)), where D(u) is the same as in the following definition; f y (t, y, y ) is bounded in a compact subset of [a, b] × R, and for (t, y) ∈ [a, b] × R, we have f (t, y, y ) = O(y 2 ), |y | → +∞.(H2)The degenerate equation f (t, u, u ) = 0 has a pair of solutions u and u R (b) = B, where t 0 ∈ (a, b), and we let

Figure 8 8 Figure 9
Figure 8 Comparison of the asymptotic and numerical solutions when k = -1.5, n = 0, and a = 7 8

Remark 9
When n = 0 and k < -1, we can try to use the left solution to figure out where the boundary layer is at the point.Since the boundary layer extension variable is ξ = x+1 ε and the position of the boundary layer is x = -1+εl, we obtain l = --1-k mh by solving the algebraic equation when substituting it into the left solution, where h = (-1k)(-1 + a)(-2 + a), and m = (-1k)(-1 + a) + (-1k)(-2 + a) + (-1 + a)(-2 + a).Then for a = 1 2 and k = -1.5, we obtain l = 5.3, which can also further determine the location of the boundary layer.

Figure 11
Figure 11First-order approximation of boundary layer