Solving Third-Order Singularly Perturbed Problems : Exponentially and Polynomially Fitted Trial Functions

For the third-order linearly singularly perturbed problems under four different types boundary conditions, we develop a weak-form integral equation method (WFIEM) to find the singular solution. In the WFIEM the exponentially and polynomially fitted trial functions are designed to satisfy the boundary conditions automatically, while the test functions satisfy the adjoint boundary conditions exactly. The WFIEM provides accurate and stable solutions to the highly singular third-order problems.


Introduction
In this paper we propose a new numerical method for solving the boundary value problem of third-order ordinary differential equation whose third-order derivative term is multiplied by a small perturbing parameter, which is one of the so-called singularly perturbed boundary value problems (SPBVPs).The SPBVPs exhibit boundary layer behavior, which is a narrow region with the solution varying rapidly.
The arrangement of this paper is given as follows.In Section 2 we introduce a third-order linearly singularly perturbed problem.In Section 3 we develop a weak-form method for both the constant and varying coefficients third-order linearly singularly perturbed problems.In Section 4 we expand the trial solution in terms of exponentially and polynomially fitted trial functions as bases, which are designed to satisfy the boundary conditions automatically, and a numerical algorithm based on the resultant weak-form integral equation method (WFIEM) is developed in Section 5.The numerical experiments are given in Section 6 to validate the performance and effectiveness of the present method.Some conclusions are drawn in Section 7.

Third-order Linearly Singularly Perturbed Problems
Let us consider a third-order linearly singularly perturbed ordinary differential equation, describing a boundary layer problem under a varying source (Valarmathi & Ramanujam, 2002;Valanarasu & Ramanujam, 2007): where ε > 0 is a small parameter.There are four types boundary conditions: (5) For cases (a) and (b) we suppose that the boundary values α and β satisfy α 2 +β 2 > 0. For a given problem if α = β = 0 we can transform it by u(x) = y(x) + x to a new problem for u(x) with the new boundary conditions u(0) = 0 and u(ℓ) = ℓ > 0.
For cases (c) and (d) we suppose that γ 0.
The classification of high-order singularly perturbed problems to the type of reaction-diffusion and convection-diffusion depends on how the order is affected if ε = 0.If the order is reduced by one, i.e., a 0, the problem is of convectiondiffusion type, and of reaction-diffusion type if the order is reduced by two, i.e., a = 0 and b 0.

A Weak Form Integral Equation Method
Applying a test function v(x) to Equation (1) and integrating it from x = 0 to x = ℓ, we can obtain a weak-form integral equation: Integration by parts we can see that the first term gives Since in the boundary conditions both y ′′ (ℓ) and y ′′ (0) are not given, we must take the adjoint boundary conditions: v(0) = v(ℓ) = 0 for the test function in order to eliminate y ′′ (ℓ)v(ℓ) and y ′′ (0)v(0) on the right-hand side of Equation ( 7), of which the simplest ones are the sine functions: Through a similar procedure for other terms in Equation ( 6) and using v(0) = v(ℓ) = 0 we can derive a weak-form integral equation for y(x): The adjoint boundary conditions for v are derived by eliminating the boundary terms of y which are not given in the boundary conditions.According to this principle, by matching the adjoint boundary conditions we can derive the following test functions for the four types boundary conditions in Equations ( 2)-( 5): Equation ( 8) is easily extended to the non-autonomous ODE of Equation ( 1) with a, b and c being functions of x.For example if b is a function of x we can write where the conditions v(0) = v(ℓ) = 0 were used.Hence, the term − ∫ ℓ 0 bv ′ (x)y(x)dx on the left-hand side of Equation ( 8) Therefore, through some manipulations we can obtain the following weak-form integral equation for the non-autonomous third-order ODE: where we have imposed the conditions v(0) = v(ℓ) = 0, and for the above four types boundary conditions, other terms can be eliminated further by using Equations ( 9)-( 12).However, we do not write them in detail.
Some authors (Valanarasu & Ramanujam, 2007;Miller, O'Riordan & Wang, 2000;Roos & Zarin, 2002;Farrel, Hegarty, Miller, O'Riordan & Shishkin, 2004) have considered the singularly perturbed problems with discontinuous source terms, which lead to a weak interior layer in the singular solution.It is interesting that the present weak-form integral equation method can handle this type problem without inducing any difficulty because the influence of the source term f (x) on the solution is given through the integral term

The Trial Functions
In this section we describe a simple method to solve y(x) in the integral equations ( 8) and ( 13).We suppose that the solution y(x) can be expanded by a set of trial functions ϕ j (x): where a j are unknown coefficients to be determined, and the scales s j are determined by the equilibrated method (Liu, 2012b).
In the weak-form formulation of differential equations it is known that the selection of trial functions is very important, of which we suppose that the set of trial functions is complete, the trial functions must be linearly independent, and the trial functions should satisfy the boundary conditions exactly (Fletcher, 1984).The most often used trial functions are trigonometric polynomials, Chebyshev polynomials and Legendre polynomials.We will give a different set of trial functions which are not used in the open literature to treat the third-order singular problems.
To prompt the introduction of the exponential type trial function, let us consider a simple second-order singular problem: which has a closed-form solution: It can be seen that the singular solution is of the exponential type function.
There are trial functions which automatically satisfy ϕ j (0) = 0 and ϕ j (ℓ) = 1: Similarly, we have They are two basic elements to generate the trial functions which automatically satisfy the given boundary conditions in Equations ( 2)-( 5): In above, besides the exponential terms the low order polynomials of x were used to adjust the trial functions, such that all the required boundary conditions are satisfied exactly.Due to the above properties ϕ j (x) are called exponentially and polynomially fitted trial functions.For the special case j = 0 we can derive ϕ 0 (x) by using the L'Hosptial rule for each function ϕ j (x) in the above.

Numerical Algorithm of WFIEM
In this section we describe the resultant algorithm of the weak-form integral equation method (WFIEM) to solve y(x).
Inserting the expansion of Equation ( 14) into Equations ( 8) and ( 13) with v replaced by v k and letting k = 1, . . ., n q we can derive a linear system: to determine the expansion coefficients a := {a j } whose number is denoted by n = m 1 + m 2 + 1.On the other hand we impose an extra moment equation: to guarantee that the trial solution satisfies the given boundary conditions exactly.The dimension of A is (n q + 1) × n, such that Equation ( 23) is an over-determined system with n q + 1 > n.
Instead of Equation ( 23), we can solve a normal linear system: where The algorithm of conjugate gradient method (CGM) for solving Equation ( 25) is summarized as follows.
(i) Give an initial c 0 and then compute r 0 = Dc 0 − b 1 and set p 0 = r 0 .

Numerical Algorithm of WFIEM
Example 1.In this example we test the performance of the weak-form integral equation method (WFIEM) by solving which is a reaction-diffusion problem with the boundary conditions of type (a), and of which a uniform asymptotic approximate solution is given by For this problem a = c = 0, b = −1 and f (x) = 1.The numerical procedure for this problem to generate A in Equation ( 23) is inserting Equation ( 19) for ϕ j (x), j = −m 1 , . . ., m 2 and v k (x) in Equation ( 9) into the left-hand side of Equation ( 8) and integrations for different k = 1, . . ., n q , of which we can obtain the all components A k j of A. Similarly, inserting v k (x) into the right-hand side of Equation ( 8) and integrations for different k = 1, . . ., n q we can obtain the components e k , k = 1, . . ., n q of the vector e in Equation ( 23).Then we apply the linear equations solver of CGM as specified in Section 5 to solve Equation ( 23), from which we can find all the expansion coefficients a j , j = −m 1 , . . ., m 2 .Then the numerical solution in Equation ( 14) is obtained.
For this problem a = 2, b = −4, c = 2 and f (x) is given by Equation (31).The numerical procedure for this problem to generate A in Equation ( 23) is inserting Equation ( 21) for ϕ j (x), j = −m 1 , . . ., m 2 and v k (x) in Equation ( 11) into the left-hand side of Equation ( 8) and integrations for different k = 1, . . ., n q , from which we can obtain the all components A k j of A. Similarly, inserting Equation ( 11) for v k (x) into the right-hand side of Equation ( 8) and integrations for different k = 1, . . ., n q we can obtain the components e k , k = 1, . . ., n q of the vector e in Equation ( 23).Then we solve Equation ( 23) to determine all the expansion coefficients a j , j = −m 1 , . . ., m 2 .Then the numerical solution in Equation ( 14) is obtained.
For the purpose of comparison we also apply the fourth-order Runge-Kutta (RK4) method to integrate Equation ( 30) by starting from the initial conditions y(0) = 1, y ′ (0) = 0 and y ′′ (0) = 115.5543325with a time stepsize ∆x = 0.001, which can match the final time value y ′ (1) = 0 with an error 8.6 × 10 −10 .Under the following parameters ε = 0.01, m 1 = 60, m 2 = 0, n q = 61, we solve this problem.From Fig. 2(a) we can see that the numerical solution y(x) is close to that obtained by the RK4 with the maximum difference being 4.47 × 10 −4 .This problem shows that the present method is highly accurate.As shown in Fig. 2(b) for the comparison of the derivatives they are also close very well.Example 3. In this example we consider (Valanarasu & Ramanujam, 2007) where This problem is a convection-diffusion problem with a weak interior layer, while the boundary conditions are of type (d).
For this problem a = 2e x , b(x) = − cos(4πx), c = 1 + x and f (x) is given by Equation (33).The numerical procedure for this problem to generate A in Equation ( 23) is inserting Equation ( 22) for ϕ j (x), j = −m 1 , . . ., m 2 and v k (x) in Equation ( 12) into the left-hand side of Equation ( 13) and integrations for different k = 1, . . ., n q , of which we can obtain the all components A k j of A. Similarly, inserting Equation ( 12) for v k (x) into the right-hand side of Equation ( 13) and integrations for different k = 1, . . ., n q we can obtain the components e k , k = 1, . . ., n q of the vector e in Equation ( 23).
Then we apply the CGM to solve Equation ( 23) to find all the expansion coefficients a j , j = −m 1 , . . ., m 2 .Then the numerical solution in Equation ( 14) is obtained.
Under the following parameters ε = 0.01, m 1 = 60, m 2 = 0, n q = 30, we solve this problem.From Fig. 3(a) we can find that the numerical solution y(x) is close to that obtained by the RK4 with the maximum difference being 1.49 × 10 −3 .This problem shows that the present method is highly accurate.To compare the derivatives in Fig. 3(b) they match very well, no matter in the boundary layer and in the interior layer near to x = 0.5.

Conclusions
In this paper we have developed a weak-form integral equation method (WFIEM) for the non-autonomous third-order linearly singularly perturbed boundary value problems under varying source.The numerical algorithm based on the WFIEM and the exponentially and polynomially fitted trial solutions was developed.When the trial functions were

Figure 1 .
Figure 1.For a singular problem of example 1, (a) showing the convergence iterations and (b) Figure 1.For a singular problem of example 1, (a) showing the convergence iterations and (b) comparing numerical solution and uniform approximation and showing the difference.

Figure 2 .
Figure 2.For a singular problem of example 2, (a) comparing numerical solution and the RK4 Figure 2.For a singular problem of example 2, (a) comparing numerical solution and the RK4 approximation and showing the difference, and (b) comparing the derivatives.

Figure 3 .
Figure 3.For a non-autonomous singular problem of example 3, (a) comparing numerical solutionFigure 3.For a non-autonomous singular problem of example 3, (a) comparing numerical solution and the RK4 approximation and showing the difference, and (b) comparing the derivatives.