Solution of two point boundary value problems, a numerical approach: parametric difference method

There has been a great deal of interest in developing techniques for the accurate numerical solution of two-point boundary value problems. In this article we have considered a parametric difference method for the numerical solution of the boundary value problems with the form

u″(x)=g(x,u,u′),a<x<b$$\begin{array}{} \displaystyle u''(x)=g(x,u,u'),\quad a \lt x \lt b \end{array}$$(1)

subject to the boundary conditions

u′(a)=αandp1u(b)+q1u′(b)=β,$$\begin{array}{} \displaystyle u'(a)=\alpha \quad \text{and} \quad p_1u(b)+q_1u'(b)=\beta , \nonumber \end{array}$$

or

p0u(a)+q0u′(a)=α0andu′(b)=β0,$$\begin{array}{} \displaystyle p_0u(a)+q_0u'(a)=\alpha_0 \quad \text{and}\quad u'(b)=\beta_0 , \nonumber \end{array}$$

where α, β, p₀ etc. are real constants and known. We assume that g is continuous function in ([a, b] × ℜ², ℜ).

In many branches of engineering, humanity and sciences ordinary differential equations occur as mathematical models. In general, these equations have solutions that can be expressed in restricted form only. So it is essential to consider approximate solutions by means of numerical techniques. To solve these problems numerically, we have many accurate numerical methods for instance, shooting-projection [1], collocation [2], finite difference methods [3] available in the literature. However the relatively less advanced method, the method of parameter differentiation is known in the literature and some historical development can be found in [4]. The existence and uniqueness of the solution for the problem (1) are assumed. The specific assumption on g(x, u, u′) to ensure the existence and uniqueness will not be considered in this article however I may refer literature in [5, 6, 7, 8].

In this article we shall develop a hybrid parametric difference method for the approximate numerical solution of problems (1). We have performed a numerical experiment to demonstrate the effectiveness of the method. We have achieved accuracy and bound on the error between the analytical and numerical approximation solution of the problem. To the best of our knowledge, no method similar to the proposed method for the numerical approximate solution has been discussed in the literature to date.

We have presented our work in this article as follows. In the next section we will discuss the parametric difference method and in Section 3 we will derive our propose method. In Section 4 and 5, respectively, we have estimated truncation error and discussed convergence under appropriate condition. The applications of the proposed method to the model problems and numerical results have been produced to show the efficiency in Section 6. Discussion and conclusion on the performance of the method are presented in Section 7.

Following the ideas in [4], re-formulate the problem (1) either by introducing a physical parameter or identify a physical parameter that appears either in the differential equation or in the boundary conditions as independent variable. Let us introduce λ as physical parameter in differential equation (1), so we have obtained,

u″(x)=g(x,u,u′)≡f(x,λ,u,u′),a<x<b$$\begin{array}{} \displaystyle u''(x)=g(x,u,u')\equiv f(x,\lambda,u,u'),\quad a \lt x \lt b \end{array}$$(2)

and boundary conditions are

u′(a)=αandp1u(b)+q1u′(b)=β,$$\begin{array}{} \displaystyle u'(a)=\alpha \quad \text{and} \quad p_1u(b)+q_1u'(b)=\beta , \nonumber \end{array}$$

or

p0u(a)+q0u′(a)=α0andu′(b)=β0,$$\begin{array}{} \displaystyle p_0u(a)+q_0u'(a)=\alpha_0 \quad \text{and}\quad u'(b)=\beta_0 , \nonumber \end{array}$$

Now let introduce new variable y(x) = u̇ = dudλ$\begin{array}{} \displaystyle \frac{du}{d\lambda} \end{array}$, a differentiation of solution of problem (1) with respect to λ. Differentiate problem (2) with respect to λ. Thus we have reduced the problem (2) into following system of equations,

y″=F(∂f∂λ,y,y′),u˙=y(x)$$\begin{array}{} \displaystyle y''=F(\frac{\partial f}{\partial \lambda},y,y'), \nonumber \\ \displaystyle\dot{u}=y(x) \end{array}$$(3)

subject to the boundary conditions,

y′(a)=0andp1y(b)+q1y′(b)=0,$$\begin{array}{} \displaystyle y'(a)=0 \quad \text{and} \quad p_1y(b)+q_1y'(b)=0 , \end{array}$$

or

p0y(a)+q0y′(a)=0andy′(b)=0,$$\begin{array}{} \displaystyle p_0y(a)+q_0y'(a)=0 \quad \text{and}\quad y'(b)=0 , \end{array}$$

and y(x) is known for some λ.

We wish to determine the approximate numerical solution u(x) of the problem (2) for different values of λ and λ₀ ≤ λ ≤ λ_M. We define M – 1 numbers of nodal points in [λ₀, λ_M], as λ₀ < λ₁ < λ₂ < …… < λ_M using a uniform step length H such that λ_i = λ₀ + jH, j = 0, 1, .., M. Also, we define N – 1 numbers of nodal points in [a, b], the domain in which the solution of the problem (1) is desired, as a ≤ x₀ < x₁ < x₂ < …… < x_N = b using a uniform step length h such that x_i = x₀ + ih, i = 0, 1, .., N. To simplify the expression we denote the numerical approximation of u(x) at the node x = x_i as u_i. Let us denote F_i as the approximation of the theoretical value of the source function F( ∂f∂λ$\begin{array}{} \displaystyle \frac{\partial f}{\partial \lambda} \end{array}$, y(x), y′(x)) at node x = x_i, i = 0, 1, 2, ……, N. Let write a system of equations (3) at nodal point x_i in notational form as,

yi″=Fiu˙i=yi$$\begin{array}{} \displaystyle y''_i=F_i \\ \dot{u}_i=y_i \end{array}$$

subject to the boundary conditions

y0′=0andp1yN+q1yN′=0,$$\begin{array}{} \displaystyle y'_0=0 \quad \text{and} \quad p_1y_N+q_1y'_N=0 , \nonumber \end{array}$$

or

p0y0+q0y0′=0andy′(b)=0,$$\begin{array}{} \displaystyle p_0y_0+q_0y'_0=0 \quad \text{and}\quad y'(b)=0 , \nonumber \end{array}$$

Following the idea in [9], we propose our difference method for the numerical solution of problem (3) as follows,

yi−<!-- - -->yi−<!-- - -->1=hyi−<!-- - -->1′+h26(2Fi−<!-- - -->1+Fi),i=1,2,..,N,yi′−<!-- - -->yi−<!-- - -->1′=h2(Fi+Fi−<!-- - -->1),i=1,2,..,N,uj−<!-- - -->1−<!-- - -->uj=−<!-- - -->H(yi−<!-- - -->1+((λ<!-- ? -->0−<!-- - -->a)+((j−<!-- - -->1)H−<!-- - -->(i−<!-- - -->1)h))yi−<!-- - -->1′),j=1,2,..,M.$$\begin{array}{} \displaystyle \qquad\qquad\qquad\qquad\qquad\quad~\, y_i-y_{i-1}=hy'_{i-1}+\frac{h^2}{6}(2F_{i-1}+F_i), \qquad i=1,2,..,N, \nonumber \\ \displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad~~~\,\! y'_i-y'_{i-1}=\frac{h}{2}(F_i+F_{i-1}), \qquad i=1,2,..,N, \nonumber \\ \displaystyle u_{j-1}-u_j= -H(y_{i-1}+((\lambda_0-a)+((j-1)H-(i-1)h))y'_{i-1}), \qquad j=1,2,..,M. \end{array}$$(4)

which is a system of linear/ nonlinear equations depending on forcing function g(x, u, u′). We solve system of equations (4) by applying an appropriate iterative method.

In this section we outline the derivation and development of the proposed finite difference method. We approximate first problem in (3) as a linear combination of y(x), y′(x) and F( ∂f∂λ$\begin{array}{} \displaystyle \frac{\partial f}{\partial \lambda} \end{array}$, y, y′),

a0yi−1+a1yi+a2hyi−1′+a3hyi′+h2(b0Fi−1+b1Fi)=0$$\begin{array}{} \displaystyle a_0y_{i-1}+a_1y_i+a_2hy'_{i-1}+a_3hy'_i+h^2(b_0F_{i-1}+b_1F_i)=0 \end{array}$$(5)

where a₀, ⋯, b₁ are constant. To determine these constants, we expand each term in the above expression in a Taylor series about a node x_i and compare the coefficients of h^p, p = 0, ⋯, 3 in the both sides of the expression. So we obtained a system of linear equations in a₀, ⋯, b₁. Solving the system of equations, we got

(a0,a1,a2,a3,b0,b1)=(−1,1,−1,0,−13,−16)$$\begin{array}{} \displaystyle (a_0,a_1,a_2,a_3,b_0,b_1)=(-1,1,-1,0,-\frac{1}{3},-\frac{1}{6}) \end{array}$$(6)

Substituting these constants from (6) in (5), we obtained

yi−yi−1−hyi−1′−h26(2Fi−1+Fi)=0$$\begin{array}{} \displaystyle y_i-y_{i-1}-hy'_{i-1}-\frac{h^2}{6}(2F_{i-1}+F_i)=0 \end{array}$$(7)

Following the similar idea as above, we derive following equations,

yi′−yi−1′=h2(Fi+Fi−1),uj−1−uj=−H(yi−1+((λ0−a)+((j−1)H−(i−1)h))yi−1′)$$\begin{array}{} \displaystyle y'_i-y'_{i-1}=\frac{h}{2}(F_i+F_{i-1}),\\ \displaystyle\quad\! u_{j-1}-u_j= -H(y_{i-1}+((\lambda_0-a)+((j-1)H-(i-1)h))y'_{i-1}) \end{array}$$(8)

The local truncation error in the difference method (4) using the exact arithmetic may be calculated by using Taylor series expansion at the nodal points x = x_i, i = 1, 2, .., N. Thus the truncation error T_i in the first equation of the method (4) may be written as:

Ti−1=−h424yi(4),1≤i≤N.$$\begin{array}{} \displaystyle T_{i-1}=-\frac{h^4}{24}y^{(4)}_i ,\qquad1\leq i\leq N. \end{array}$$(9)

Similarly we can obtain truncation error in second equation of method (4) at x = x_i, i = 1, 2, .., N as given below,

Ti′=−h33yi(4),1≤i≤N.$$\begin{array}{} \displaystyle T'_i=-\frac{h^3}{3}y^{(4)}_i , \quad 1\leq i \leq N. \end{array}$$(10)

The truncation error in final equation of method (4) at λ_j–1, j = 1, 2, .., M and x_i–1, i = 1, 2, .., N is,

T∗=H2u¨j−1−((λ0−a)+((j−1)H−(i−1)h))22yi−1″.$$\begin{array}{} \displaystyle T_{*}=\frac{H}{2}\ddot{u}_{j-1}-\frac{((\lambda_0-a)+((j-1)H-(i-1)h))^2}{2}y''_{i-1}. \end{array}$$(11)

Thus we have obtained a truncation error at each node of order O(H + h²) if λ₀ = a.

Let us assume the following conditions [10],

u(x) and u′(x) are finite

Under above conditions problem (1) will posses unique solution [7]. However the forcing function F( ∂f∂λ$\begin{array}{} \displaystyle \frac{\partial f}{\partial \lambda} \end{array}$, y, y′) in (3) may be written as,

F(∂f∂λ,y,y′)≡∂f∂uy+∂f∂u′y′+∂f∂λ$$\begin{array}{} \displaystyle F(\frac{\partial f}{\partial \lambda},y,y')\equiv \frac{\partial f}{\partial u}y+ \frac{\partial f}{\partial u'}y'+\frac{\partial f}{\partial \lambda} \end{array}$$(12)

Let us assume Mi=(∂f∂y)i,Mi′=(∂f∂y′)i$\begin{array}{} \displaystyle M_i=(\frac{\partial f}{\partial y})_i, M'_i=(\frac{\partial f}{\partial y'})_i \end{array}$ and using (12), we can write difference method (4) as follows,

−(1+h23Mi−1)yi−1+(1−h26Mi)yi−h26Mi′yi′−h(1+h3Mi−1′)yi−1′=h26(2(∂f∂λ)i−1+(∂f∂λ)i),−h2Mi−1yi−1−h2Miyi+(1−h2Mi′)yi′−(1−h2Mi−1′)yi−1′=h2((∂f∂λ)i−1+(∂f∂λ)i),i=1,2,..,N$$\begin{array}{} \displaystyle -(1+\frac{h^2}{3}M_{i-1})y_{i-1}+(1-\frac{h^2}{6}M_i)y_i-\frac{h^2}{6}M'_iy'_i-h(1+\frac{h}{3}M'_{i-1})y'_{i-1} =\frac{h^2}{6}(2(\frac{\partial f}{\partial \lambda})_{i-1}+(\frac{\partial f}{\partial \lambda})_i), \\ \displaystyle -\frac{h}{2}M_{i-1}y_{i-1}-\frac{h}{2}M_iy_i+(1-\frac{h}{2}M'_i)y'_i-(1-\frac{h}{2}M'_{i-1})y'_{i-1} =\frac{h}{2}((\frac{\partial f}{\partial \lambda})_{i-1}+(\frac{\partial f}{\partial \lambda})_i), \qquad i=1,2,..,N\qquad\quad \end{array}$$(13)

Let Y_i and y_i be respectively an exact solution and approximate solution of the system of equations (3). Let us define ε_i = y_i – Y_i, an error between the exact analytical and numerical approximation solution of the problem (3) and similarly we can define εi′=yi′−Yi′.$\begin{array}{} \displaystyle \varepsilon'_i=y'_i-Y'_i. \end{array}$ Let us write (13) in matrix form

Dε+T=0,$$\begin{array}{} \displaystyle \mathbf{D}\boldsymbol{\varepsilon}+\mathbf{T} = \boldsymbol{0}, \end{array}$$(14)

where matrix D is of order 2(N – 1) × 2(N – 1), moreover we have partitioned the D in the following manner,

D=A11⋮A12⋯⋮⋯A21⋮A22$$\begin{array}{} \displaystyle \mathbf{D}=\left( \begin{array}{ccc} \mathbf{A}_{11} &\vdots & \mathbf{A}_{12} \\ \cdots &\vdots &\cdots \\ \mathbf{A}_{21}&\vdots & \mathbf{A}_{22} \\ \end{array} \right) \end{array}$$

where each matrix A is of order (N – 1) × (N – 1) and these matrices are,

A11=−1−h23M01−h26M10−1−h23M11−h26M2⋱⋱−1−h23MN−21−h26MN−10−1−h23MN−1N−1×N−1$$\begin{array}{} \displaystyle \mathbf{A}_{11}=\left( \begin{array}{cccc} -1-\frac{h^2}{3}M_0 & 1-\frac{h^2}{6}M_1 & & 0 \\ & -1-\frac{h^2}{3}M_1 & 1-\frac{h^2}{6}M_2 & \\ & \ddots & \ddots & \\ & & & \\ & & -1-\frac{h^2}{3}M_{N-2} & 1-\frac{h^2}{6}M_{N-1} \\ 0 & & & -1-\frac{h^2}{3}M_{N-1}\\ \end{array} \right)_{N-1 \times N-1} \end{array}$$

an upper diagonal matrix with bandwidth of 1,

A12=−h26M1′0−h−h23M1′−h26M2′⋱⋱−h−h23MN−2′−h26MN−1′0−h−h23MN−1′−h26MN′N−1×N−1$$\begin{array}{} \displaystyle \mathbf{A}_{12}=\left( \begin{array}{cccc} -\frac{h^2}{6}M'_1 & & & 0 \\ -h-\frac{h^2}{3}M'_1 & -\frac{h^2}{6}M'_2 & & \\ \ddots & \ddots & &\\ & & & \\ & -h-\frac{h^2}{3}M'_{N-2} & -\frac{h^2}{6}M'_{N-1}& \\ 0 & &-h-\frac{h^2}{3}M'_{N-1} & -\frac{h^2}{6}M'_N\\ \end{array} \right)_{N-1 \times N-1} \end{array}$$

a lower diagonal matrix with bandwidth of 1,

A21=−h2M0−h2M10−h2M1−h2M2⋱⋱−h2MN−2−h2MN−10−h2MN−1N−1×N−1$$\begin{array}{} \displaystyle \mathbf{A}_{21}=\left( \begin{array}{cccc} -\frac{h}{2}M_0 & -\frac{h}{2}M_1 & & 0 \\ & -\frac{h}{2}M_1 & -\frac{h}{2}M_2 & \\ & \ddots & \ddots & \\ & & & \\ & & -\frac{h}{2}M_{N-2} & -\frac{h}{2}M_{N-1} \\ 0 & & & -\frac{h}{2}M_{N-1}\\ \end{array} \right)_{N-1 \times N-1} \end{array}$$

an upper diagonal matrix with bandwidth of 1 and finally,

A22=1−h2M1′0−1−h2M1′1−h2M2′⋱⋱−1−h2MN−2′1−h2MN−1′0−1−h2MN−1′1−h2MN′N−1×N−1$$\begin{array}{} \displaystyle \mathbf{A}_{22}=\left( \begin{array}{cccc} 1-\frac{h}{2}M'_1 & & & 0 \\ -1-\frac{h}{2}M'_1 & 1-\frac{h}{2}M'_2 & & \\ \ddots & \ddots & &\\ & & & \\ & -1-\frac{h}{2}M'_{N-2} & 1-\frac{h}{2}M'_{N-1}& \\ 0 & &-1-\frac{h}{2}M'_{N-1} & 1-\frac{h}{2}M'_N\\ \end{array} \right)_{N-1 \times N-1} \end{array}$$

a lower diagonal matrix with bandwidth of 1,

ε=(ε0,ε1,..,εN−1,ε1′,..,εN′)t,$$\begin{array}{} \displaystyle \boldsymbol{\varepsilon}=(\varepsilon_0,\varepsilon_1,..,\varepsilon_{N-1},\varepsilon'_1,..,\varepsilon'_N)^t , \end{array}$$

and

T=(T0,T1,..,TN−1,T1′,..,TN′)t.$$\begin{array}{} \displaystyle \mathbf{T}=(T_0,T_1,..,T_{N-1},T'_1,..,T'_N)^t . \end{array}$$

are 1 × 2(N – 1) matrices.

Let us assume that either Mi′<0orM′≠2hfori=1,2,..,N.$\begin{array}{} \displaystyle M'_i \lt 0 ~\text{or}~ M'\neq\frac{2}{h} ~\text{for}~ i=1,2,..,N. \end{array}$ These either assumptions ensure that diagonal block matrices A₁₁ and A₂₂ are invertible. Moreover A11−1<0$\begin{array}{} \displaystyle \mathbf{A}^{-1}_{11} \lt 0 \end{array}$ a negative definite but A22−1>0$\begin{array}{} \displaystyle \mathbf{A}^{-1}_{22} \gt 0 \end{array}$ a positive definite if |Mi′|<2h.$\begin{array}{} \displaystyle |M'_i| \lt \frac{2}{h}. \end{array}$ Let us assume

|A12A22−1A21A11−1|<1$$\begin{array}{} \displaystyle |\mathbf{A}_{12}\mathbf{A}^{-1}_{22}\mathbf{A}_{21}\mathbf{A}^{-1}_{11}| \lt 1 \end{array}$$(15)

hence block matrix D is invertible [11].

Let us define,

a2∗=∥A12A22−1∥∞anda∗1=∥A21A11−1∥∞.$$\begin{array}{} \displaystyle a^*_2=\|\mathbf{A}_{12}\mathbf{A}^{-1}_{22}\|_\infty \quad\text{and} \quad a_{*1}=\|\mathbf{A}_{21}\mathbf{A}^{-1}_{11}\|_\infty . \end{array}$$

and

Si∗=−h26(2Mi−1+Mi),1≤i≤N−1−1−h23Mi−1,i=N$$\begin{array}{} \displaystyle S^*_i=\begin{cases} -\frac{h^2}{6}(2M_{i-1}+M_i) ,& \quad 1\leq i\leq N-1 \\ -1-\frac{h^2}{3}M_{i-1} ,& \quad i= N \end{cases} \end{array}$$

and

S∗i=1−h2Mi′,i=1−h2(Mi−1′+Mi′),2≤i≤N$$\begin{array}{} \displaystyle S_{*i}=\begin{cases} 1-\frac{h}{2}M'_i ,& \quad i=1 \\ -\frac{h}{2}(M'_{i-1}+M'_i),& \quad 2\leq i \leq N \end{cases} \end{array}$$

sum of the elements in a row of the matrices A₁₁ and A₂₂ respectively. Let Smin∗=min{Si∗},1≤i≤N$\begin{array}{} \displaystyle S^*_{min}=min\{S^*_i\}, 1\leq i\leq N \end{array}$ and S_∗max = max{S_∗i}, 1 ≤ i ≤ N. Thus

∥A11−1∥∞≤1|Smin∗|and∥A22−1∥∞≤1|S∗max|.$$\begin{array}{} \displaystyle \|\mathbf{A}^{-1}_{11}\|_\infty\leq \frac{1}{|S^*_{min}|}\quad \text{and}\quad \|\mathbf{A}^{-1}_{22}\|_\infty\leq \frac{1}{|S_{*max}|}. \end{array}$$

So, we have

1Smin∗=max{1|Smin∗|,1|S∗max|}$$\begin{array}{} \displaystyle \frac{1}{S^*_{min}}= \text{max}\{\frac{1}{|S^*_{min}|},\frac{1}{|S_{*max}|}\} \end{array}$$(16)

But we know in [11] that

∥D−1∥∞≤max{∥A11−1∥∞,∥A22−1∥∞}(1+a2∗)(1+a∗1)(1+a2∗)+(1+a∗1)−(1+a2∗)(1+a∗1)$$\begin{array}{} \displaystyle \|\mathbf{D}^{-1}\|_\infty\leq\frac{\text{max}\{\|\mathbf{A}^{-1}_{11}\|_\infty,\|\mathbf{A}^{-1}_{22}\|_\infty\}(1+a^*_2)(1+a_{*1})}{(1+a^*_2)+(1+a_{*1})-(1+a^*_2)(1+a_{*1})} \end{array}$$(17)

Each term in (17) is well defined, so from (14)

∥ε∥∞≤max{∥A11−1∥∞,∥A22−1∥∞}(1+a2∗)(1+a∗1)(1+a2∗)+(1+a∗1)−(1+a2∗)(1+a∗1)∥T∥∞$$\begin{array}{} \displaystyle \|\boldsymbol{\varepsilon}\|_\infty\leq \frac{\text{max}\{\|\mathbf{A}^{-1}_{11}\|_\infty,\|\mathbf{A}^{-1}_{22}\|_\infty\}(1+a^*_2)(1+a_{*1})}{(1+a^*_2)+(1+a_{*1})-(1+a^*_2)(1+a_{*1})}\|\mathbf{T}\|_\infty \end{array}$$(18)

It follows from (13), (14) and (18) that ∥ε∥ → 0 as h → 0. Thus we conclude that method (4) will converge.

In this section we have considered model problems to perform the numerical experiment. In these model problems, we consider different number of noddle points for both x_i and λ_j. In computation of maximum absolute error MAE between the analytical solution u(λ_j) and computed numerical solution u_j of the problem, we have used the following formula,

MAE=max1≤j≤Mu(λj)−uj.$$\begin{array}{} \displaystyle MAE= \max_{\substack {1\leq{j}\leq {M}}}{\left\vert u(\lambda_j)-u_j\right\vert }. \end{array}$$

We have respectively applied Gauss-Seidel and Newton-Raphson method to solve the system of linear and nonlinear equations those arise from the method (4). The solutions are computed on different values of N and M. The iteration is continued until either the maximum difference between two successive iterates is less than 10^–8 or the number of iterations reached 10³. All computations were performed on a Windows 2007 Ultimate operating system in the GNU FORTRAN environment version 99 compiler (2.95 of gcc) on Intel Core i3-2330M, 2.20 Ghz PC.

A parametric finite difference method to find the numerical solution of two point boundary value problems with uniform mesh has been developed and discussed. The method discretizes the problem (1) at the discrete nodal points and is transformed into a system of algebraic equations given by (4). The propose method produces an approximate numerical solution of the model problems with the uniform step size. We have discussed the convergence of the proposed method, but we have not estimated exact rate of convergence. The numerical results for the model problems showed that the proposed method is convergent for large values of N. The idea presented in this article leads to the possibility to develop efficient computational method for the numerical solution of higher order boundary value problems. Works in these directions are in progress.