Second International Conference Modelling and Development of Intelligent Systems Sibiu – Romania, September 29 – October 02, 2011 Several methods of… [626147]

Second International Conference
Modelling and Development of Intelligent Systems
Sibiu – Romania, September 29 – October 02, 2011
Several methods of approxi mation for second order nonlinear
boundary value proble mwith boundary conditions at infinity
Daniel N. Pop, Radu T. Trˆ ımbit¸a¸s
Abstract
Consider the problem:/braceleftbiggy/prime/prime(x)+f(x, y)=0 ,0<x<∞,
y(0) =∞,y(∞)=0
where f(x, y)∈C([0,∞]×R),y(x)∈C1(0,∞).This is not a classical two-points boundary value
problem since y(0) =∞,y(∞) = 0. To solve this kind of proble ms we need to know the values in two
inner points a, b∈(0,∞),a/negationslash=b.The aim of this work is to present three approxi mation procedures:
1. A com binedmethod using collocation method on B-splines of order ( k+2) witha( k+1)order
Runge-Kutta method.
2. A pseudospectral collocation method with Chebychev extreme points co mbined with a Runge-
Kuttamethod.
3. MATLAB function bvp4ccombined with a Runge-Kutta method.
Then we give a nu merical exam ples and com pare the costs (ti me U.C) using MATLABfunctions
tic-toc.
1 Introduction
Consider the problem (PVP):
y/prime/prime(x)+f(x,y)=0,x∈(0,∞)( 1)
y(a)=α (2)
y(b)=β, a,b ∈(0,∞),a<b . (3)
wheref(x,y)∈C((0,∞)×R),a,b,α,β ∈R.
We try to solve this problem using three approxi mation methods:
1. A com bined method based on collocation with B-splines of order ( k+2)anda Runge-Kutta method
order (k+1).
2. A pseudospectral collocation method with Tchebychev extreme points co mbined with a Runge-
Kutta method.
3. MA TLAB function bvp4ccombined with a Runge-Kutta method.
121

/g3

Several methods of approximation for second order nonlinear boundary value problem with boundary
conditions at infinity
/g3/g3
/g3

The methods are new in this context since the conditions are stated at the interior points of the interval
(0,1). Proble ms of type (1)+(2)+ (3) occurs in practice. Exa mples are in sem iclassical description of
the charge density in ato mso fh i g ha t o mic nu mber ( Thomas-Fermi equation) [19, pp.155-156], reaction-
diffusion equation [9], frequency do main equation for the vibrating string (Greengard-Rokhlin proble m)
[11], electro magnetic self interaction theory [4, pp.336-337], the model of the steady concentration of a
substrate in an enzy me-catalyzed reaction (Michaelis-Men ten kinetics)[19, page 145].
We also consider the problem (BVP):
y/prime/prime(x)+f(x,y)=0,x∈[a,b]( 4)
y(a)=α (5)
y(b)=β, (6)
Also it is shown that the Runge-Kutta method does not degrade the accuracy provided by the col-
location method for the (BVP) proble m[4,Theorem5.73 pp219, Theorem5.140 pp253]. To apply the
collocation theory we need to have an isolated solution of (BVP) problem and this occurs if the above
linearized problem fory(x) is uniquely solvable.
R.D. Russel andL.F. Shampine [18] study the existence and the uniqueness of the isolated solution.
Our methods consists into deco mposition of the proble m(1)+(2)+(3) into three proble ms:
1. A (BVP) proble mon [a,b]( p r o b l e m(4)+(5)+(6)).
2.Two (IVPs) on (0,a ]a n d[b,+∞).
For the existence and uniqueness of an (IVP), see [14, pp: 112-113].
If the problem (BVP) has the unique solution, the require menty(x)∈C1(0,+∞) ensure the existence
and the uniqueness of the solution of the problem (PVP). Our choice to use these methods is based on
the following reasons:
1. We write the code using the function spcolinMatlab-Spline Toolbox [15] and the functions
cebdif, cebint, cebdifft contained in dmsuite [21].
2.Theoretical results on the convergence of collocation method are given in ( [12], [13]).
3.The accuracy of spectral method is superior to finite elem ents method (FEM) and finite difference
methods (FDM) (the rate of convergence associated with this proble ms with s mooth conditions are
O(exp(−CN)) or O(exp(C2√
N)w h e r eNis the num ber of degrees of freedo min the expansions).
4. For each Newtoniteration,the resulting linear algebraic syste mof equations (after using Newton
method with quasilinearization) is solved using method given in [8].
2A c o mbinedmethod using B-splines and Runge-Kutta meth-
ods
First we solve the (BVP) proble musing the collocation method with B-splines of order ( k+2 )
presented in [17, Section 2].
Consider the mesh of [a,b]:
Δ:a=x0<x1<···<xN=b, (7)
where the multiplicity of aandbis (k+ 2) and the multiplicity of inner points is k.So the dimension
of spline space is n=Nk+2.Also we construct the collocation points ξj,j=1,2,…,n−2 like in [17,
Section 1] and [2].
We wish to find an approxim ate solution of the (BVP) problem , having the following for m:
122

/g3
/g3
Daniel N. Pop, Radu T. Trîmbi/g288a/g250
/g3
/g3

u
Δ(x)=n−1/summationdisplay
i=0ciBi,k+1(x), (8)
whereBi,k+1(x) is the B-spline of order ( k+2).
Our approximation method is inspired fro m[7, Chapter 2,5]. We impose the conditions:
c1The approxim ate solution (8) satisfies the differential equation (4) at collocation points:
ξj,j=1,2,…n−2.
c2The solution satisfies u
Δ(a)=α, u
Δ(b)=β.
The above conditions yield a nonlinear syste mwithnequations:
⎧
⎨
⎩/summationtextn−1
i=0ciBi,k+1(a)=α,/summationtextn−1
i=0ciB/prime/prime
i,k+1(ξj)+f(ξj,/summationtextn−1
i=0ciBi,k+1(ξj)) = 0,j=1,2,…,n−2,/summationtextn−1
i=0ciBi,k+1(b)=β,
with unknowns ci,i=0,…,n−1. IfF=[F0,F1,…,FN−1] are the functions defined by the equations
of the nonlinear system, using the quasilinearization of Newton method [4, pp: 52-55], we find the next
approximation by means of
c(k+1)=c(k)−w(k),
wherec(k)is the vector of unknowns obtained at the k−th step and w(k)is the solution of the linear
syste m
F/prime(c(k))w=F(c(k)).
To solve the (BVP) proble mwe use the method presented in [20] and the initial approximation
u(0)∈C1[0,1] is required. The successful stopping criterion [3] is:
/vextenddouble/vextenddouble/vextenddoubleu(k+1)−u(k)/vextenddouble/vextenddouble/vextenddouble≤abstol+/vextenddouble/vextenddouble/vextenddoubleu(k+1)/vextenddouble/vextenddouble/vextenddoublereltol,
whereabstolandreltolis the absolute and the relative error tolerance, respectively and the nor mis the
usual unifor mconvergencenor m.Thereliabilityofthe error-estimationprocedurebeing used forstopping
criterion was verified in [8]. For the solution of two IVPs on (0 ,a]a n d[b,+∞)w eu s ea Runge-Kutta
method of appropriate order, this need good approxi mation ofy/prime(a)a n dy/prime(b), which could be obtained
with noadditional effort during the collocation method.
The stability and convergence of Runge-Kutta method are guaranteed in [10, Theore m5.3.1 page
285, Theorem 5.3.2 page 288]. A ( k+ 1) order explicit Runge-Kutta method is consistent and stable,
so is convergent. The convergence and accuracy of our co mbined method to whole interval (0 ,+∞)w a s
proved in [17, Section 3, Theorem3.1] and the total costs of this method was studied in [17, Section 4].
3A c o mbinedmethod using a pseudospectral collocation with
Tchebychev extre me points and Runge-Kutta methods
Consider the grid:
Δ:0=x−q< … < x −1<a=x0<x1< … < x N=b<xN+1< … < x N+p. (9)
Oursecond method is a co mbined apseudospectral method forthe (BVP) proble mandaRunge-Kutta
method for the two IVPs on (0 ,a]a n d[b,+∞).The approxim ate solution of (BVP) problem follow the
ideas presented in [5]. Let y(x) of this problem and considering the Lagrange basis (l k)w eh a v e :
y(x)=N/summationdisplay
k=0lk(x)y(xk)+(RNy)(x),x∈[a,b]
123

/g3

Several methods of approximation for second order nonlinear boundary value problem with boundary
conditions at infinity
/g3/g3
/g3

where :
(RNy)(x)=y(N+1)(ξ)
(N+1)!(x−x0)…(x−xN)
is the remainder of Lagrange interpolation. Since y(x) fulfills the differential equation (4) we obtain:
N/summationdisplay
k=0l/prime/prime
k(x)y(xk)+(RNy)/prime/prime(x)=−f(xi,y(xi)),i=1,2,…,N−1.
Settingy(xk)=ykand ignoring the rest, one obtains the nonlinear syste m:
N/summationdisplay
k=0l/prime/prime
k(x)y(xk)=−f(xi,y(xi)),i=1,2,…,N−1, (10)
with unknowns yk,k=1,…,N−1,herey0=y(a)=αandyN=y(b)=β.The approxim ate solution
(that is the collocation polyno mial for (BVP) proble m), is the Lagrange interpolation polynom ial at
nodes{xk},k=0,1,2,…N:
yN(x)=N/summationdisplay
k=0lk(x)y(xk). (11)
The nonlinear syste m(10) can be rewritten as:
AYN=F(YN)+bN
where:
A=[aik],aik=l/prime/prime
k(xi),k , i=1,2,…,N−1,
F(YN)=⎡
⎢⎢⎢⎣−f(x
1,y1)
−f(x2,y2)
…
−f(xN−1,yN−1)⎤
⎥⎥⎥⎦,b
N=⎡
⎢⎢⎢⎣−αl/prime/prime
0(x1)−βl/prime/prime
N(x1)
−αl/prime/prime
0(x2)−βl/prime/prime
N(x2)
…
−αl/prime/prime
0(xN−1)−βl/prime/prime
N(xN−1)⎤
⎥⎥⎥⎦.
If the nodes {x
k},k=0,1,…Nare sy mmetric with respect of ( a+b)/2,Ais centro-symmetric [6, for
proof], so nonsingular. So we choose the nodes given by:
xi=(b−a)cosπi
N+b+a
2,i=1,2,…,N, (12)
i.e. the Chebyshev -Lobatto nodes. We introduce :
G(Y)=A−1F(Y)+A−1bN. (13)
Tos o l v en u merically (PVP) problem on Δ given by (9) we apply pseudo-spectral collocation method
at points [ a,b]a n daRunge-Kutta method to other points. To apply the Runge-Kutta method for the
solution of two (IVP) on (0 ,a]a n d[b,+∞) we need the derivatives y/prime(a)a n dy/prime(b),this can be computed
by deriving the form ula (11). In work [5], the authors prove the existence of unique solution of the syste m
(10) which can be calculated by successive approxi mation method:
Y(N+1)=G(Y(N)),n∈N∗,
withY(0)fixed and Ggiven by (13), also they esti mate the error:
/bardblY−YN/bardbl≤/vextenddouble/vextenddoubleA−1/vextenddouble/vextenddouble/bardblR/bardbl
1−/bardblA−1/bardblL
124

/g3
/g3
Daniel N. Pop, Radu T. Trîmbi/g288a/g250
/g3
/g3

where
Y=[y(x1),y(x2),…,y(xN−1)]T,
y(x) is the exact solution of (BVP) proble m,
YN=[y1,y2,…,yN−1]T,
yiare the values of approximated solution at xicomputed by (12),
R=[−(RNy)/prime/prime(x1),−(RNy)/prime/prime(x2),…,−(RNy)/prime/prime(xN−1)]T,
andLis theLipschitz constant. Co mbining these results with the stability and convergence of Runge-
Kutta methods in [16, Theorem2.3] the authors prove the convergence of this method and occurs:
|yN(a)−y/prime(a)|=O(hk),
|yN(b)−y/prime(b)|=O(hk),
and for each points xiin Δ giving by (9):
|yN(xi)−yi|=O(hk),i=−q,…,N+p.
If theRunge-Kutta method is stable and has the order k, then the final solution has the sam e accuracy.
4M A T L A B s o l v e r bvp4cand Runge-Kutta methods
MATLAB solver bvp4cis a strong solver based on collocation. It allows a flexible description of the
ODEs, various kind of boundary conditions, para meters and options (jacobian, tolerances, vectorization
and so on). It requires a guess solution. As a side effect it provides an approxi mation of the derivative
of the solution. This allows us to co mbine bvp4cwith an IVP solver.
5N umerical exam ples
For the both methods we im plemented the ideas in MATLAB 2010a [15], for the first method our
code use Matlab Spline Toolbox,thefunction spcolallowsus to co mpute easily the collocation matrix
and for (IVP) proble ms the solver ode23tb works fine (when the problem is stiff). To avoid the error
propagation, we choose for (BVP) proble mB-splines of order 4 (degree 3) or order 5 (degree 4), in this
we implemented the function polycalnlinRK.
The second method was implemented in MATLABusing the functions cebdif,cebint and cebdiff
contained in dmsuite and described in [21], we write the function solvepolylocalceb who solve the
nonlinear system and call the Runge-Kutta solver ode23tb. The derivatives at aandbwere com puted
by calling cebdifft. The third method following the idea given by [19].
In order to co mpare the costs (run-times) experi mentally we use Matlabfunctions ticandtoc.
•T h efi r s te x a mple is the case where we know the exact solution
Consider the (BVP) proble m:
⎧
⎨
⎩y/prime/prime(x)+2π2exp(−y(x)) = 0,0<x<∞,
y(1/10) = 2lnsin π/10,
y(9/10) = 2lnsin9 π/1 0=2l ns i n π/10.(14)
The exact solution of this proble mis:
y(x)=2l n ( s i n ( πx)),
125

/g3

Several methods of approximation for second order nonlinear boundary value problem with boundary
conditions at infinity
/g3/g3
/g3

and we see that the exact solution is a periodic function of the period T= 1. Using step control
algorith m[1] we deter mined that problem (14) has singularities in x=0a n di n x=1.We
have set the tolerance to ε=1 0−10,we took N= 1025 and maximumnumber of iterations
NMAX =5 0.The start solution are obtained using the Lagrange interpolation polyno mial with
nodes: 1/4,5/24,1/6.The results we have obtained after 10 iterations, run ti mes are:
Tolerance
1st Method
2nd Method
3rd Method
10−10
0.777144
7.204652
1.350691
The graph of approxi mate solution are presented in Figure 1, and the errors in sem i-logarith mic
scale for the first method in Figure 2(a), for the second in Figure 2(b) and for the third in Figure
2(c), respectively.
0 0.2 0.4 0.6 0.8 1−12−10−8−6−4−20Approximate solution
solution
d=1/10,e=9/10
Figure 1: Approx-solution
0 0.2 0.4 0.6 0.8 110−1610−1410−1210−1010−810−610−410−2
xln |y(x)−yΔ(x)|Error in semilogarithmic scale
(a) B-splines+Runge-Kutta0 0.2 0.4 0.6 0.8 110−1610−1410−1210−1010−810−610−410−2
xln |y(x)−yΔ(x)|Error in semilogarithmic scale
(b) C.C+Runge-Kutta0 0.2 0.4 0.6 0.8 110−1610−1410−1210−1010−810−610−410−2
xln |y(x)−yΔ(x)|Error in semilogarithmic scale
(c)bvp4c+Runge-Kutta
Figure 2: Errors
126

/g3
/g3
Daniel N. Pop, Radu T. Trîmbi/g288a/g250
/g3
/g3

•The second exa mple is the case of unknown exact solution
y/prime/prime(x)=x−1/2y3/2, (15)
with boundary conditions:
y(0) = +∞,y(∞)=0. (16)
This (BVP) arises in a sem iclassical description of the charge density in ato ms of high ato mic
number. There are difficulties at both end points. These difficulties are discussed at length in
Davis(1962) and in Bender and Orszag (1999).Davisdiscusses series solutions for:
y(x)a sx→0.
It is clear that there are fractional powers in the series. That is because, with y(0) = 1,ODE
requires:
y/prime/prime(x)∼x−1/2asx→0,
and hence there be a term4
3×3/2in series for y(x).Of course, there must also be lower-order term s
so as to satisfy the boundary condition at x=0.Bender and Orszag discuss the asy mptotic
behavior of y(x),x→0.Verify that trying a solution of the for m
y(x)∼axα,
yields for the start solution:
y0(x) = 144x−3.
We use for inner points a=0.015,andb=5 9.The results we have obtained after 12 iterations, run
times are:
Tolerance
1st Method
2nd Method
3rd Method
10−10
0.493448
1.204652
0.838735
The graph of approxi mate nonlinear solution of Fermi-Thomas proble mis presented in Figure 3.
0 10 20 30 40 50 60 7000.10.20.30.40.50.60.70.80.91Fermi−Thomas equation BS+RK
Figure 3: The charge density in ato mso fh i g ha t o mic nu mber
127

/g3

Several methods of approximation for second order nonlinear boundary value problem with boundary
conditions at infinity
/g3/g3
/g3

6 Conclusions
The running tim e for B-spline collocation is the shortest, because its collocation matrix is banded.
Tchebychev collocation has the longest ti me, since its collocation matrix is full. The bvp4csolver has an
inter mediary position, since it has a different type collocation. Nevertheless, further tests are necessary.
References
[1] O.Agratini, I.Chiorean, G.Co man, R. Trˆımbit¸a¸s,Analiz˘an u m e r i c ˘ a¸si Teoria Aproxim˘ arii, Presa
Universitar˘ a Clujean˘ a, Cluj-Napoca, 2002.
[2] U.Ascher, S.Pruess and R.D. Russel, On spline basis selection for solving differential equations,
SIAM J. Numer. Anal. 20 nr:1 , 1983.
[3] U.Ascher, Solving boundary value proble ms with spline collocation code, Journal of co mputational
physics, 34 , Nr:3, pp: 401-415, 1980.
[4] U. Ascher, R.M. Mattheij, and R.D. Russel, Numerical Solution of Boundary Value Problems for
Ordinary Differential Equations, SIAM, 1997.
[5] F.A.Costabile, A.Napoli, A coolocation method for global approximation of general second order
BVPs,Computing Letters 3, pp: 23-24, 2007.
[6] F.A.Costabile, F.Dell-Accio, On the nonsingularity of a special class of centrosymmetric matrices
arising in spectral method in BVPs, Applied Mathematics and Computation , 206, pp: 991-993,2008.
[7] C.de Boor, A practical guide to splines , Springer-Verlag, Berlin, Heidelberg, New-York, 1978.
[8] C.de Boor, R.Weiss, Solveblock: A package for solving almost block diagonal linear systems, with
applications on the numerical solution of odinary differential equations ,M C R # TSR 1625, Madi-
son,Wisconsin, 1976.
[9] M.A. Di minescu, V. Nistorean, Gospod˘arirea durabil˘ a a resurselor de ap˘ a,Conf. Internat ¸. Energie-
Mediu,CIEM , Bucure¸ sti, 2005.
[10] W.Gautschi, Numerical Analysis, An Introduction , Birkhauser Verlag, Basel, 1997.
[11] L.Greengard, V.Rohlin On the Numerical Solution of Two-Point Boundary Value Problems ,
Comm.Pure Appl.Math.XLIV , pp:419-452,1991.
[12] W.Heinrichs, Strong convergence estimates for pseudospectral methods,Applications of Mathemat-
ics, 37(6), pp:24-37, 1992.
[13] E.Houstis, A collocation method for syste ms of nonlinear ordinary differential equations, J. Math.
Anal. Appl. , 62, pp:24-37, 1978.
[14] P. Henrici, Discrete Variable Methods in Ordinary Differential Equations , John Willey and Sons,
New York, 1962.
[15] Mathlab 2010a, Set of manuals, , MathWorks, Inc. , Natick, MA, 2010.
[16] Ion P˘ av˘aloiu, Daniel N.Pop, Radu T.Trˆımbit¸a¸s, Solution of polylocal problem with a pseodospec-
tral method,Annals of ”Tiberiu Popoviciu” Seminar of functional equations, approximation and
convexity, vol 8, pp: 53-65, Cluj-Napoca, 2010.
[17] Daniel N.Pop, Radu T.Trˆımbit¸a¸s, An approximation methods for second order nonlinear value
polylocal proble ms using B-splines and Runge-Kutta metods, N.A.A. T. 2010, Studia Universitas
”Babes-Bolyai”, Series Mathematica , Vol.LVI, Nu mber 2, pp:515-527, June 2011.
[18] R.D.Russel, L.F.Sha mpine, A Collocation Method for Boundary Value Proble ms,Numer. Math. ,
19, Springer-Verlag, pp: 1-28, 1972.
[19] L. F. Sham pine, I. Gladwell, S. Thompson, Solving ODEs with MA TLAB, Ca mbridge University
Press, 2003.
[20] R.Stover, Collocation methods for solving linear-algeb raic boundary value proble m,Numer.Math.
88, pp:771-795, 2001.
[21] J.A.C Weidem an, S.C Reddy, A MA TLAB Differentiation Matrix Suite, ACM Trans. on Math.
Software, 26, pp: 465-519, 2000.
128

/g3
/g3
Daniel N. Pop, Radu T. Trîmbi/g288a/g250
/g3
/g3

Daniel N.Pop Radu T. Trˆ ımbit¸a¸s
Romanian-German University Sibiu “Babe¸s-Bolyai“ University Cluj-Napoca
Faculty of economic engineering Faculty of Mathematics and Computer
in electric, energy, electronic Science
Calea Dumbravii street nr: 28-32 Mihail Kogalniceanu street nr. 1
Romania Romania
E-mail: popdaniel31@yahoo.com E-mail: tradu@math.ubbcluj.ro
129