 American Journal of Computational Mathematics, 2011, 1, 318-327 doi:10.4236/ajcm.2011.14038 Published Online December 2011 (http://www.SciRP.org/journal/ajcm) Copyright © 2011 SciRes. AJCM A New Fourth Order Difference Approximation for the Solution of Three-Dimensional Non-Linear Biharmonic Equations Using Coupled Approach Ranjan Kumar Mohanty1, Mahinder Kumar Jain2*, Biranchi Narayan Mishra3# 1Department of Mat hem at i cs, Faculty of Mathematical Sciences, University of Delhi, Delhi, India 2Department of Mat hem at i cs, Indian Institute of Technology, New Delhi, Indi a 3Department of Mat hem at i cs, Utkal University, Bhubaneswar, India E-mail: rmohanty@maths.du.ac.in, bn_misramath@hotmail.com Received October 13, 2011; revised November 2, 2011; accepted November 9, 2011 Abstract This paper deals with a new higher order compact difference scheme, which is, O(h4) using coupled ap- proach on the 19-point 3D stencil for the solution of three dimensional nonlinear biharmonic equations. At each internal grid point, the solution u(x, y, z) and its Laplacian 2u are obtained. The resulting stencil al- gorithm is presented and hence this new algorithm can be easily incorporated to solve many problems. The present discretization allows us to use the Dirichlet boundary conditions only and there is no need to discre- tize the derivative boundary conditions near the boundary. We also show that special treatment is required to handle the boundary conditions. Convergence analysis for a model problem is briefly discussed. The method is tested on three problems and compares very favourably with the corresponding second order approxima- tion which we also discuss using coupled approach. Keywords: Three-Dimensional Non-Linear Biharmonic Equation, Finite Differences, Fourth Order Accuracy, Compact Discretization, Block-Block-Tridiagonal, Tangential Derivatives, Laplacian, Stream Function, Reynolds Number 1. Introduction We are interested to develop a new algorithm for solving the three-dimensional non-linear biharmonic partial dif- ferential equation 44444444442222 2222 22,,2,,,, ,,,,,,0,,1xyzx y zuuuuxyz xyzuuuxy yz zxfxyzuu uuuuuuxyx  (1) defined in the solution region ,, 0,,1xyz xyz with boundary , where 222222,, uuuxyzmensional Laplacians of u in . We assume that the solution u(x, y, z) is smooth enough to maintain the order and accuracy of the algorithm as high as possible under consideration. The Dirichlet boundary conditions are given by  2122,, ,,, , ,,uufxyzfxyz xyzn. (2) The nonlinear biharmonic equation is a fourth order elliptic partial differential equation which occurs in many areas of physics and applied mathematics, especially in elasticity theory and stokes flow problems. It occurs with fluid flowing over an obstacle or with the movement of a natural or artificial body. Common examples are the flows past an airplane, a submarine, an automobile, or wind blowing past a high-rise building. Till about three and half decades ago, second order accuracy was consid- ered to be sufficient for most biharmonic problems. In particular, the central difference schemes have been the most popular ones because of their straightforwardness 2uxyz represents the three di- *Present address: 4076, C/4, Vasant Kunj, New Delhi, India. #Present address: Department of Mathematics, Rajasunakhala College, Nayagarh, Orissa, India. R. K. MOHANTY ET AL.319 in application. Though for problems having well behaved solutions, the solution may be of poor quality for con- vection dominated flows, or for high Reynolds number if the mesh is not sufficiently refined. Again, higher order discretization is generally associated with non-compact stencil which increase the band-width of the resultant coefficient matrix. Both mesh refinement and increased matrix band-width invariably lead to a large number of arithmetic operations. Thus neither a lower order accu- rate method on a fine mesh nor a higher order accurate one on a non-compact stencil seems to be computation- ally cost-effective. This is where higher order compact (HOC) finite difference methods become important. A compact finite difference schemes is one which utilizes grid points located only directly adjacent to the node about which differences are taken. In addition, if the scheme has an accuracy greater than two, it is termed as HOC method. The higher order accuracy of the HOC methods combined with the compactness of the differ- ence stencils yields highly accurate numerical solutions on relatively coarser grids with greater computational efficiency. A compact difference scheme is one that is restricted to the patch of cells immediately surrounding any given node and does not extend further. Most stan- dard difference schemes such as the central difference scheme for second order elliptic partial differential equa- tions are compact. The high order compact method con- sidered here is different in that the governing differential equation is used to approximate the lower order deriva- tive terms with the imbedding technique. The scheme is difficult to develop due to the need for extensive alge- braic manipulation, especially for non-linear problems. However, once high-order method developed, it can be incorporated easily in application. Various approaches for the numerical solution of 2D biharmonic problems has been discussed in the literature. Not many authors have tried to solve the three-dimensional biharmonic problems. The reason is that for high order approxima- tion it is difficult to discretize the nonlinear biharmonic equation as well as the associated boundary conditions using a single computational cell. Another reason of ne-cessity is that three-dimensional problems require large computing power and place huge amount of memory requirements on the computational systems. Such com- puting power has only recently begun to become avail- able for academic research. Smith  and Ehrlich [2,3] have used coupled approach and solved 2D biharmonic equation using second order accurate finite difference equations. Bauer and Riess  have used block iterative method to solve block 5-diagonal matrices for the solu- tion of 2D biharmonic problem of first kind. Glowinski and Pironneau  have developed a stable lower order numerical method for the first biharmonic problem. Later, kwon et al. , Stephenson , and Mohanty et al. [8-12] have developed certain second- and fourth-order finite difference methods for the solution of two and three di-mensional biharmonic problems using single compact cells. Recently, using more grid points Singh et al  and Khattar et al.  have discussed the fourth order numerical methods for the solution of 2D and 3D bihar- monic problems of second kind. Further, using coupled approach Mohanty  has derived a new nine point fourth order finite difference method for the solution of non-linear biharmonic problems of second kind. In most recently, using nine point compact cell Mohanty et al.  have discussed fourth order finite difference method for the solution of nonlinear 2D triharmonic problems. Ear- lier, Mohanty et al. [17-19] have developed fourth order compact difference schemes for the solution of multi- dimensional non-linear elliptic partial differential equa-tions. Fourth order compact difference schemes have become quite popular as against the other lower order accurate schemes which require high mesh refinement and hence are computationally inefficient. On the other hand, the higher order accuracy of the fourth order com-pact methods combined with the compactness of the dif-ference stencil yields highly accurate numerical solutions on relatively coarser grids with greater computational efficiency. A conventional approach for solving the 3D biharmonic equation is to discretize the differential Equa- tion (1) on a uniform grid using 125-point approxima- tions with truncation error of order . This approxima- tion connect the values of central point in terms of 124 neighbouring values of u in 55 cubic grid. We note that the central value of u is connected to grid points two grids away in each direction from the central point and the difference approximations needs to be modified at grid points near the boundaries. There are serious computational difficulties with solution of the linear and non-linear systems obtained through 125-point discreti-zation of the 3D biharmonic equation. Approximations using compact cells avoid these difficulties. The compact approach involves discretizing the biharmonic equations using not just the grid values of the unknown solution u but also the values of the derivatives 2h5xx, uyy and uzz at selected grid points (see Mohanty and Pandey ). It is required to solve the system of four equations to obtain the solutions of u,xxuu, yy and uzzu. In this paper, we write the original differential equation in a coupled man-ner and introduce new concepts to handle the boundary conditions without discretizing them. For fourth order approximations we use only 19-point compact cell (see Figure 1). No approximations for derivatives need to be carried out at the boundaries and the given Dirichlet boundary conditions are exactly satisfied. The main ad- vantage of this work is that we require to solve a system of two equations, whereas in our previous work , we were required to solve a system of four equations to ob- Copyright © 2011 SciRes. AJCM R. K. MOHANTY ET AL. Copyright © 2011 SciRes. AJCM 320 (xl, ym+1,zn+1)(xl+1, ym, zn+1) (xl-1, ym, zn+1) (xl-1, ym, zn-1) (xl+1, ym, zn-1) (xl, ym-1, zn-1) (xl+1, ym+1, zn) (xl-1, ym+1, zn) (xl, ym-1, zn+1) (xl-1, ym, zn) (xl-1, ym-1, zn) (xl+1, ym-1, zn) (xl, ym, zn-1) (xl, ym+1, zn-1) (xl, ym, zn) (xl, ym, zn+1) (xl, ym-1, zn) (xl+1, ym, zn) (xl, ym+1, zn) Figure 1. 19-point 3D single computational cell. tain the numerical solution of u(x, y, z). Thus the pro- posed method requires less algebraic operations as com- pared to the method discussed in . In Section 2, we give the formulation of the method. The complete deri- vation of the method is carried out in Section 3. In order to validate the proposed fourth order method and test its robustness, we solve three problems in Section 4 and also discussed stability analysis briefly for a model pro- blem. Concluding remarks are presented in Section 5. 2. Formulation of the Fourth Order Discretization Consider the three-dimensional solution region , which is replaced by a set of grid points (xl, ym, zn), where h > 0 is the mesh sizes in x-, y- and z- directions, and grid points are defined by xl = lh, ym = mh, zn = nh; l, m, n = 0 (1) N + 1 with (N + 1) h = 1. Let ,,lmn and ,,lmn be the approximate and exact solution values of u(x, y, z) at the grid point (xl, ym, zn ), respectively. uUThe Dirichlet boundary conditions are given by (2). Since the grid lines are parallel to coordinate axes and the values of u are exactly known on the boundary, this implies, the successive tangential partial derivatives of u are known exactly on the boundary. For example, on the plane y = 0, the values of u(x, 0, z) and are known, i.e., the values of xux , , , ,··· etc are known on the plane (,0,)yyuxz,)z(,0zux(,0 ,)z(,0,)xxux z(,0,)zzux zy = 0. This implies, the values of u(x, 0, z) and 2, 0,, 0,, 0,,0,xx yy zzux zux z uxz ux z  2u are known on the plane y=0. Similarly the values of u and are known on all plane sides of the cubic region . Let us denote . Then we re- formulate the boundary value problems (1) and (2) in a coupled manner as 2(, ,)(, ,)uxyz vxyz 2222222,,,, ,,,uuuuxyz vxyzxyzxyz, (3a) 2222222,,,,,,,, ,,, ,,,,xxyyzzvvvvxyz xyzfxyzuvu v uvu vxyz. (3b) subject to the Dirichlet boundary conditions prescribed by  ,, ,,, ,,,uaxyzv cxyzxyz. (4) Let at the grid points (xl, ym, zn), the approximate and exact solution values of v(x, y, z) be denoted as and , respectively. ,,lmnv,,lmnIn order to obtain fourth order approximations on the 19-point compact cell for the system of non-linear dif- ferential Equations (3a) and (3b), we need the following approximations: V,, 1, ,1, ,12xlm nlmn lmnUUUh, (5a) ,, 1, ,1, ,12xl m nlmn lmnVVVh, (5b) ,, ,1, ,1,12yl m nlm nlm nUUUh, (6a) ,, ,1, ,1,12ylm nlm nlm nVVVh, (6b) ,, ,, 1,, 112zl m nlmn lmnUUUh, (7a) R. K. MOHANTY ET AL.321 ,, ,, 1,, 112zl m nlmn lmnVVVh, (7b) 1, ,1, ,, ,1, ,1342xlm nlmnlmn lmnUUUUh , (8a) 1, ,1, ,,,1, ,1342xlm nlmnlmn lmnVVVVh , (8b) ,1, 1, 1,1, 1,12xl mnlmnlmnUUUh , (9a) ,1,1, 1,1, 1,12xl mnlmnlmnVVVh, (9b) ,, 11,,11, ,112xlm nlmn lmnUUUh , (10a) ,, 11,,11, ,112xlm nlmn lmnVVVh , (10b) 1, ,1,1,1, 1,12ylm nlmnlmnUUUh , (11a) 1, ,1, 1,1, 1,12ylm nlmnlmnVVVh , (11b) ,1, ,1,,, ,1,1342yl mnlm nlmnlmnUUUUh , (12a) ,1, ,1,,, ,1,1342yl mnlm nlmnlm nVVVVh , (12b) ,, 1,1,1 ,1,112yl m nlmnlm nUUUh , (13a) ,, 1,1,1 ,1,112ylm nlm nlm nVVVh , (13b) 1, ,1,,11, ,112zlm nlmnlmnUUUh, (14a) 1, ,1, ,11, ,112zlm nlmnlmnVVVh , (14b) ,1, ,1,1,1,112zl mnlmnlm nUUUh , (15a) ,1,,1,1 ,1,112zl mnlmnlm nVVVh , (15b) ,, 1,, 1,,,, 11342zlmnlmnlmn lmnUUUUh , (16a) ,, 1,, 1,,,, 11342zl m nlmnlmn lmnVVVVh , (16b) ,1, 1,1,,1,1,1,212xxl mnlmnlmn lmnUUUUh , (17a) ,1, 1, 1,, 1,1, 1,212xxl mnlmnlmnlmnVVVVh , (17b) ,, 11,, 1,, 11,, 1212xxl m nlmnlmn lmnUUUUh , (18a) ,, 11,, 1,, 11,, 1212xxl m nlmnlmn lmnVVVVh , (18b) 1, ,1,1,1, ,1,1,212yylm nlmnlmn lmnUUUUh , (19a) 1, ,1,1,1,,1,1,212yylm nlmnlmn lmnVVVVh , (19b) ,, 1,1,1,,1 ,1,1212yylm nlm nlmnlm nUUUUh , (20a) ,, 1,1,1,,1 ,1,1212yylm nlm nlmnlm nVVVVh , (20b) 1, ,1,,11, ,1,,1212zzlm nlmnlmn lmnUUUUh , (21a) 1, ,1, ,11,,1, ,1212zzlm nlmnlmnlmnVVVVh , (21b) ,1,,1,1,1, ,1,1212zzl mnlm nlm nlm nUUUUh  , (22a) ,1, ,1,1,1, ,1,1212zzl mnlm nlm nlmnVVVVh  . (22b) Then we need to evaluate 1, ,1, ,1, ,1, ,1, ,1, ,1, ,11,,1,,,,,, ,,,,,,lmnxlmn xlmnylmnylmnzlmn zlmnlmnlmnlmnFfxyzUVUVUVUV  (23) ,1,,1, ,1,,1,,1,,1, ,1,1,1,,1,,,,, ,,,,,,lm nxlm nxlmnylm nylm nzlm nzlmnlmnlmnlmnFfxyzUVUVUVUV  (24) ,, 1,,1,,1,,1,, 1,, 1,, 11,,1,,1,, ,,,,,,,,lmn xlmn xlmnylmnylmnzlmn zlmnlmnlmn lmnFfxyzUVUVUVUV  (25) Further, we define ,, ,,1,,1,,1,,1,,1, ,1, ,12 1212xlmnxlmnyyl mnyyl mnzzlmnzzl mnlmn lmnhh hUUVVUUU U   (26a) 1, ,1, ,,, ,,1,,1,,1,,1,,12 1212lmn lmnxlmnxlmnyyl mnyyl mnzzlmnzzl mnhh hVVFFVVVV   (26b) Copyright © 2011 SciRes. AJCM R. K. MOHANTY ET AL. 322 ,,,,, 1,, 1,, 1,, 1,,1,,1,12 1212ylmnylmnxxlm nxxlm nzzlm nzzlm nlm nlmnhh hUUVVU UU U   (27a) ,1,,1,,,,,, 1,, 1,,1,, 1,12 1212lm nlm nylmnylmnxxlm nxxlmnzzlmnzzlm nhh hVVFFVVVV    (27b) ,,,,,, 1,, 1,, 1,, 1,, 1,,112 1212zlmn zlmnxxlmnxxlmnyylmnyylmnlmn lmnhh hUUVVU UUU   (28a) ,, 1,, 1,,,,,, 1,, 1,, 1,, 112 1212lmn lmnzlmn zlmnxxlmnxxlmnyylmnyylmnhh hVVFFVVVV    (28b) Finally, let ,, ,,,,,,,,,, ,,,, ,,,,,, ,,,,,,lmnxlmn xlmnylmn ylmnzlmn zlmnlmn lmnlmnFfxyzUVUVUVUV (29) Then at each internal grid point (xl, ym, zn) of the solu-tion region , the given system of differential Equa-tions (3) are discretized by ,1, 11,, 1,, 11,, 1,1, 11,1,,1,1,1,1,,, ,1,,1,1,,1,1,1,,1,11, ,1, ,11, ,1,1,2224 222lm nlmnlmnlmnlm nlm nlm nl m nl mnlmnl mnl m nlm nlm nlm nlmnlmnlmnlmLUU UUUU UUUUUUUUU U UUU U      21261,,1,, ,1, ,1,,,1,,1,,6()2,, 11nlmnlmn lmn lmnlmnlmnlmnhVVVVVVVOhlmn N (30a) ,1, 11,, 1,, 11,, 1,1, 11,1,,1,1,1,1,,,,1,,1, 1,, 1,1, 1,, 1,11,,1,,11,,1, 1,2222 2lmnlmnlmnlmnlmnlmnlmnlmnlmnlml mnlmnlmnlmnlmnlmnlmnlmnlmLVV VVV V VVVVVVVVV V VVV V     224n121,,1,, ,1, ,1, ,,1,,1,,,,6,2 ,,11nlmn lmn lmn lmn lmnlmnlmnlmnhFFFFFFFTlmn N  (30b) where 6,,lmnTOh. Note that, the approximations (30a) and (30b) require only 19-grid points with a single computational cell. Incorporating the Dirichlet boundary conditions given by (4) into the difference methods (30a) and (30b), we obtain the sparse system of tri-block-block diagonal matrix equations in coupled form, which can be solved by appropriate iterative methods (see [20-24]). 3. Derivation of the Fourth Order Approximations At the grid point ,,lmnxyz, let us denote (1)(2) (1)(2),,,, ,,,,,, ,,,, ,,(1) (2),, ,,,, ,,,,,,,ijkijklmnlmn lmnlmnijk xlm nxl m nyl m nylm nlm nlmn lmnzl mnzlmnUfffUUVUVxy zffUV ,f (31) Further, at the grid point ,,lmnxyz, we define ,,,, ,,,,,,,,,,,,,,,,,, ,,,,,,lmnlm nlmn lmnxlmnxlmnylmnylmnzlmnzlmnFfxyzUVUVUVUV (32) Using Taylor expansion about the grid point ,,lmnxyz, we first obtain 261,,1,,,1,,1, ,,1,,1,,62lmn lmn lmn lmnlmnlmnlmnhLVFFFFFFFOh   (33) Copyright © 2011 SciRes. AJCM R. K. MOHANTY ET AL. Copyright © 2011 SciRes. AJCM 323 Now by the help of the approximations (8a)-(16b), from (23)-(25), we obtain 231, ,1, ,16lmn lmn hFFTOh, (34) 23,1,,1, 26lm nlm nhFFTOh, (35) 23,, 1,, 136lmn lmnhFFTOh. (36) where (1)(2) (1)1300, ,300, ,030, ,(2)(1) (2)030, ,003, ,003,,22lmnlmn lmnlmnlmn lmnTUV UVUV  , (1)(2 )(1)2300 ,,300 ,,030 ,,(2)(1) (2)030, ,003, ,003, ,22lmn lmnlmnlmnlmnlmnTU VUVUV (1) (2)(1), 3300, ,300, ,030,,(2)(1) (2)030, ,003, ,003, ,22lmn lmnlmnlmnlmnlmnTUV UVUV . Let us consider the linear combination ,,,, 111, ,1, ,1, ,1, ,121, ,1,,13xlmn xlmnlmn lmnyylm nyylm nzzlm nzzlm nUUhaVVha UUha UU (37a) 1, ,1, ,,,,, 111, ,1, ,121, ,1,,13lmn lmnxlmn xlmnyylm nyylm nzzlm nzzlm nVVhbF Fhb VVhb VV , (37b) ,,,, 21, 1,,1,,1,,1,22,1,,1,23ylm nyl m nlm nlm nxxl mnxxl mnzzl mnzzl mnUUhaVVha UUha UU , (38a) ,1,,1,,,,, 21,1,,1,22,1,,1,23lm nlm nylm nylm nxxl mnxxl mnzzl mnzzl mnVVhbFFhb VVhb VV , (38b) ,,,, 31, ,1, ,1,, 1,, 132,, 1,, 133zl mnzlmnlmn lmnxxlm nxxlm nyylm nyylm nUUhaVVha UUha UU , (39a) ,, 1,, 1,,,, 31,, 1,, 132,, 1,,133lmn lmnzl m nzl m nxxlm nxxlm nyylm nyylm nVVhbF Fhb VVhb VV  (39b) where and i, j = 1, 2, 3 are parameters to be determ. Now using the approximations (5a)-(7b), (1and (34)-(36), from (37a)-(39b), we obtain ijainedijb; 7a)-(22b) 4,, ,, 46xlm nxl m nUU TOh, (40a) 2h24,, ,,56xlm nxlm nhVV TOh , (40b) 24,, ,, 66ylm nylm nhUU TOh, (41a) 24,, ,, 76ylm nyl m nhVV TOh , (41b) 24,, ,, 86zl m nzl m nhUU TOh , (42a) . (42b) 24,, ,, 96zl m nzl m nhVV TOhwhere 411 3001112 1201113102112 1212TaUaaaaU  , U11 3001112 1201113 102112 1212TbVbbVbb , V5621 03021222102123 012112 1212TaUaaaaU U, 72103021 2221021 23012112 1212TbVbbVbb V, 831 00331322013133 021112 1212TaUaaaaU U, 931003313220131 33021112 1212TbVbbVbb VBy the help of the approximations (40a)-(42b), from (29), we get . 24,, ,, 106lmn lmn hFFTO  h (43 wly, by the help of the relations (34)-(36) and (43), 0b) and (33), we obtain the local truncation error )here 10T(1)(2) (1)(2)4,, 5,, 6,,7,,(1) (2)8,,9,,lmn lmnlmnlmnlmn lmnTT T TTT . Finalfrom (346,, 123 1036lmn hTTTTTO h (44) The proposed difference methods (30a) and (30b) are to be of 4Oh , the coefficient of 4h in (44) must be zero and obtain a relation we R. K. MOHANTY ET AL. 324 Substituting the values of and in (45), e val123 1030TTT T  (45) 1T, 23,TT 10Twe obtain thues of parameters 11112iab for i = 1, 2, 3 and i1ab for 12ij iji = 1, 2, th3 and j = 2, 3, and e local truncation error (44) reduces to 6,lmTOh. We may summarize the results as fo,nllows: int finite dimethods (30a) and (30b) with the approximations listed of the non-li- t boundarco (46) the abTheorem 1: The compact 19-pofference in (5a)-(29) are of 4for the numerical solutions of Oh,,uxyzand its Laplacian near biharmonic Equation (1) with Diriy 2,,uxyzchlenditions (2). 4. Stability Consideration and Results of Computational Experiments Let us consider the test equation 4,,,, ,0,,1uxyz Gxyzxyz  Applying the proposed method (30a) and (30b) toove equation, we obtain , 1,11,,1,,2lm nlmnlmUU U 1 1,,1,1,11,1,,1,1,1,,1,1 1,,1,, 11,, 1,1, 1222422nl mnlm nlmnlm nl m nlm nlmnlmn lmnlmnUUU UUUUUUUUU      1, ,,,1, ,1,1,,1, 1,1,2lmnlmnlmn lmnlm nl m nUU UU  21, ,1, ,,1,,1,,, 1,,1,,,, 11lmnlmnlmnlmnlmn lmnlmnhVVVVVV VlmnN  (47a) ,1, 11,, 1,,11,, 1, 1,11, 1,, 1,1, 1,1, ,,,1, ,1,1,,1,1,1,,1,1 1,,1,, 11,, 1,1, 122224222lm nlmnlmnlmnlmnlmnlmn lmnlmnlmnlmn lmnlm nl m nlmnl mnlmn lmn lmnVV VVVV VVVVVVVV V VVV V        21, ,1, ,,1,,1,,, 1,, 1,,26,, 11lmnlmnlmnlmlmn lmnlmnhGGGGGG Glmn N  n (47b) where , ... etc. An ) and (47b) can be written as ,, ,,lmnlm nGGxyziterative method for (47a 2k+1 kk24= 2hIuAuBvRHU k+1 kk24= Iv0u+ AvRHV  kk (48b) where are solution vectors and RHU, RHV are righe vectors consists of boundary and ho- mogenous function values. Above iterative method in matrix form can be written as ,uv t hand sidk+1 kk+1 kUUGRVV H (49) where 21,RH .2h24 RHUABG= We denote I = [0, 1, 0] as the Nth order and H = [1, 2, 1] and Q = [2, 0, 2] as Nth onal matrices, and C = [I, H, I] and D = [HNth order block-tridiagonal matrices, whete RHV0Aidentity matrix rder tridiago-, Q, H] as the re, in general, we denobc,,NNabcabcab ab c0 0as Nth given byorder tridiagonal matrix whose eigen values are π2cos ,1,2,1jbacj NN . (50) The eigenvalues of I, H and Q are 1(N-times), 22cosπkhand 4cos πkh, 11kN respectively, where (N + 1)h = 1. The eigenvalues of C and D are given by  22cosπcos π;, 11jk jhkh jkN  (51a) and 4cos πcos π)cosπcos π,11jk jhkhjh khjk N(51b) The matrix A = [C, D, C] associated with the iteration matrix G is a Nth order block-block-tridiagonal matrix whose eigenvalues are given by  2cosπ;, ,11ijk jkjkih ijkN or,   4cos πcosπcosπcos πcos π;ih jhkhkh ihπcosπcosπcos πcosijkkhih jhjh (48a) ,, 11ijkN (52) Copyright © 2011 SciRes. AJCM R. K. MOHANTY ET AL. Copyright © 2011 SciRes. AJCM 325Similarly, the eigenvalues of the matrix B as iteration matrix G are given by where Gis the spectral radius of G. sociated The cteristic equation of the matrix G is given by with the62cosπcos πcos π;,, 11ijk ihjh khijkN  harac2124 481 (53) 24 ijk0Thus the eigenvalues of G are given bydet0;,,11hijkijk ijkN  (54) The iterative method (49) is stable as long as1G,   πcos πcos πcos πcos πcos π;ihjhjhkhkh ih (55) The maximum values of all eigenvalues of Gccur at . Hence, 11246 cos πcosπcos πcosijkijk ihjh kh   o1jk i  cos πmax .1cosπ12ijkhh G, (56) which is satisfied for all variable angles πh. Hence the itethe system of differential Equations (3a) and (3b) are straightforward and can be written in a coupled manner (57a) rative method (48a)-(48b) is stable. The second order approximations for 24,,1,1,1,,,,1,, ,1, ,,1,,6,,, 11lmnlmn lmnlmnlmn lmn lmnlmnUUUUUUU hVOhlmn N ,, ,,,, 11lmn lmnlmnlmn N,,1,1,1,,,,1,, ,1, ,,12 4,, ,,,,,,,,,,6,,,, ,,,,,,lmnlmn lmnlmn lmnlmn lmnxlmnxlmnylmnylmnzlmn zlmnVVV VVVVhf xyzUVUVUVUVOh (57b) Note that, the second order approximations (57a) and (57b) require only 7-grid points on a single cotional cell (see Figure 1). By combining the difference equations at each internal grid points, we obtain a large sparse system of matrix to solve. At each interior mesh point, we have two un- knowns u and , that is, the number of bands with non-zero enincreased, and so is the size of the final matrix for the same mesh size. However, by this new method, the value of the Laplacian, which is often of interest, is also computed. Whenever mputa- 2uvtries is ,,,,,,,,, ,xxyyzzfxyzuvu v uv u v is linear (or, non-linear) in ,,, ,xxuvu v ,,yyzuvu andzv, the difference Equations (30a) and (30b) form a linear (or, non-linear) system. To solve such a system or indeed to de double precision arithmetic. Test Problem 1: (Biharmonic problem) calculated maximum absolute errors (l-norm) for dif- ferent grid sizes. All computations were performed using 4444 44u444 2222222,, ,0,,1uuuu uxyz xyyzzxGxyzxyz    (58) The exact solution is  ,, 1cos2π1cos2π1cos2πuxyzx y z. The maximum absolute errors are tabulated in Table 1. Table 1. Test Problem 1: The maximum absolute errors. h 4-methodOh 2-methodOh 1/8 U 2u 0.1065(–01) 0.6145(+00) 0.1089(–01) 0.2932(+02) 1/16 U 2u 0.6659(–03) 0.3776(–01) 0.2614(+00) 0.7121(+0monstrate the existence of a solution, we use iterative methods. In this section, we solve the following three test problems in the region 0 < x, y, z < 1, whose exact solu- tions are known. The Dirichlet boundary conditions and right hand side homogeneous functions are obtained by using the exact solutions. We have also compared the numerical results obtained by proposed fourth order ap- proximations (30a) and (30b) with the numerical results ob o1) Utained by cximations (57a) and (57b). In allered (0) u0 as the initial approximation and the iterations were stopped when the absolute error tolerance  rresponding second order appro cases, we have consid 11210kkuu was achieved. In all cases, we have 1/32 2u  0.2349(v02) 0.1767(+01) 1/64 U 2u0.4157(–04) 0.6469(–01)  0.2586(–05) 0.1466(–03)0.1611(–01) 0.4410(+00) R. K. MOHANTY ET AL. 326 blems) 2Test Problem 2: (Variable coefficient pro422221111sin1sin,, ,xxx xyy xzzxxy yyy yzzxxz yyz zzz21sin0,,1xyzuxuuuyu u uzu u uxuyzu Gxyz   u (59a) The exact solution is xyz,, sinπsin πsin πuxyzx y z. 4244(1 cos)()(1 co(1 cos)()(1) (1(,,1xxx xyy xzzyyy yzzxxz yyz zzz2s)( xxy24))) (1,),0 ,xyzuxuuuu uzu u uyuxuy uGxy z  zuz xy  (59b) The exact solution is u(x, y, z) =exyz . re tabulaThe mTest Problem 3: (Navier-Stokes model equation in terms of stream function ) (see ) 41()((, ,),0, ,1)yxxxxyyx xxyyyyeRGxyz xyz   (60) The exact solution is ,, esinπsin πxxyzy z. The maximum absolute errors are tabulated in Table 3 for various values of Reynolds number . eR5. Concluding Remarks In this paper, using coupled approach we discuss a new fourth order 19-point compact finite difference approxi- mation for the solution of 3D non-linear biharmonic el- liptic partial differential equations. The method is de- rived on a single computational cell using the values of u and 2uas the unknowns. We have obtained the nu- merical solution of 2u as a by-product, which is quite often of interest in many physical problems. Our method is applied to solve several problems including Navier Stokes model equation in terms of stream function  and enables us to obtain high accuracy solutions with great efficiency. Numerical results confirm that the pro- axm absoae 2. imulute errors ted in Tabl em 2: The maximum absolute errors. Problem (59a) Problem (59b) Table 2. Test Problh 4-MethodOh 2-MethodOh 4-Met ho dOh 2-Met ho dOh 1/4 u 0.7778(–02) 0.1161(+00) 0.1078(+00) +01) 0.1464(–04) 0.6745(–03) 0.5088(–02) 0.3249(–01) 2u 0.1562(1/8 U 0.4655(–03) 0.6952(–02) 0.2578(–01) 0.3812(+00) 0.8861(–06) 0.4195(-04) 0.1444(–02) 0.1016(–01) 1/16 U 0.2877(–04) 0.4551(–03) 0.6377(–02) 0.9474(–01) 0.5529(–07) 0.2626(–05) 0.3699(–03) 0.2577(–02) 0.1797(–05) 0.2836(–04) 0.1590(–02) 0.2393(–01) 0.3456(–08) 0.1666(–06) 0.9353(–04) 0.6467(–03) 2u 2u 1/32 U 2u Table 3. Test Problem 3: The maximum absolute errors. 4-Met ho dOh 2-Met ho dOh h 210eR 46810 ,10 ,10eR 246810 ,10 ,10 ,10eR  2 0.1808(–03) 0.3332(–02) 0.1880(–03) 0.3524(–02) Over Flow 1/8 1/16  2 0.2079(–00.1134(–04) 3) 0.1202(–04) 0.2253(–03) Over Flow 1/32  2 0.7095(–06) 0.1299(–04) 0.7545(–06) 1413(–04) Over Flow 1/64 0. 2 0.4432(–07) 0.8127(–00.4705(–07) 0.8808(–06) w 6) Over Flo Copyright © 2011 SciRes. AJCM R. K. MOHANTY ET AL.327 posed fourth order method produces oscillation free so-lution for large Reynolds number, whereas the second or- der method is unstable. 6. References  J. Smith, “The pled Equation Approach to the Nu-merical Solution of the Bihaation by FDifferences,” SIAM Journal on Analysis, V7, No. 1, 1970, p04-111. doi:10.1137/0707005Cou rmonic EquNumerical inite ol. p. 1  L. W. Solving the Buation as Cpled Finite Difference Equations,” SIAM Journal on Nu- merical Analysis, Vol. 8, No. 2doi:10.1137/0708029Ehrlich, “iharmonic Eqou- , 1971, pp. 278-287.  L. W. EhrlichPoint and Block SOR Applied to a Cou- pled Set of Difference Equations,” Computing, Vol. 12, No. 3, 1974, pp. 181-194. doi:10.1007/BF02293104, “  L. Band E. L. Riess, “Bloconal Matrices and thst Numerical Solutionrmonic Equa- tion,” Mathematics of Computation, Vol. 26, No. 118, 1972, pp. 311-326. doi:10.1 90/S00250312751-9uer ae Fak Five Diag of the Biha0 -5718-1972-  R. Glowinski and O. Pironneau, the First Biharmonic EquationsTwo-Dimen- sional Stokes Problems,” SIAM l. 21, No. 2, 1979, pp. 167-212. doi:10.1137/1021028“Numerical Methods for and for the Review, Vo  Y. Kwon, R. Manohar and J. W., “Single Cell Fourth Order Methodn,” Con- armonic Problems,” Journal of Com-, No. 1, 1984, pp. 65-80. 21-9991(84)90015-9 Stephensonarmonic Equatios for the Bihgress Numerantium, Vol. 34, 1982, pp. 475-482.  J. W. Stephenson, “Single Cell Discretization of Order Two and Four for Bihputational Physics, Vol. 55 doi:10.1016/00 y and P. K. Pandey, “Difference Methods of  R. K. MohantOrder Two and Four for Systems of Mildly Non-linear Biharmonic Problems of Second Kind in Two Space Di- mensions,” Numerical Methods for Partial Differential Equations, Vol. 12, No. 6, 1996, pp. 707-717. doi:10.1002/(SICI)1098-2426(199611)12:6<707::AID-NUM4>3.0.CO;2-W  R. K. Mohanty, M. K. Jain and P. K. Pandey, “Finite Dif- ference Methods of Order Two and Four for 2D Nonear Biharmonic Probllin- ems of First Kind,” International Journal of Computer Mathematic s, Vol. 61, No. 1-2, 1996, pp. 155-163. doi:10.1080/00207169608804507 R. K. Mohanty and P. K. Pandey, “Families of Acc urate Discretizations of Order Two and Four for 3D Mildly Nonlinear Biharmonic Problems of Second Kind,” Inter- national Journal of Computer Mathematics, Vol. 68, No. 3-4, 1998, pp. 363-380. doi:10.1080/00207169808804702  D. J. Evans and R. K. Mohanty, “Block Iterative Methods for the Numerical Solution of Two-Dimensional Nonlin- ear Biharmonic Equations,” International Journal of C om- puter Mathematics, Vol. 69, No. 3-4, 1998, pp. 371-390. doi:10.1080/00207169808804729  R. K. Mohanty, D. J. Evans and P. K. Pandey, “Block Ite- rative Methods for the Numerical Solution of Three Di-First Kind,” International Journal of Computer Mathe- matics, Vol. 77, No. 2, 2001, pp. 319-332. doi:10.1 010880506 mensional Mildly Nonlinear Biharmonic Problems of 080/0020716 8  S. Singh and R. K. A New Cou- pled Ap Accuracy Nethod for the Solutioninear Biharions,” Neural Parallel and, 2009, pp. 239-256] D. Khattar, S. Singh and R. K. A New Cou- pled Approachcy Numethod for the Solution-Linear Bihuations,” Ap- plied Mnd Compu 215, No. 8, 2009, pp. 3036-3044. doi:10.1016/j.amc.2009.09.052, D. KhattarMohanty, “proach High umerical Mof 2D Nonl Scientific Computationsmonic Equat, Vol. 17. [14 Mohanty, “ High Accuraerical M of 3D Nonathematics aarmonic Eqtations, Vol. w High Accuracy Finite Difference Discretization for the Solution of 2D Non-Linear Bihar- monic Equations Using Coupled Approach,” Numerical Methods for Partial Differential Equations, Vol. 26, No. 4, 2010, pp. 931-944. doi:10.1002/num.204605 R. K. Mohanty, “A Ne Mohanty, M. K. Ja. Mishra, “A Com- t Discretization of O(h) for Two-Dimensional Non- ear Triharmonic Equations,” Physica Scripta, Vol. 84, . 2, 2011, pp. 025002. /0031-8949/84/02/025002 R. K.pacin and B. N4linNodoi:10.1088 K. Mohanty and S. Deyell Fourth Order Dif- ence Approximations for (u/x), (u/y) and (u/z) the Three Dimensional Quasi-Linear Elliptic Equation,” merical Methods for fferential Equations, , No. 5, 2000, pp. 417-425.  R. , “Single Cferof NuVol. 16Partial Didoi:10.1002/1098-2426(200009)16:5<417::AID-NUM1>3.0.CO;2-S  R. K. Mohanty, S. Karaa and U. Arora, “Fourth Order Nine artial Differential Equations,” Nume- Point Unequal Mesh Discretization for the Solution of 2D Non-Linear Elliptic Partial Differential Equations,” Neu- ral Parallel and Scientific Computations, Vol. 14, 2006, pp. 453-470.  R. K. Mohanty and S. Singh, “A New Highly Accurate Dis-cretization for Three Dimensional Singularly PerturbedNon-linear Elliptic Prical Methods for Partial Differential Equations, Vol. 22, No. 6, 2006, pp. 1379-1395. doi:10.1002/num.20160  L. A. Hageman and D. M. Young, “Applied Iterative Me- thods,” Dover Publications, New York, 2004.  M. K. Jain, “Numerical Solution of Differential Equa- tions,” 2nd Edition, John Wiley, New Delhi, 1984.  C. T. Kelly, “Iterative Methods for Linear and Non-Linear Equations,” SIAM Publications, Philadelphia, 1995.  Y. Saad, “Iterative Methods for Sparse Linear Systems,” SIAM Publications, Philadelphia, 2003. doi:10.1137/1.9780898718003  G. Meurant, “Computer Solution of Large Linear Systems,” North-Holland, Amsterdam, 1999.  W. F. Spotz and G. F. Carey, “High Order Compact Scheme for the Steady Stream-Function Vorticity Equations,” In- ternational Journal for Numerical Methods in Enginee- ring, Vol. 38, No. 20, 1995, pp. 3497-3512. doi:10.1002/nme.1620382008 Copyright © 2011 SciRes. AJCM