Submit manuscript...
Open Access Journal of
eISSN: 2575-9086

Science

Research Article Volume 2 Issue 6

Solving numerically a sixth order differential equation as coupled finite difference equations approach

Pramod Kumar Pandey

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

Download PDF

Abstract

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

Introduction

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+1b 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),

ui1 + 2uiui+1 =h2u"i1+h312(17ui+5ui+1),   i=1  (2.4)

ui1 + 2uiui+1 =h2u"i2+h312(47ui+23ui+1),   i=2

ui2  3ui1+3uiui1 =h33(7fi+ui+1),  3iN

5ui 8ui+13ui+2 =2h2u'i1+h33(7f+4fi+1),   i=1  (2.5)

υi1_3υi3υ+1+υi+2=h32(fi+4fi+1),    2iN1

υi1_3υi3υ+1+υi+2=h32(fi+4fi+1),    2iN1

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+a1ui1+aoui+h2b1u"i1+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

(vi1+2vivi+1)+h2u"i1+h312(17fi5fi+1)+Tui =0 (3.3)

Where Tui  is local discretization error and equal to 43h3120u(5)i . Similarly we can derive the following equations

(vi1+2vivi+1)+h2u"i1+h312(17fi5fi+1)+Tui =0  (3.4)

ui23ui1+3uiui+1=h32(3ui+ui+1)+Tui,     3iN

Where local discretization error  Tui are respectively equal to, 313h5120u(5)i,i=2 , and 146h5120u(5)i,3iN

5υi8i+1+3υi+2=2hυ'i1+h33(7fi+4fi+1)+Tυi,  i=1  (3.5)

υi13υi3υi+1=h32(fi+fi+1)+Tυi, 2iN1 

υi1+4υi3υ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 2iN1 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)ba+ba6(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=(α1h2α2h2α20β1 + h32β22hα3+h33(7f1+4f2)h33(f2+f3)h32(fN1+fN)3β22hβ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                          01  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           01     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+2m2N1h560υ(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 C11 (υ)=(ki,j) N×N , where

ki,j={(Ni+1)(N+N),                                                                                        j = 1, i  j Ni(Nj+1)(Nj+2)2(N+N),                                                                      j1,ijN,(Ni)(Ni+1)(N2j(j1)1)((Ni)21)(Nj+1)(N+j),2(N+1)       j<i   

and C2 (υ) =(li,j)N×N

li,j={(i1)23a2,j(i4)23a1,j,       i < j (Nj+1)24(N+N),                     j=1,ijN,(Ni+1)2(2jN+8)4(N+1)      1 < j, j  i N   

a2,j={N2+2N6)4(N2j+2)((N+1)(N2j)(N3))8(N+N),      2j<N2(N2+2N6)8(N+N),                                                               jN2(N2+2N6)2(5N16)(N2j)2+8,8(N+1)                          N2<jN   

a1,j={N2(N1)(N2j)(N2j+2)8(N+N),      2j<N2N(N+2)8,                                          j=N2(N+1)(N+2Nj2j2),4(N+1)                      N2<jN   

It is easy to prove that matrices C1 (u)   and C2 (υ)  are positive. Let us define following terms,15

υup=C2(u)C12(υ),           υlow=C1(υ)C11(υ)

M*= (1 + υup) and M* = (1 + υlow).

Let us assume

M*M*  <M*+M*    and    M=maxp=1,2C1p(u/υ)

then matrix J is invertible15 and moreover

j1MM*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=TJ1TMM*M*M*+M*M*M*  (4.6)

and E is bounded. It is easy to prove E tends to zero as h0 . 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 = max1iN|u(xi)ui|

MAEV = max1iN|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 106 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) = (x21) 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)

Conclusion

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.

Acknowledgements

None.

Conflict of interest

The author declares there is no conflict of interest.

References

  1. Chandrasekhar S. Hydrodynamics and Hydromagnetic Stability. The Clarendon Press, Ox-ford, 1961.
  2. Agarwal RP. Boundary Value Problems for Higher Order Differential Equations. World Scientific, Singapore, 1986.
  3. Chawla MM, Katti CP. Finite Difference Methods for Two Point Boundary Value Problems Involving Higher Order Differential Equations. BIT. 1979;19:27–33.
  4. Pandey PK. Fourth order finite difference method for sixth order boundary value problems. Computational Mathematics and Mathematical Physics. 2013;53(1):57–62.
  5. Deacon AG, Osher S. Finite Element Method for Boundary Value Problem of Mixed Type. SIAM J Numer Anal. 1979;16:756–778.
  6. El-Gamel M, Cannon JR, Latour J. Sinc-Galerkin method for solving linear sixth order boundary-value problems. Mathematics of Computation. 2003;73(247):1325–1343.
  7. Siddiqi SS, Twizell EH. Spline Solutions of Linear Sixth Order Boundary Value Problems. Int J Comput Math. 1996;60(3-4):295–304.
  8. Loghmani GB, Ahmadinia M. Numerical Solution of Sixth Order Boundary Value Problems with Sixth Degree B-Spline Functions. Appl Math Comput. 2007;186:992–999.
  9. Che CHH, Kilicman A. On the solutions of nonlinear higher-order boundary value problems by using deferential transformation method and Adomian decomposition method. Mathematical Problems in Engineering. 2011;ID 724927:19.
  10. Ghazala A, Hamood UR. Solutions of a Class of Sixth Order Boundary Value Problems Using the Reproducing Kernel Space. Abstract and Applied Analysis. 2013;ID 560590:8.
  11. Wazwaz AM. The Numerical Solution of Sixth Order Boundary Value Problems by the Modified Decomposition Method. Appl Math Comput. 2001;118:311–325.
  12. Noor MA, Mohyud-din ST. Homotopy perturbation method for solving sixth-order boundary value problems. Computers & Mathematics with Applications. 2008;55(12):2953–2972.
  13. Pandey PK. The Numerical Solution of Third Order Differential Equation Containing the First Derivative. Neural Parallel & Scientific Comp. 2005;13:297–304.
  14. Pandey PK. Solving third-order Boundary Value Problems with Quartic Splines. Springer- plus. 2016;5(1):1–10.
  15. Gil MI. Inevitability Conditions for Block Matrices and Estimates for Norms of Inverse Matrices. Rocky Mountain Journal of Mathematics. 2003;33(4):1323–1335.
  16. Khalid M, Sultana M, Zaidi F. Numerical Solution of Sixth-Order Differential Equations Arising in Astrophysics by Neural Network. International Journal of Computer Applications. 2014;107(6):1–6.
Creative Commons Attribution License

©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.