Cubic Hermite Collocation Method for Solving Boundary Value Problems with Dirichlet, Neumann, and Robin Conditions

Abstract

Cubic Hermite collocation method is proposed to solve two point linear and nonlinear boundary value problems subject to Dirichlet, Neumann, and Robin conditions. Using several examples, it is shown that the scheme achieves the order of convergence as four, which is superior to various well known methods like finite difference method, finite volume method, orthogonal collocation method, and polynomial and nonpolynomial splines and B-spline method. Numerical results for both linear and nonlinear cases are presented to demonstrate the effectiveness of the scheme.

1. Introduction

Since time immemorial mathematicians are ardently pursuing the solution of linear or nonlinear two point boundary value problems (BVPs) of the type: 𝑑²𝑦𝑑𝑥²+𝛼₁(𝑥)𝑑𝑦𝑑𝑥+𝛼₂(𝑥)𝑦=𝑓(𝑥),𝑥∈[𝑎,𝑏] (1) subject to Dirichlet, Neumann, and Robin’s boundary conditions.

Such BVPs have wide application in astronomy, biology, boundary layer theory, deflection in cables, diffusion process, electromagnetism, heat transfer, and other topics. It is well known that closed form analytical solution of such problems cannot be obtained in many cases; therefore, numerical techniques such as collocation method [1, 2], B-spline interpolation [3], Hermite cubic collocation [4–6], finite difference method [7–9], nonlinear shooting method [10], geometric Hermite interpolation [11], quintic B-spline collocation method [12], polynomial and nonpolynomial spline approaches [13–15], quartic spline solution [16], cubic spline collocation method [17], and finite volume element method [18] are frequently used.

In this study, cubic Hermite collocation method (CHCM) involves cubic Hermite basis function to reduce mathematical complexity. Different linear and nonlinear differential equations are solved subject to Dirichlet, Neumann, and Robin boundary conditions using the present method. Moreover, the decoupling technique [4] used to solve elliptic problems with Neumann and Dirichlet conditions is a particular case of present technique. In this paper, linear and nonlinear boundary value problems reported in recent papers [2, 3, 7, 9, 13, 14] are solved using CHCM. It is worth mentioning that CHCM is giving better results than finite difference method, finite element method, finite volume method, B-spline method, and polynomial and nonpolynomial spline approach with fourth order of convergence.

The paper comprises five sections. Section 1 deals with general introduction of the problem. Section 2 gives brief description of cubic Hermite collocation method. In Section 3, symbolic solution of (1) is presented. Seven numerical examples are discussed in Section 4 and finally overall conclusions are given in Section 5.

2. Proposed Technique

In the present method, the domain is divided into finite elements and then orthogonal collocation method with cubic Hermite as basis function is applied within each element. The cubic Hermite interpolant of the function 𝑓 relative to partition 𝑎=𝑥₁<𝑥₂<⋯<𝑥_𝑁+1=𝑏 is a function 𝑠 that satisfies the following:(1)on each subinterval [𝑥_𝑘,𝑥_𝑘+1],  𝑠 coincides with a cubic polynomial 𝑠_𝑘(𝑥),(2)𝑠 interpolates 𝑓 and 𝑓⁽¹⁾ at 𝑥₁,𝑥₂,…,𝑥_𝑁+1,(3)𝑠 and 𝑠⁽¹⁾ are continuous on [𝑎,𝑏].

The cubic Hermite interpolant of 𝑓 and its first derivative at 𝑥=𝑥_𝑘 requires that 𝑠_𝑘𝑥_𝑘=𝑓𝑥_𝑘,𝑠_𝑘⁽¹⁾𝑥_𝑘=𝑓⁽¹⁾𝑥_𝑘. (2) Combining the continuity of 𝑠 and 𝑠⁽¹⁾ at 𝑥=𝑥_𝑘+1 with interpolation of 𝑓 and its first derivative at 𝑥=𝑥_𝑘+1, one gets 𝑠_𝑘𝑥_𝑘+1=𝑠_𝑘+1𝑥_𝑘+1=𝑓𝑥_𝑘+1,𝑠_𝑘⁽¹⁾𝑥_𝑘+1=𝑠⁽¹⁾_𝑘+1𝑥_𝑘+1=𝑓⁽¹⁾𝑥_𝑘+1. (3) Hence, 𝑠_𝑘(𝑥) is a third degree polynomial that interpolates both 𝑓 and 𝑓⁽¹⁾ at 𝑥=𝑥_𝑘 and at 𝑥=𝑥_𝑘+1. Therefore, 𝑠_𝑘(𝑥) can be written as 𝑠_𝑘(𝑥)=𝐻^𝑘₁(𝑥)𝑓𝑥_𝑘+𝐻^𝑘₂(𝑥)𝑓⁽¹⁾𝑥_𝑘+𝐻^𝑘₃(𝑥)𝑓⁽¹⁾𝑥_𝑘+1+𝐻^𝑘₄(𝑥)𝑓𝑥_𝑘+1, (4) where 𝐻^𝑘_2𝑝−1(𝑥)=𝑥−𝑥_𝑘−1ℎ_𝑘−1²23−𝑥−𝑥_𝑘−1ℎ_𝑘−1;𝑥∈𝑥_𝑘−1,𝑥_𝑘1−𝑥−𝑥_𝑘ℎ_𝑘²21−𝑥−𝑥_𝑘ℎ_𝑘;𝑥∈𝑥_𝑘,𝑥_𝑘+1𝐻0;otherwise,_𝑘2𝑝(𝑥)=−ℎ_𝑘−1𝑥−𝑥_𝑘−1ℎ_𝑘−1²1−𝑥−𝑥_𝑘−1ℎ_𝑘−1;𝑥∈𝑥_𝑘−1,𝑥_𝑘ℎ_𝑗1−𝑥−𝑥_𝑘ℎ_𝑘²𝑥−𝑥_𝑘ℎ_𝑘;𝑥∈𝑥_𝑘,𝑥_𝑘+10;otherwise. (5) Here, 𝑝=1,2 and 𝑘=1,2,…,𝑁.

The grid points, 𝑥_𝑘, are often called the “knots” of the piecewise polynomial since they are points where polynomials are “tied together.” The Hermite polynomials do not require the subsidiary condition to make first derivative continuous. This fact reduces the number of equations by (𝑁−1), where 𝑁 is the number of elements.

The global variable 𝑥 varies in the 𝑘th element, where 𝑘=1,2,…,𝑁. A new variable 𝑢=(𝑥−𝑥_𝑘)/ℎ_𝑘 is introduced in 𝑘th element in such a way that as 𝑥 varies from 𝑥_𝑘 to 𝑥_𝑘+1, 𝑢 varies from 0 to 1. Orthogonal collocation is applied on local variable 𝑢.

Approximation of function 𝑦(𝑢) in the 𝑘th element is given as [6] 𝑦(𝑢)=⁴_𝑖=1𝑎_{𝑖+2𝑘−2}𝐻_𝑖(𝑢). (6)

To apply the collocation method, one must evaluate the trial function (6) and its derivatives at two internal collocation points 𝑢=𝑢_𝑗(𝑗=1,2). These are given by 𝑦𝑢_𝑗=⁴_𝑖=1𝑎_{𝑖+2𝑘−2}𝐻_𝑖𝑢_𝑗,𝑑𝑦𝑢𝑑𝑢^𝑗=1ℎ_𝑘⁴_𝑖=1𝑎_{𝑖+2𝑘−2}𝐴_𝑗𝑖,𝑑²𝑦𝑑𝑢²𝑢_𝑗=1ℎ²_𝑘⁴_𝑖=1𝑎_{𝑖+2𝑘−2}𝐵_𝑗𝑖, (7) where the Hermite polynomials and their first and second derivatives are defined as 𝐻₁𝑢_𝑗=1+2𝑢_𝑗1−𝑢_𝑗²,𝐻₂𝑢_𝑗=𝑢_𝑗1−𝑢_𝑗²ℎ_𝑘,𝐻₃𝑢_𝑗=𝑢²_𝑗3−2𝑢_𝑗,𝐻₄𝑢_𝑗=𝑢²_𝑗𝑢_𝑗−1ℎ_𝑘,𝐴_𝑗1𝑢_𝑗=6𝑢²_𝑗−6𝑢_𝑗,𝐴_𝑗2𝑢_𝑗=1−4𝑢_𝑗+3𝑢²_𝑗ℎ_𝑘,𝐴_𝑗3𝑢_𝑗=6𝑢_𝑗−6𝑢²_𝑗,𝐴_𝑗4𝑢_𝑗=3𝑢²_𝑗−2𝑢_𝑗ℎ_𝑘,𝐵_𝑗1𝑢_𝑗=12𝑢_𝑗−6,𝐵_𝑗2𝑢_𝑗=6𝑢_𝑗−4ℎ_𝑘,𝐵_𝑗3𝑢_𝑗=6−12𝑢_𝑗,𝐵_𝑗4𝑢_𝑗=6𝑢_𝑗−2ℎ_𝑘, (8) where 𝑢_𝑗’s are the zeros of shifted orthogonal Legendre polynomial 𝑃₂^(0,0)(𝑢) with 𝑢₁=0.2113248654 and 𝑢₂=0.7886751346, as shown in Figure 1.

Figure 1 

Subdivision of mesh points on the global domain. The four coefficients, in each 𝑘 element, are estimated by using four collocation points 𝑢0, 𝑢1, 𝑢2, and 𝑢3.

3. Symbolic Solution

Equation (1) can be discretized using (7) as follows: 1ℎ²_𝑘⁴_𝑖=1𝑎_{𝑖+2𝑘−2}𝐵_𝑗𝑖+𝛽₁(𝑢)1ℎ_𝑘⁴_𝑖=1𝑎_{𝑖+2𝑘−2}𝐴_𝑗𝑖+𝛽₂(𝑢)×⁴_𝑖=1𝑎_{𝑖+2𝑘−2}𝐻_𝑗𝑖=𝑓(𝑢), (9) where, in the 1st element, 𝑎_𝑖’s vary from 1 to 4, 2nd element varies from 3 to 6, and so on. Symbolic solution is given subject to Neumann boundary conditions, for 𝑎=0 and 𝑏=1: 𝑦⁽¹⁾(0)=𝑐₀,𝑦⁽¹⁾(1)=𝑐_𝐿, (10)(11) where 𝑐₀ and 𝑐_𝐿 are finite real constants.

In the first element, 𝑢_𝑗=0⇒𝐴_𝑗1=𝐴_𝑗3=𝐴_𝑗4=0; therefore, (10) becomes 1ℎ_𝑘⁴_𝑖=1𝑎_𝑖𝐴_𝑗𝑖𝑢_𝑗=𝑐₀⟹𝑎₂=𝑐₀. (12) In the last element, 𝑢_𝑗=1⇒𝐴_𝑗1=𝐴_𝑗2=𝐴_𝑗3=0; therefore, (11) becomes 1ℎ_𝑘⁴_𝑖=1𝑎_{𝑖+2𝑁−2}𝐴_𝑗𝑖𝑢_𝑗=𝑐_𝐿⟹𝑎_2𝑁+1=𝑐_𝐿. (13) A system of (2𝑁+2) equations is obtained from (9) to (13). It includes all parameters of the system and the dependent variables at the boundaries. The support of each Hermite cubic basis function spans at most two subintervals; therefore, a band matrix is obtained with bandwidth two (Figure 2). Of these, two unknowns are found using boundary conditions and rest 2𝑁 are from discretized system of (9), using Mathematica. After substituting the appropriate values of 𝑎’s in (7), the result can be obtained for any element.

Figure 2 

Pattern of nonzero elements of banded matrix arising from CHCM for 𝑁=6, where ×’s represent nonzero elements and 0’s are represented by dots.

3.1. Error Analysis

Suppose the function 𝑓(𝑥) and its four derivatives are continuous on [𝑎,𝑏] and there is a positive constant which satisfies |𝑓⁴(𝑥)|≤𝑀 for all 𝑥∈[𝑎,𝑏]. If 𝐻(𝑥) is the cubic Hermite interpolant of 𝑓 at 𝑎 and 𝑏, then according to de Boor [19], |𝑓(𝑥)−𝐻(𝑥)|≤𝜀ℎ⁴, where ℎ=𝑏−𝑎, 𝜀=𝑀/384.

Also the placement of the collocation points plays a critical role in obtaining the 𝑂(ℎ⁴) estimate [20, 21]. For Gauss Legendre roots, as collocation points, an error estimate of 𝑂(ℎ⁴) is obtained whereas for other choices of the collocation points, only second-order accuracy is obtained.

3.2. Algorithm of the Method

Step 1. Divide the domain 0≤𝑥≤1 into a mesh 0=𝑥₁<𝑥₂<⋯<𝑥_𝑁+1=1.

Step 2. Transform the global variable 𝑥 into local variable 𝑢=(𝑥−𝑥_𝑘)/ℎ_𝑘.

Step 3. Approximate the solution at 𝑢=𝑢_𝑗.

Step 4. Obtain the trial function 𝑦(𝑢_𝑗) and the derivatives (𝑑𝑦/𝑑𝑢)(𝑢_𝑗), (𝑑²𝑦/𝑑𝑢²)(𝑢_𝑗).

Step 5. Carry out discretization of the model using Step 4.

Step 6. Evaluate Step 5 at collocation points 𝑢₁=0.2113248654 and 𝑢₂=0.7886751346.

Step 7. The obtained system in Step 6 is solved using any software.

4. Numerical Examples and Discussion

In this section, seven examples demonstrate the efficiency and accuracy of the method. Following formulae are used for estimation of error in this study.

Rate of convergence of CHCM is calculated using 𝜌=(ln(𝐿^∞(𝑁2)/𝐿^∞(𝑁1)))/(ln(𝑁2/𝑁1)).

Relative error is obtained by (𝑦_ex−𝑦_nm)/𝑦_ex.

Max norm is found by 𝐿^∞=max^𝑛_𝑖=1|𝑦(𝑢_𝑖)−𝑦(𝑢_𝑖)|.

Example 1. Solve (1) for 𝛼₁(𝑥)=0,𝛼₂(𝑥)=−𝐻²,𝑓(𝑥)=0, (14) subject to Robin’s boundary conditions 𝑦(0)=1,𝑦⁽¹⁾(1)=0, (15) which has exact solution [9], for planer geometry ()𝑦=cosh𝐻1−𝑥.cosh𝐻 (16)
The problem is solved for different values of dimensionless parameter 𝐻 by taking 10 to 40 elements. The exact and numeric results are plotted in Figure 3. The results reported by [9] for 10 elements using finite difference method are matching with exact solution up to 3 decimal places, whereas using CHCM the results are matching up to 9 decimal places. This shows the superiority of cubic Hermite collocation method over the finite difference method. Relative error between CHCM and exact values is presented in Table 1.

Table 1 

Relative error for different values of parameter 𝐻 for Example 1.

Figure 3 

Temperature profiles in a rectangular fin for dimensionless heat transfer coefficient 𝐻.

Example 2. Solve (1) for 𝛼₁(𝑥)=−1,𝛼₂(𝑥)=0,𝑓(𝑥)=−𝑒^𝑥−1−1, (17) subject to Dirichlet boundary conditions 𝑦(0)=𝑦(1)=0,0<𝑥<1, (18) which has exact solution [3, 7] 𝑦(𝑥)=𝑥1−𝑒^𝑥−1. (19)
The exact and numeric results are plotted in Figure 4. The CHCM results are matching up to 11 decimal places with the exact ones in Table 2. On comparing present results with the results of [3, 7], shown in Table 3, a big difference of errors between CHCM with finite difference method, finite element method, finite volume method, and B-spline method is observed. This indicates the supremacy of the present method.

Table 2 

Relative error between CHCM and exact values for Example 2.

Table 3 

Max norm of errors for five methods with respect to exact solution.

Figure 4 

Comparison of CHCM with exact result for Example 2.

Example 3. Solve (1) for 𝛼₁(𝑥)=0,𝛼₂(𝑥)=−𝜙²,𝑓(𝑥)=0, (20) subject to Dirichlet boundary conditions 𝑦(0)=0.1,𝑦(1)=1, (21) which has exact solution [2], for 𝜙²=10 and ()𝑦=𝑔sinh𝜙𝑥+𝑓sinh𝜙1−𝑥.sinh𝜙 (22)
The exact and numeric results are plotted in Figure 5. The results are matching up to 9 decimal places as shown in Table 4. From Table 5, order of convergence is found to be 4.

Table 4 

Relative error between CHCM and exact values for Example 3.

Table 5 

Rate of convergence of CHCM for Example 3.

Figure 5 

Comparison of CHCM with exact result for Example 3.

Example 4. Solve (1) for 𝛼₁(𝑥)=0,𝛼₂(𝑥)=1,𝑓(𝑥)=−1, (23) subject to Neumann boundary conditions 𝑦⁽¹⁾(0)=1−cos(1)sin(1),𝑦⁽¹⁾(1)(1)=−1−cos(1),sin (24) which has exact solution [13, 14] as 𝑦(𝑥)(1)=cos𝑥+1−cossin(1)sin𝑥−1. (25)

Example 5. Solve (1) for 𝛼₁(𝑥)=0,𝛼₂(𝑥)𝑓(𝑥)==𝑥,3−𝑥−𝑥²+𝑥³sin𝑥+4cos𝑥, (26) subject to Neumann boundary conditions 𝑦⁽¹⁾(0)=−1,𝑦⁽¹⁾(1)=2sin(1), (27) which has exact solution [13, 14] as 𝑦(𝑥)=𝑥²−1sin𝑥. (28)
For different number of elements 𝑁=8, 16, 32, 64, and 128, Examples 4 and 5 are solved using the present method. Maximum absolute errors of the numerical solutions are calculated and compared with those reported by [13, 14], for different methods, in Tables 6 and 7. The present approach is giving much more accurate results than the others. From Table 8, the order of convergence is found to be 4.

Table 6 

Maximum absolute errors obtained by different methods in Example 4.

Table 7 

Maximum absolute errors obtained by different methods in Example 5.

Table 8 

Rate of convergence of CHCM for Examples 4 and 5.

Example 6. Consider a nonlinear problem, that is, solving (1) for 𝛼₁(𝑥)=0,𝛼₂(𝑥)𝑓(𝑥)=−𝑦,=2𝜋²cos(2𝜋𝑥)−sin⁴(𝜋𝑥), (29) subject to Neumann boundary conditions 𝑦⁽¹⁾(0)=0,𝑦⁽¹⁾(1)=0, (30) which has exact solution [13] 𝑦(𝑥)=sin²(𝜋𝑥). (31)
The exact and CHCM results are showing good agreement in Figure 6 for 𝑁=50.

Figure 6 

Comparison of CHCM with exact result for Example 6.

Example 7. Consider a nonlinear problem, that is, solving (1) for 𝛼₁(𝑥)=0,𝛼₂(𝑥)=𝑒^−2𝑦𝑦,𝑓(𝑥)=0, (32) subject to Neumann boundary conditions 𝑦⁽¹⁾(0)=1,𝑦⁽¹⁾(1)=12, (33) which has exact solution [13] 𝑦(𝑥)=ln(1+𝑥). (34)
An excellent matching is found between the exact and CHCM results in Figure 7 for 𝑁=50. In Table 9, the order of convergence is again found to be 4 for Examples 6 and 7.

Table 9 

Rate of convergence of CHCM for Examples 6 and 7.

Figure 7 

Comparison of CHCM with exact result for Example 7.

5. Conclusion

In this paper, cubic Hermite collocation method is tested for seven problems. The numerical results obtained are quite satisfactory and comparable with the existing solution available in the literature. The superiority over the finite difference method, finite element method, finite volume method, B-spline method, and polynomial and nonpolynomial spline approach shows the strength of this method. The convergence of the CHCM technique is of order 4.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by NBHM, Mumbai, and UGC, New Delhi, India, in the form of Research Projects 2/48(14)/2009/R&D-II/2806 and 41-786/2012(SR), respectively. The JRF provided to Mr. I. A. Ganaie by NBHM, Mumbai, is thankfully acknowledged.