On Solution of Regular Singular Ordinary Boundary Value Problem On Solution of Nonlinear Singular Boundary Value Problem Luma N. M. Tawfiq Hiba W. Rasheed Dept. of Mathematics/College of Education for Pure Science (Ibn Al-Haitham,)/ Baghdad University . Received in : 30 September 2012 , Accepted in : 3 February 2013 Abstract This paper is devoted to the analysis of nonlinear singular boundary value problems for ordinary differential equations with a singularity of the different kind. We propose semi - analytic technique using two point osculatory interpolation to construct polynomial solution, and discussion behavior of the solution in the neighborhood of the singular points and its numerical approximation. Two examples are presented to demonstrate the applicability and efficiency of the methods. Finally, we discuss behavior of the solution in the neighborhood of the singularity point which appears to perform satisfactorily for singular problems. Kay ward : Singular boundary value problems, ODE, BVP. 331 | Mathematics @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ÚÓ‘Ój�n€a@Î@Úœäñ€a@‚Ï‹»‹€@·rÓ:a@Âig@Ú‹©@Ü‹1a26@@ÖÜ»€a@I3@‚b«@H2013 Ibn Al-Haitham Jour. for Pure & Appl. Sci. Vol. 26 (3) 2013 1. Introduction In the study of nonlinear phenomena in physics, engineering and other sciences, many mathematical models lead to singular two-point boundary value problems (SBVP) associated with nonlinear second order ordinary differential equations (ODE). In mathematics, a singularity is in general a point at which a given mathematical object is not defined, or a point of an exceptional set where it fails to be well-behaved in some particular way, such as many problems in varied fields as thermodynamics, electrostatics, physics, and statistics give rise to ordinary differential equations of the form : y″ = f(x, y, y') , a < x < b , (1) On some interval of the real line with some boundary conditions(BC). A two-point BVP associated to the second order differential equation (1) is singular if one of the following situations occurs: a and/or b are infinite; f is unbounded at some x0∈ [0,1] or f is unbounded at some particular value of y or y′ [1] . There are two types of a point x0 ∈ [0,1]: Ordinary Point and Singular Point. Also, there are two types of Singular Point : Regular and Irregular Points, A function y(x) is analytic at x0 if it has a power series expansion at x0 that converges to y(x) on an open interval containing x0. A point x0 is an ordinary point of the ODE (2), if the functions P(x) and Q(x) are analytic at x0. Otherwise x0 is a singular point of the ODE, i.e. y′′ + P(x)y′ + Q(x)y = 0 , (2) P(x) = P0 + P1(x-x0) + P2(x-x0)2+…….. = i i i xx )( 0 0 −Ρ∑ ∞ = , (3) Q(x) = Q0 + Q1(x-x0) + Q2(x-x0)2+………. = i i i xxQ )( 0 0 −∑ ∞ = , (4) On the other hand, if P(x) or Q(x) are not analytic at x0 then x0 is said to be a singular [2]. There are four kinds of singularities : • The first kind is the singularity at one of the ends of the interval [0,1] ; • The second kind is the singularity at both ends of the interval [0,1] • The third kind is the case of a singularity in the interior of the interval; • The forth and final kind is simply treating the case of a regular differential equation on an infinite interval. In this paper, we focus on the first kind. 2. Solution of Second Order Nonlinear SBVP In this section , we suggest semi analytic technique to solve second order nonlinear SBVP as following, we consider the SBVP : xm y''+ f( x, y, y' ) = 0 , (5a) gi( y(0), y(1), y'(0), y'(1) ) = 0 , i = 1 , 2 , m ϵ N , (5b) where f , g1, g2 are in general nonlinear functions of their arguments . The simple idea behind the use of two-point polynomials is to replace y(x) in problem (5), or an alternative formulation of it, by a P2n+1 which enables any unknown boundary values or derivatives of y(x) to be computed . Therefore, the first step is to construct the P2n+1 and to do so, we need the Taylor coefficients of y (x) at x = 0 : 332 | Mathematics @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ÚÓ‘Ój�n€a@Î@Úœäñ€a@‚Ï‹»‹€@·rÓ:a@Âig@Ú‹©@Ü‹1a26@@ÖÜ»€a@I3@‚b«@H2013 Ibn Al-Haitham Jour. for Pure & Appl. Sci. Vol. 26 (3) 2013 http://en.wikipedia.org/wiki/Mathematics http://en.wikipedia.org/wiki/Set_(mathematics) http://en.wikipedia.org/wiki/Well-behaved y = a 0 + a 1 x + ∑ ∞ =2i a i x i , (6) where y(0)= a0, y'(0)= a1, y"(0) / 2! =a2 ,…, y(i)(0) / i! = ai , i= 3, 4,… then insert the series forms (6) into (5a) and equate coefficients of powers of x to obtain a2 . Also we need Taylor coefficient of y(x) about x = 1: y = b 0 + b 1 (x-1) + ∑ ∞ =2i b i (x-1) i , (7) where y(1) = b0, y'(1) =b1, y"(1) / 2! =b2, … , y(i)(1) / i! =bi , i = 3,4,… then insert the series form (7) into (5a) and equate coefficients of powers of (x-1) to obtain b2, then derive equation (5a) with respect to x and iterate the above process to obtain a3 and b3 , now iterate the above process many times to obtain a4, b4, then a5, b5 and so on, that is, we can get ai and bi for all i ≥ 2( the resulting equations can be solved using MATLAB to obtain ai and bi for all i ≥ 2 ), the notation implies that the coefficients depend only on the indicated unknowns a0, a1, b0, b1, we get two of these four unknown by the boundary condition. Now, we can construct a P2n+1(x) from these coefficients ( aisۥ and bisۥ ) by the following : P2n+1 = ∑ = n i 0 { ai Qi(x) + (-1)i bi Qi(1-x) } , (8) where ( x j / j!)(1-x) 1+n ∑ − = jn s 0       + s sn xs = Q j (x) / j! it can be seen that (8) have only two unknowns from a0, b0, a1 and b1 to find this, we integrate equation (5a) on [0, x] to obtain : xmy'(x) – mxm-1y(x) + m(m–1) ∫ x 0 xm-2y(x) dx + ∫ x 0 f(x, y, y') dx = 0 , (9a) and again integrate equation (9a) on [0, x] to obtain: xmy(x) –2m ∫ x 0 xm-1y(x) dx +m(m-1)∫ x 0 (1-x)xm-2y(x)dx+∫ x 0 (1-x)f(x,y,y') = 0 , (9b) Putting x = 1 in (9) then gives : b1 – mb0 + m(m-1) ∫ 1 0 xm-2 y(x) dx + ∫ 1 0 f(x, y, y') dx = 0 , (10a) and b0–2m∫ 1 0 xm-1y(x)dx +m(m-1)∫ 1 0 (1-x)xm-2 y(x) dx + ∫ 1 0 (1-x)f(x, y, y')dx = 0 , (10b) Use P2n+1 as a replacement of y(x) in ( 10 ) and substitute the boundary conditions (5b) in (10) then, we have only two unknown coefficients b1, b0 and two equations (10) so, we can find b1, b0 for any n by solving this system of algebraic equations using MATLAB, so insert b0 and b1 into (8), thus (8) represents the solution of (5). Extensive computations have shown that this generally provides a more accurate polynomial representation for a given n . 3. Examples In this section, we introduce two examples of second order SBVP, non homogenous, nonlinear ordinary differential equations with two point BC to assess the performance of the proposed method. Example 1 : Emden's equation Emden's equation arises in modeling a spherical body of gas. The PDE of the model is reduced by symmetry to the ODE : 333 | Mathematics @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ÚÓ‘Ój�n€a@Î@Úœäñ€a@‚Ï‹»‹€@·rÓ:a@Âig@Ú‹©@Ü‹1a26@@ÖÜ»€a@I3@‚b«@H2013 Ibn Al-Haitham Jour. for Pure & Appl. Sci. Vol. 26 (3) 2013 x ∈ [0, 1] , the coefficient 2/x is singular at x = 0, but symmetry implies the BC: y'(0) = 0. With this boundary condition, the term (2/x)y'(0) is well-defined as x approaches 0. For the boundary condition , y(1) = 2 3 , this nonlinear SBVP has the analytical solution : We solve this problem by using semi-analytic technique by the following; from equation (8), we have: P9 = –0.0005654571561x9 +0.001454732112x8+0.00321137204x7 –0.01349180792 x6 + 0.0004165633765x5 + 0.04166666759x4 – 0.1666666687x2 + 1.000000002 Now, increase n, to get higher accuracy ,let n = 5, 6,respectively ,i.e., P11 = 0.000030043522x11+0.0002663457696x10 –0.00234783883x9+0.00531993067 x8 – 0.0008040441181x7–0.01143903456x6+0.04166666759x4–0.1666666687x2+1.000000002 P13 = 0.00002406664878 x13 – 0.0002794397776 x12 + 0.001201282562 x11 – 0.00224058548x10 +0.0007139586051x9 + 0.003149932497x8 + 0.00003026184015x7 – 0.01157407444x6 + 0.04166666759x4 – 0.1666666687x2 + 1.000000002 The standard bvp4c syntax was implemented in [3] to solve the current problem, the example evaluates the analytical solution at 100 equally-spaced points and plots it along with the numerical solution computed using bvp4c and the results given in Table 1 with the results of suggested method for n = 4, 5, 6, i.e., P5, P7, P9 and the results of [3] . Also Table 2 gives the accuracy of suggested method and result of [3]. Figures 1, 2, 3, illustrate Emden problem, and suggested method for n = 4, 6 respectively . It can be seen from Table2 that the maximum error of method in [3] is 0.0000018522, while the maximum error of suggested method is 2.763537509942182e-009 Example 2 Consider the following nonlinear SBVP : , b ϵ R , 0 ≤ x ≤ 1 5 53 4 )45(5 x bxexx y + −−− = y')1+ x 1 + ( y" with BC : y'(0) = 0 , y'(1) = -1. The exact solution for this problem is: y = - ln( x5 + 4) This problem is an application of oxygen diffusion , we solve this problem by suggested method and we take b=1, by applying equation (8) we have ( for n = 7 ) : P15 = - 0.1142318896 x15 + 0.7852891723 x14 - 2.23317271 x13 + 3.413936838 x12 - 3.087794733x11+1.697854523x10-0.4991720457x9+ 0.06414729306x8 - 0.25x5 -1.386294361. For more details ,table(3) gives the results for different nodes in the domain, for n = 7 and figures (4) illustrate suggested method for n = 7. Abukhaled et al in [4] applying L̛ Hopital ̛̛s rule to overcome the singularity at x = 0 and then the modified spline approach is used and got maximum error 7.79e-4 and resolution this problem using finite difference method then gives maximum error 1.46e-3, from Table 3 we have maximum error 9.399395723974635e-007 Therefore, the proposed method provides superiority results. 4. Behavior of the solution in the neighborhood of the singularity x= 0 Our main concern in this section will be to study the behavior of the solution in the neighborhood of singular point . Consider the following SIVP : y′′(x) + ((N − 1) / x ) y′(x) = f(y) , N ≥ 1 , 0 < x < 1 , (11) y(0) = y0 , limx→0+ x y′(x) = 0 , (12) 334 | Mathematics @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ÚÓ‘Ój�n€a@Î@Úœäñ€a@‚Ï‹»‹€@·rÓ:a@Âig@Ú‹©@Ü‹1a26@@ÖÜ»€a@I3@‚b«@H2013 Ibn Al-Haitham Jour. for Pure & Appl. Sci. Vol. 26 (3) 2013 where f(y) is continuous function . As the same manner in [5], let us look for a solution of this problem in the form : y(x) = y0 − C xk (1 + o(1)) , (13) y′(x) = − C k xk−1(1 + o(1)) , y′′(x) = − C k (k − 1) xk−2(1 + o(1)) , x → 0+ where C is a positive constant and k > 1. If we substitute (13) in (11) we obtain : C = (1/ k) (f(y0) /N )k−1 , (14) In order to improve representation (13) we perform the variable substitution : y(x) = y0 − C xk (1 + g(x)) , (15) we easily obtain the following result which is similar to the results in [5]. Theorem 2 For each y0 > 0, problem (11), (12) has a unique solution in the neighborhood of x = 0 that can be represented by: y(x, y0) = y0 − C xk (1 + g xk + o(xk) ) , where k, C and g are given by (14) and (15), respectively. We see that these results are in good agreement with the ones obtained by the method in [5], they are also consistent with the results presented in [4]. In order to estimate the convergence order of the suggested method at x = 0, we have carried out several experiments with different values of n and used the formula : cy0 = − log2 ( |y0n3 − y0n2| / |y0n2− y0n1 | ) , (16) where y0ni is the approximate value of y0 obtained with ni = 1,2, 3, 4,… Now, we apply the formula in equation (16) to the example 1, as following : Let y0i is the approximate value of y0 evaluated by suggested method with n = i, i = 2, 3, 4, 5 ,6 . y0i P2n+1 1.000048243794667 P5 1.000009487495905 P7 1.000000801904571 P9 1.0000000507769192 P11 1.0000000024600459 P13 First ,we take i= 2, 3, 4, i.e., || || log 0203 0304 20 yy yy C y − − −= 005 006 20 043788756298762.3 45396855133393.8 log − − −= e e C y = 2.157734820323959 . The value of Cy0 illustrates that the convergence order estimate of this case is close to 2. Now, if we take i= 3, 4, 5, i.e., Cy0 || || log 0304 0405 2 yy yy − − −= 006 007 2 45396855133393.8 458845112765185.7 log − − −= e e = 3.531494058780649. The value of Cy0 illustrates that the convergence order estimate of this case is close to 4. Now if we take i= 4, 5, 6 , i.e., || || log 0405 0506 20 yy yy C y − − −= 007 008 2 458845112765185.7 085288316873169.4 log − − −= e e , =3.958459111546360 The value of Cy0 illustrate that the convergence order estimate of this case is close to 4, and so on, if we increase n, we see that the order of convergence also increases. 335 | Mathematics @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ÚÓ‘Ój�n€a@Î@Úœäñ€a@‚Ï‹»‹€@·rÓ:a@Âig@Ú‹©@Ü‹1a26@@ÖÜ»€a@I3@‚b«@H2013 Ibn Al-Haitham Jour. for Pure & Appl. Sci. Vol. 26 (3) 2013 References [1] Robert, L. B., and Courtney, S. C. (1996) " Differential Equations A Modeling perspective" , United States of America. [2] Rachůnková, I., Staněk, S., and Tvrdý, M. (2008) " Solvability of Nonlinear Singular Problems for Ordinary Differential Equations", New York, USA. [3] Hamdan, S. M. (2010) "Singular Two Points Boundary Value Problem", MSc thesis in Computational Mathematics, An-Najah National University, Faculty of Graduate Studies, Nablus, Palestine. [4] Abukhaled, M., Khuri, S. A. and Saufy, A. (2011) "A Numerical Approach for Solving A Class of Singular Boundary Value Problems Arising in Physiology", International Journal of Numerical Analysis and Modeling, Vol.8, No.2, PP:353– 363. [5] Morgado, L., and Lima, P. (2009) " Numerical methods for a singular boundary value problem with application to a heat conduction model in the human head", Proceedings of the International Conference on Computational and Mathematical Methods in Science and Engineering, CMMSE. Table 1:The result of the method for n = 4,5,6, and result in [3] for Example1 Table 2:The accuracy of the method for n= 6 ,i.e, P13 and result in [3] for Example1 Table 3 :The exact and suggested solution for n = 7 of Example 2 xi exact solution y(x) y1(x) using numerical solution Osculatory interpolation P9 Osculatory interpolation P11 Osculatory interpolation P13 0.000 1.000000000000000 1.0000005109 1.000000801904571 1.000000050776919 1.000000002460046 0.125 0.997405961908059 0.9974064705 0.997406758716985 0.997406012221857 0.997405964341573 0.250 0.989743318610787 0.9897438277 0.989744154146833 0.989743371348405 0.989743321151803 0.500 0.960768922830523 0.9607707751 0.960769739668355 0.960768977298348 0.960768925594060 0.750 0.917662935482247 0.9176644937 0.917663076945894 0.917662942247223 0.917662935723868 1.000 0.866025403784439 0.8660254037 0.866025403784439 0.866025403784439 0.866025403784439 xi exact solution y(x) Error│y(x) - y1(x) │ Error │ y(x) - P13│ 0.000 1.000000000000000 0.0000005109 2.460045944729927e-009 0.125 0.997405961908059 0.0000005086 2.433513945909738e-009 0.250 0.989743318610787 0.0000005091 2.541016064228074e-009 0.500 0.960768922830523 0.0000018522 2.763537509942182e-009 0.750 0.917662935482247 0.0000015582 2.416213895628516e-010 1.000 0.866025403784439 0.0000000000 0.0000000000 S.S.E=2.612609927853621e-017 P15 1.386294361119891- ao 1.6094379124341- bo Errors | y(x) - P15 | P15 exact solution y(x) 4.440892098500626e-016 1.386294361119891- 1.386294361119891- 0 2.812807764485115e-010 -1.386296860835485 1.386296861116765- 0.1 2.834240797611187e-008 1.386374329577653- 1.386374357920061- 0.2 336 | Mathematics @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ÚÓ‘Ój�n€a@Î@Úœäñ€a@‚Ï‹»‹€@·rÓ:a@Âig@Ú‹©@Ü‹1a26@@ÖÜ»€a@I3@‚b«@H2013 Ibn Al-Haitham Jour. for Pure & Appl. Sci. Vol. 26 (3) 2013 Figure 1: Emden problem – SBVP Figure 2: Comparison between the exact and semi-analytic solution P9 of example1 Figure 3 :Comparison between the exact and semi-analytic solution P13 of example1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1 1.02 1.04 the solution at n=4 x y y P9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.86 0.88 0.9 0.92 0.94 0.96 0.98 1 1.02 1.04 the solution at n=6 x y exact P13 2.489597183963355e-007 -1.386901427706747 1.386901676666466- 0.3 7.099738041915771e-007 -1.388850379927776 -1.388851089901581 0.4 9.399395723974635e-007 1.394075561622373- 1.394076501561946- 0.5 6.23686427614345e-007 -1.405547194355350 -1.405547818041777 0.6 1.862599887658689e-007 1.427452912675769- 1.427453098935757- 0.7 1.679306937951708e-008 -1.465031584864205 -1.465031601657275 0.8 1.136124527789661e-010 -1.523986772073695 1.523986772187307- 0.9 2.220446049250313e-016 1.609437912434100- -1.609437912434100 1 S.S.E = 1.874293078482109e-012 337 | Mathematics @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ÚÓ‘Ój�n€a@Î@Úœäñ€a@‚Ï‹»‹€@·rÓ:a@Âig@Ú‹©@Ü‹1a26@@ÖÜ»€a@I3@‚b«@H2013 Ibn Al-Haitham Jour. for Pure & Appl. Sci. Vol. 26 (3) 2013 Figure 4:Comparison between the exact and suggested solution P15 of example2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 -1.65 -1.6 -1.55 -1.5 -1.45 -1.4 -1.35 -1.3 the solution at n=7 x y p15 exact 338 | Mathematics @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ÚÓ‘Ój�n€a@Î@Úœäñ€a@‚Ï‹»‹€@·rÓ:a@Âig@Ú‹©@Ü‹1a26@@ÖÜ»€a@I3@‚b«@H2013 Ibn Al-Haitham Jour. for Pure & Appl. Sci. Vol. 26 (3) 2013 غیر الخطیة مسائل القیم الحدودیة الشاذةحول حل لمى ناجي محمد توفیق ھبة ولید رشید جامعة بغداد / )ابن الھیثم للعلوم الصرفة ( كلیة التربیة/ قسم الریاضیات 2013شباط 3في ، قبل 2012أیلول 30أستلم البحث في : المستخلص للمع�ادالت التفاض�لیة غی�ر الخطی�ة ع�رض دراس�ة تحلیلی�ة لمس�ائل الق�یم الحدودی�ة الش�اذة بح�ثال االھدف من ھذ النقطت�ین للحص�ول عل�ى الح�ل يالتقنیة شبھ التحلیلیة باس�تخدام االن�دراج التماس�ي ذننا نقترح إذ إاالعتیادیة وبأنواع مختلفة وسھولة أداء الطریق�ة المقترح�ة و أخی�را ناقش�نا س�لوك الح�ل ف�ي والكفایةلتوضیح الدقة ، مثالینكمتعددة حدود، كذلك ناقشنا إیجاد الحل التقریبي لھا بشكل مرضي فیما یخص المسائل الشاذة .، وجوار النقاط الشاذة 339 | Mathematics @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ÚÓ‘Ój�n€a@Î@Úœäñ€a@‚Ï‹»‹€@·rÓ:a@Âig@Ú‹©@Ü‹1a26@@ÖÜ»€a@I3@‚b«@H2013 Ibn Al-Haitham Jour. for Pure & Appl. Sci. Vol. 26 (3) 2013 Figure 2: Comparison between the exact and semi-analytic solution P9 of example1 Figure 3 :Comparison between the exact and semi-analytic solution P13 of example1 Figure 4:Comparison between the exact and suggested solution P15 of example2