Research Article Volume 2 Issue 6
Associate Professor, Dyal Singh College (University of Delhi), India
Correspondence: Pramod Kumar Pandey, Associate Professor, Dyal Singh College (University of Delhi), Lodhi Road, New Delhi 110003, India, Tel 9910206270,
Received: December 21, 2017 | Published: November 28, 2018
Citation: Pandey PK. Solving numerically a sixth order differential equation as coupled finite difference equations approach. Open Access J Sci. 2018;2(6):381-385. DOI: 10.15406/oajs.2018.02.00115
In this article we have considered sixth order boundary value problem. We have proposed the finite difference method for solving problem as coupled equations. We get numerical approximation of third derivative of solution of the problem as a byproduct of the proposed method. We have established theoretically the convergence of the proposed method under appropriate conditions. We have tested the proposed method on model problems for the numerical result. Experimental result approves the theoretical result.
Keywords: boundary value problem, finite difference method, quadratic order convergence, sixth order differential equation, splitting method
In the present article we consider a technique for the numerical solution of the sixth order boundary value problems of the following form:
u6(x)= f(x,u), a<x<b,u6(x)= f(x,u), a<x<b, (1.1)
Subject to the boundary conditions
u(a)=α1,un(a)=α2,u4(a)=α3,u(b)−β1,u3(b)=β2 and u4(b)=β2u(a)=α1,un(a)=α2,u4(a)=α3,u(b)−β1,u3(b)=β2 and u4(b)=β2
Where α1,α2,α3 β1,β2 and α1,α2,α3 β1,β2 and are real constant.
The occurrence of sixth order deferential equation and corresponding boundary value problems in detail studied and discussed in.1 Let us assume f(x, u)f(x, u) be smooth to ensure the existence and uniqueness of the solution of the problem (1.1). However the detail analytical concepts on existence, uniqueness of solution of problem (1.1) discussed in.2 The emphasis in present article will be on the development of a numerical technique for the numerical solution of the problem (1.1). Several numerical methods for solving the sixth order boundary value problem have been reported in literature. Some of these reported methods such as finite difference3,4 finite element5,6 spline solution7,8 differential transformation and adomian decomposition,9 reproducing kernel,10 modified decomposition11 and homotopy perturbation12 and other references therein have been developed for numerical solution of sixth order boundary value problems. Some progress has been made in recent years in development of numerical technique for the solution of third order boundary value problems. These techniques are very satisfactory and yield a highly accurate solution.13 Hence, the purpose of this article is to incorporate this development in developing numerical technique for numerical solution of sixth order boundary value problems (1.1). We have developed finite difference method for numerical solution of sixth order boundary value problem by splitting method, a system of boundary value problems. We hope that others may and the proposed method an improvement to those existing finite difference methods for sixth order boundary value problems. We have presented our work in this article as follows: In Section 2 the finite difference method, in Section 3 we have derived a finite difference method. In Section 4, we have discussed convergence of the proposed method under appropriate condition. The experiment with the proposed method on model problems and numerical results in Section 5. A discussion and conclusion on the overall performance of the proposed method are presented in Section 6.
The difference method
Let us assume u(x)u(x) is solution of the considered sixth order differential equation and following third order boundary value problem,
u(3)(x)=v(x), a<x<bu(3)(x)=v(x), a<x<b (2.1)
And the boundary conditions are
u(a)=α1,un(a)=α2 and u(b)=β1u(a)=α1,un(a)=α2 and u(b)=β1
So problem (2) represents a splitting of system (1) into two coupled problems. Where v(x)v(x) is some regular differentiable function in [a; b]. Finally we have following problem, a splitting of problem (1.1) into coupled third order boundary value problem,
u(3)(x)=f(x,u), a<x<bu(3)(x)=f(x,u), a<x<b (2.2)
And the boundary conditions are
υ(b)=β2,υ'(a)=α3 and υ'(b)=β2
If u"(b) is known instead of u3(b) , so we will have to approximate the value of v (b). Thus the six order boundary value problem (1.1) has been reduced to a system of cubic order boundary value problems (2.1)-(2.2). We partition the interval [a, b] on which the solution of problem (1.1) is desired to introduce finite number of nodal points into a subintervals a ≤ x0 < x1 x2<......<xN+1≤b by using uniform step length h such that xi = a + ih; i =0,1,2,....,N+1 . We wish to determine the numerical solution u(x) of the problem (1.1) at these nodal points xi . Let ui denote the numerical approximation of u(x) at these node x= xi, i=1,2,....N .Let fi denote the approximate value of the source function f(x, u(x)) at these node x= xi, i=0,1,2,....N+1 . Thus the boundary value problem (1.1) replaced by the system of boundary value problems (2.1)-(2.2) may be written at this node x= xi, i=0,....N+1 as,
u(3)i=ui (2.3)
u(3)i=fi
In this case, we have two cubic order differential equations with different boundary conditions. Following the ideas in,13,14 we propose our finite difference method for a numerical solution of problem (2.3),
−ui−1 + 2ui−ui+1 =−h2u"i−1+h312(−17ui+5ui+1), i=1 (2.4)
−ui−1 + 2ui−ui+1 =−h2u"i−2+h312(−47ui+23ui+1), i=2
ui−2 − 3ui−1+3ui−ui−1 =h33(7fi+ui+1), 3≤i≤N
5ui −8ui+1−3ui+2 =2h2u'i−1+h33(7f+4fi+1), i=1 (2.5)
−υi−1_3υi−3υ+1+υi+2=h32(fi+4fi+1), 2≤i≤N−1
−υi−1_3υi−3υ+1+υi+2=h32(fi+4fi+1), 2≤i≤N−1
If the forcing function f(x, u) in problem (1.1) is linear then the system of equations (2.4) will be linear otherwise we will obtain nonlinear system of equations.
Derivation of the difference method
In this section we outline the derivation of the proposed method, we have followed the same approach as given in.13 Let us write a linear combination of solution u(x) and v(x) at nodes xx+1,xi
a2ui+1+a1ui−1+aoui+h2b1u"i−1+h3(b2u'i+1+bou'1)=0 (3.1)
Where a0,a1,a2,b0,b1 and b2 are constants to be determined. Expanding each term on the left hand side of (3.1) in Taylor series about the point xi and using method of undetermined coefficients, we get
(a0,a1,a2,b0,b1,b2)=(2,−1,−1,1712,−512,1) (3.2)
Thus from (3.1)-(3.2), we have
−(vi−1+2vi−vi+1)+h2u"i−1+h312(17fi−5fi+1)+Tui =0 (3.3)
Where Tui is local discretization error and equal to −43h3120u(5)i . Similarly we can derive the following equations
−(vi−1+2vi−vi+1)+h2u"i−1+h312(17fi−5fi+1)+Tui =0 (3.4)
ui−2−3ui−1+3ui−ui+1=h32(−3ui+ui+1)+Tui, 3≤i≤N
Where local discretization error Tui are respectively equal to, −313h5120u(5)i,i=2 , and −146h5120u(5)i,3≤i≤N
5υi−8i+1+3υi+2=−2hυ'i−1+h33(7fi+4fi+1)+Tυi, i=1 (3.5)
−υi−1−3υi−3υi+1=h32(fi+fi+1)+Tυi, 2≤i≤N−1
−υi−1+4υi−3υi+1=−2hυ'i+1h36(3fi+fi+1)+Tυi, i=N
Where local discretization error Tυi are respectively equal to −13h5120u(5)i,i=1 , −2h515u(5)i 2≤i≤N−1 and −1h560u(5)i,i=N .
From (3.3)-(3.5), we have obtained local discretization error of the order O(h5) .Thus by neglecting the O(h5) term in (3.3)-(3.5), we will get our proposed difference method for the numerical solution of the problem (1.1). If we need to find v (b) corresponding to u"(b) the given boundary condition. We approximate v(b) for the boundary condition using u"(a),u"(b),u4(b) i.e.
υ(b)=u"(b)−u"(a)b−a+b−a6(2υ(4)(b)+u4(a)) (3.6)
But in this case the order of truncation error will not be O(h5)
Convergence analysis
We will consider following test equation for convergence analysis of the proposed method (2.4-2.5).
u6(x)=f(x,u), a<x<b. (4.1)
u(a) =α1; u"(a) =α2; u(4)(a) =α3; u(b) =β1; u(3)(b) =β2 and u(4)(b) =β3 Let w be the approximate solution of difference method (2.4-2.5) for numerical solution of the problem (4.1), we can write in the matrix form
Jw=Rh (4.2)
Where J is coefficient matrix, w=[u,v]T and Rh=[rh1,rh2]T these matrixes are
Rh=(α1−h2α2−h2α20⋮β1 + h32β2−2hα3+h33(7f1+4f2)h33(f2+f3)⋮h32(fN−1+fN)3β2−2hβ3+h36(3fN+fN+1))2N×1,w=(u1...uNυ1...υN)2N×1,
and let us define the coefficients matrix J in terms of block matrix,
J=(C1(u)…C2(v)⋮…⋮ C1(u)…C2(v) )2N×2N
Where
C1(u)=( 2 −1 0 0−1 2 −1 01 −3 3 −1 ⋱ ⋱ ⋱ ⋱ 1 −3 3 −1 1 −3 3 )N×N,C1(υ)=h312(17 −5 0 47 −23 18 −6 ⋱ ⋱ 18 −6 0 18)N×N
C2(υ)(5 −8 3 0−1 3 −3 1 ⋱ ⋱ ⋱ ⋱ −1 3 −3 1 −1 3 −3 0 −1 4 )N×N,
and matrix C2(υ)N×N depends on forcing function f(x,u) may be well defined. The exact solution W=[U,V]T of the difference method (2.4-2.5) satisfies the following equation
JW = Rh + T (4.3)
Where T = [Tu,Tv]T . Let T= (tm,1)2N1,
where
tm,1={43h5120u(5)m, m=1131h5120u(5)m, m=2 73h560u(5)m, 2<m≤ N 13h512υ(5)m, m=N+12h515υ(5)m, N+2≤m≤2N−1h560υ(5)m,m = 2N.
Let us define an error function the difference between approximate and exact solution of the difference method (2.4-2.5) i.e. E = w − W. Subtract (4.3) from (4.2) and substitute the above defined error into it, we will obtain
JE =T (4.4)
Let investigate the inevitability of the matrices C1 (u) and C2 (υ) . These matrices have different structure so we have to rely on computation of explicit inverse. Let explicit inverses of C1 (u) and C2 (υ) be respectively C−11 (υ)=(ki,j) N×N , where
ki,j={(N−i+1)(N+N), j = 1, i ≤ j ≤Ni(N−j+1)(N−j+2)2(N+N), j≠1,i≤j≤N,(N−i)(N−i+1)(N2−j(j−1)−1)−((N−i)2−1)(N−j+1)(N+j),2(N+1) j<i
and C2 (υ) =(li,j)N×N
li,j={(i−1)23a2,j−(i−4)23a1,j, i < j (N−j+1)24(N+N), j=1,i≤j≤N,(N−i+1)2(2j−N+8)4(N+1) 1 < j, j ≤ i ≤N
a2,j={N2+2N−6)−4(N−2j+2)((N+1)(N−2j)−(N−3))8(N+N), 2≤j<N2(N2+2N−6)8(N+N), jN2(N2+2N−6)−2(5N−16)(N−2j)2+8,8(N+1) N2<j≤N
a1,j={N2−(N−1)(N−2j)(N−2j+2)8(N+N), 2≤j<N2N(N+2)8, j=N2(N+1)(N+2Nj−2j2),4(N+1) N2<j≤N
It is easy to prove that matrices C1 (u) and C2 (υ) are positive. Let us define following terms,15
υup=∥C2(u)C−12(υ)∥, υlow=∥C1(υ)C−11(υ)∥
M*= (1 + υup) and M* = (1 + υlow).
Let us assume
M*M* <M*+M* and M=maxp=1,2∥C−1p(u/υ)∥
then matrix J is invertible15 and moreover
∥j−1∥≤MM*M*M*+M*−M*M* (4.5)
It is easy to prove that MM*M*M*+M*−M*M* is finite. Thus from (4.4) and (4.5), we have
∥E∥=∥TJ−1∥≤∥T∥MM*M*M*+M*−M*M* (4.6)
and ∥E∥ is bounded. It is easy to prove ∥E∥ tends to zero as h→0 . So we can conclude that finite difference method (2.4-2.5) converge and the order of convergence of the difference method (2.4-2.5) is at least O(h2)
Numerical Results
To test the computational efficiency of method (2.4), we have considered four model problems. In each model problem, we took uniform step size h. In Table 1 to Table 4, we have shown MAEU and MAEV the maximum absolute error in the solution u(x) and derivatives of solution v(x) of the problems (1.1) for different values of N. We have used the following formulas in computation of MAEU and MAEV:
MAEU = max1≤i≤N|u(xi)−ui|
MAEV = max1≤i≤N|u'(xi)−vi|
N |
||||
---|---|---|---|---|
16 |
32 |
64 |
128 |
|
MAEU.1 |
4035692(-2) |
.26897894(-3) |
.57788714(-6) |
.98105907(-7) |
MAEV. |
11683015(-4) |
.13530459(-6) |
.20602256(-6) |
.20949233(-5) |
Table 1 Maximum absolute error (Problem 1)
We have used Gauss Seidel and Newton-Raphson iteration method to solve respectively linear and nonlinear system of equations arised from equation (2.4). 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. The solutions are computed on N nodes and iteration is continued until either the maximum difference between two successive iterates is less than 10−6 or the number of iteration reached 103 .
Problem 1 The model linear problem given by
u(6)(x) = −xu(x) + (24 −11x + 2x2 + x3) exp(−x), 0 < x < 1
Subject to boundary conditions
u(0) = 0, u"(0) = 0, u(4)(0) = 8, u(1) = 2 exp(−1),
u(3)(1) = exp(−1) and u(4)(1) = 2 exp(−1).
The analytical solution of the problem is u(x) = x(1 + x) exp(−x). The MAEU and MAEV computed by method (2.4-2.5) for different values of N are presented in Table 1.
Problem 2 The model linear problem given by
u(6)(x) = −u(x) + 6(2xcos(x)+5sin(x)), 0<x<1
Subject to boundary conditions
u(0) = 0, u"(0) = 0, u(4)(0) = 0, u(1) = 0,
u(3)(1) = −6sin(1)+6cos(1) and u(4)(1)= −12sin(1)−8cos(1).
The analytical solution of the problem is u(x) = (x2−1) sin(x). The MAEU and MAEV computed by method (2.4-2.5) for different values of N are presented in Table 2.
N |
||||
---|---|---|---|---|
16 |
32 |
64 |
128 |
|
MAEU.1 |
.26380718(-2) |
.56681037(-3) |
.15497208(-5) |
.89406967(-7) |
MAEV. |
.71525574(-5) |
.95367432(-6) |
.23841858(-5) |
.23841858(-5) |
Table 2 Maximum absolute error (Problem 2)
Problem 3 The model nonlinear problem16 given by
u(6)(x) = u2(x) exp(x), 0 < x < 1
Subject to boundary conditions
u(0) = 1, u"(0) = 1, u(4)(0) = 1, u(1) = exp(−1),
u"(1) = exp(−1) and u(4)(1) = exp(−1).
The analytical solution of the problem is u(x) = exp(−x). The MAEU and MAEV computed by method (2.4-2.5) for different values of N are presented in Table 3.
N |
||||
---|---|---|---|---|
32 |
64 |
128 |
256 |
|
MAEU.1 |
.14109015(-2) |
.64671040(-4) |
.17881393(-6) |
.11920929(-6) |
MAEV. |
.25051057(-1) |
.25049835(-1) |
.25040478(-1) |
.25039405(-1) |
Table 3 Maximum absolute error (Problem 3)
Problem 4 The model nonlinear problem given by
u(6)(x) = 8sin(u2(x)exp(x)+cos2(x)exp(−1)), 0<x<1
Subject to boundary conditions
u(0) = 0, u"(0) = −2, u(4)(0) = 0, u(1) = sin(1) exp(−1),
u(3)(1) = 2sin(1)exp(−1)+2cos(1) exp(−1) and u(4)(1) = −4sin(1)exp(−1).
The analytical solution of the problem is u(x) = exp(−x). The MAEU and MAEV computed by method (2.4-2.5) for different values of N are presented in Table 4. The accuracy in numerical solution in considered model problems is satisfactory and increases as step size h decreases. The order of proposed method can be observed from the numerical experiment. Overall method is efficient and order of accuracy is at least quadratic. However inaccurate approximation in boundary condition affect the accuracy, it approved in numerical results obtained in considered example 3. If we do higher order approximation in boundary condition then this situation will destroy the matrix structure and proposed method may not converge. The advantage of the proposed method (2.4-2.5) is we get numerical approximation of third derivative of solution of problem as a byproduct which is otherwise useful.
N |
||||
---|---|---|---|---|
8 |
16 |
32 |
64 |
|
MAEU.1 |
.11576825(-2) |
.30554202(-3) |
.29175042(-4) |
.37556333(-7) |
MAEV. |
.46980762(-3) |
.49946433(-4) |
.80634038(-6) |
.64549158(-6) |
Table 4 Maximum absolute error (Problem 4)
To find the numerical solution of sixth order boundary value problems using finite difference method has been developed. At nodal point x = xi, i = 1, 2..,N we have obtained a system of algebraic equations given by (2.4-2.5) which is system of linear equations if source function f(x, u) is linear otherwise system of nonlinear equations. The propose method in numerical experiments has approved its efficiency and accuracy; moreover we get numerical approximation of third derivative as a byproduct. In future work, we will deal with improvement in present idea. Work in this direction is in progress.
None.
The author declares there is no conflict of interest.
©2018 Pandey. This is an open access article distributed under the terms of the, which permits unrestricted use, distribution, and build upon your work non-commercially.