Our main purpose is to construct one standard and three nonstandard finite difference schemes for the cube-root differential equation. After an analysis of the general qualitative features of the solutions to this equation and a calculation of the exact period, we study the stability of the numerical solutions for the four discretization schemes. Our general conclusion is that the standard forward-Euler method gives unstable numerical solutions, while the three nonstandard schemes provide suitable integration procedures

Ehrhardt, Matthias

Mickens, Ronald E.

DepositOnce

Discrete Models for the Cube–RootDifferential EquationMatthias Ehrhardt1, Ronald E. Mickens2ehrhardt@math.tu-berlin.de, rohrs@math.gatech.edu1 Institut für Mathematik, Technische Universität Berlin, Str. des 17. Juni136, D–10623 Berlin, Germany3 Department of Physics, Clark Atlanta University, Atlanta, GA 30314,USAAMS 2000 Subject Classification: 34C15, 34C25, 65L07, 65L12Key words: Nonlinear oscillations, periodic solutions, nonstandard finitedifferences, stability of numerical solutionsOctober 17, 2007AbstractOur main purpose is to construct one standard and three nonstan-dard finite difference schemes for the cube–root differential equation.After an analysis of the general qualitative features of the solutions tothis equation and a calculation of the exact period, we study the stabil-ity of the numerical solutions for the four discretization schemes. Ourgeneral conclusion is that the standard forward–Euler method givesunstable numerical solutions, while the three nonstandard schemesprovide suitable integration procedures.1 IntroductionThe cube–root differential equation (CRDE) is [1]ẍ + x13 = 0, (1.1)1where the ’dot’ notation is used to indicate the derivative, i.e., ẋ ≡ dx/dt,etc. The mathematical properties of solutions to this equation has been in-vestigated by Mickens [1, 2] and a summary of these results will be presentedin Section 2. Our purpose here is to construct finite difference schemes for theCRDE and check, using numerical simulations, if they reproduce the periodicsolutions known to exist for Eq. (1.1) [1, 2]. The discretizations formulatedconsist of a standard forward–Euler scheme and three nonstandard schemesbased on the methodology of Mickens [3, 4, 5]. Thus, our interest is notin determining the overall accuracy of these schemes, but rather the evalu-ation of which can provide numerical solutions having the same qualitativeproperties of the solutions to the CRDE.In Section 2, we provide a brief summary of the general properties ofthe solutions to Eq. (1.1). The four finite difference schemes are given inSection 3 with short discussions on their derivations. Section 4 provides asummary of our numerical experiments on the four finite difference schemes.Finally, in Section 5, we discuss our results and come to several conclusionsregarding the construction of discrete models such as given by Eq. (1.1).It should be noted that, for purposes of the implementation of the variousfinite difference schemes, care must be taken with regard to the evaluationof x1/3. Certain computer software give ambiguous or no numerical eval-uations for negative x–values, e.g. MATLAB returns using x1/3 a complexvalued root; for the needed real valued cubic root one has to use the MAT-LAB function nthroot(x,3). However, we have overcome this weakness byreplacing x1/3 by the following equivalent expressionx1/3 = [sign(x)] (|x|) 13 , (1.2)wheresign(x) =1, if x > 0;0, if x = 0;−1, if x < 0.In the remainder of the paper, it should be assumed that whenever x1/3kappears, the corresponding form given in Eq. (1.2) was used for the actualnumerical simulations.To avoid the repetition of writing certain expressions over several times,the following abbreviations are used in this paper:2CRDE : cube–root differential equationNSFD : nonstandard finite differenceNSFE : nonstandard forward–EulerODE : ordinary differential equationSFE : standard forward–Euler2 The Cubic–Root Differential EquationThe initial value problem for the CRDE is [6]ẍ + x13 = 0, x(0) = A, ẋ(0) = 0. (2.1)The more general set of initial conditions, x(0) = A, and ẋ(0) = B, is notrequired since it can be easily shown that the analysis of the special case, inEq. (2.1), gives the same relevant results as the general situation for A 6= 0and B 6= 0 [7].The CRDE can be written as a system of two coupled first–order ODEsdxdt= y,dydt= −x 13 , x(0) = A, y(0) = 0. (2.2)The variables, x and y, constituent a two dimensional phase–space [7]. Notethat the fixed–point or constant solution for this system is (x̄, ȳ) = (0, 0).Further, in the (x, y) phase–space, the trajectories, y = y(x), are solutionsto the first–order ODE [7]dydx= −x13y. (2.3)Since this equation is separable, the following first integral results from inte-grating Eq. (2.3)y22+(34)x43 =(34)A43 . (2.4)The integration constant on the right–side was evaluated by use of the initialconditions in Eq. (2.2). Since all of the terms in the first integral are non–negative for all values of x and y, Eq. (2.4) represents a closed curve in thephase–space. This results implies that all of the solutions of the CRDE areperiodic except for the fixed–point located at the origin [1, 6, 7].3The ODE for the trajectories in phase–space, Eq. (2.3), is invariant underthe following three transformationsT1 : x → −x, y → y;T2 : x → x, y → −y;T3 : x → −x, y → −y.They correspond, respectively, to reflection through the y–axis, reflectionthrough the x–axis and inversion through the origin. These symmetries canalso be used to show directly that all the trajectories in phase–space areclosed curves and hence all solutions to the CRDE are periodic [1, 7].Since Eq. (2.1) has (all) periodic solutions, it is of interest to determinethe period T of these solutions, i.e.,x(t + T ) = x(t). (2.5)The following expression gives the period [2]T (A) =√322∫ A0dx√A43 − x 43. (2.6)After a rather involved set of calculations, we obtainT (A) = (2√6)φA13 , (2.7)whereφ =∫10√w(1 + w)(1 − w) dw (2.8)has a known numerical value represented in terms of complete elliptic inte-grals of the first and second kinds [8]. Putting all this together givesT (A) = (5.86966)A13 . (2.9)From a physical point of view, the CRDE represents a nonlinear conser-vative oscillator for which the elastic force term is x1/3. This fact will beuseful for the construction of one of our NSFD schemes.In summary, the CRDE corresponds to a nonlinear conservative oscillator.Except for the fixed–point or equilibrium solution, all other solutions areperiodic with constant amplitudes.43 Discretizations of the CRDEThe first finite difference scheme for the CRDE is the standard forward–Eulerrepresentation [9, 10] for the system equations; see Eq. (2.2). For this case,we havexk+1 − xkh= yk,yk+1 − ykh= −x13k , (3.1)wheret → tk = hk, h = ∆t; k = 0, 1, 2, . . . ,x(t) → x(tk) = xk.Solving for xk+1 and yk+1 givesxk+1 = xk + hyk, yk+1 = yk − hx13k . (3.2)We call this the SFE scheme for the CRDE.The second scheme isxk+1 − xkh= yk,yk+1 − ykh= −x13k+1. (3.3)Solving these expressions for xk+1 and yk+1 gives the following expressionsxk+1 = xk + hyk, yk+1 = yk − h(xk + hyk)13 . (3.4)Note that for this scheme, the x1/3 term is evaluated at tk+1 rather than tk,as was the case for the SFE discretization. We denote the scheme, Eq. (3.4),by the notation NSFE–1.For the third scheme we use the discrete constructionxk+1 − xkh= yk+1,yk+1 − ykh= −x13k , (3.5)and these equations can be rewritten to the formsxk+1 = xk + hyk − h2x13k , yk+1 = yk − hy13k . (3.6)This is scheme NSFE–2.The fourth scheme is derived from formulations derived by Mickens [11],and Anguelov and Lubuma [12]. While the derivations appear to be com-pletely different, for the CRDE the same NSFD scheme is obtained. The5basic idea is to construct a discrete form for the first integral of the CRDE;see Eq. (2.4). Denoting this by Hk = H(xk, xk−1), we have [11]Hk =12(xk − xk−1h)2+34x23k x23k−1 = const. (3.7)Since Hk is constant, thenHk+1 = Hk. (3.8)After some algebraic work, Eq. (3.8) can be rewritten to the formxk+1 − 2xk + xk−1h2+3x23kx23k+1 + x13k+1x13k−1 + x23k−1x13k+1 + x13k−12 = 0 (3.9)Observe that the second–order derivative d2x/dt2 is represented by a (stan-dard) central difference scheme, while the x1/3 term has a very complexalgebraic structure. It is very clear that the scheme of Eq. (3.9) would notbe formulated within the context of standard numerical methods for ODEs[11, 12]. We denote the scheme of Eq. (3.9) by NSFE–3.4 Numerical ResultsNumerical simulations were carried out for the four finite difference schemesconstructed for the CRDE using the step size h = 0.01 until the final timeT = 100 (i.e. 10,000 steps). For all of the simulations, the initial values wereselected to bex(0) = 1, ẋ(0) = y(0) = 0; (4.1)consequently, x0 = 1 and y0 = 0. The NSFE–3 scheme is a second orderdifference equation and therefore both x0 and x1 are needed. The value forx1 was calculated by the following evaluation:x1 = x(h) = x(0) + h ẋ(0) +h22ẍ(0) + O(h3)= x0 + hy(0) +h22[−(x0)13]+ O(h3)= 1 − h22+ O(h3).(4.2)6The first two terms on the right–side of Eq. (4.2) were then used to calculatex1.Figure 1 presents the results for the standard forward–Euler scheme; seeEq. (3.2). Observe that this scheme produces numerical solutions that oscil-late with increasing amplitude. This behaviour is clearly inconsistent withthe known properties of the CRDE: all its solutions oscillate with constantamplitude, i.e., the general solution (except for the fixed–point) is periodic.Therefore, we conclude that the SFE scheme is unstable and should not beused to calculate numerical solutions for the CRDE.0 10 20 30 40 50 60 70 80 90 100−2−1.5−1−0.500.511.52txCube−Root problem: SFESFigure 1: Numerical solution to the CRDE using the SFE scheme.It should be indicated, for dynamic systems having either a linear or non-linear center fixed–point, the use of forward–Euler discretizations, in general,transforms the center into an unstable node for the corresponding differencescheme [13]. In addition to the paper of Wang et al. [13], further work onthis issue has been carried out by Mickens [14], Sanz–Serna [15], Mickens[16], and Roeger [17].Of both interest and importance, the three nonstandard finite differenceschemes gave essentially identical numerical results for the integration ofthe CRDE; the corresponding three schemes are given in Eqs. (3.4), (3.6),and (3.9). Observe, from Figure 2, that for all three nonstandard schemesthe amplitudes of the oscillations are constant and that computationally70 10 20 30 40 50 60 70 80 90 100−1−0.8−0.6−0.4−0.200.20.40.60.81txCube−Root problem: nonstandard schemeFigure 2: Numerical solution to the CRDE using one of the nonstandardschemes.the solutions are periodic. Clearly the NSFE schemes produce numericalsolutions having all of the essential mathematical features of the solutions ofthe CRDE.Finally, we tried for a comparison, the standard MATLAB one–step ODEsolver ode45 that is based on an explicit Runge-Kutta (4,5) formula, theDormand-Prince pair [18]. We prescribed a relative error tolerance of 10−4and computed the solution until the extended final time T = 300. In Figure 3,one clearly observes that the amplitude of the oscillations is decaying intime. Thus, this solution behaviour is another motivation to use nonstandardschemes like Eqs. (3.4), (3.6), and (3.9) to solve the CRDE numerically.5 DiscussionOur work has demonstrated that the use of nonstandard finite differenceschemes allows the construction of dynamic consistent [19] discretizationsfor the CRDE. However, it would be of value to gain insight as to why theSFE scheme did not work, i.e., why it produced an unstable solution. Togain this insight into a possible mechanism, we will analyze a much simpler80 50 100 150 200 250 300−1−0.8−0.6−0.4−0.200.20.40.60.81txCube−Root problem: ODE45Figure 3: Numerical solution to the CRDE using the MATLAB solver ode45.system, namely, the linear harmonic oscillator [4]ẍ + x = 0, (5.1)which can be written in the system formdxdt= y,dydt= −x. (5.2)The SFE scheme for these equations arexk+1 − xkh= yk,yk+1 − ykh= −xk. (5.3)Eliminating yk and shifting the index k down by one gives the followingsecond–order difference equationxk+1 − 2xk + xk−1h2+ xk−1 = 0. (5.4)To find the approximating differential equation [4] make an expansion of xk−1,xk−1 = x(tk − h)= x(tk) − hẋ(tk) + O(h2),(5.5)9substitute this expression into Eq. (5.4), and keep terms to order h. Carryingout this procedure gives the ODEẍ − hẋ + x = 0. (5.6)Note that the solutions to this modified harmonic oscillator equation oscillatewith increasing amplitude [4]. Thus, the conclusion is that the SFE schemefor the harmonic oscillator equation changes the center fixed–point, for theODE, into an unstable node for the difference scheme.We believe that the same mechanism is at work when the SFE scheme isapplied to the CRDE. This clearly is indicated by the elimination of yk inEq. (3.1) to obtainxk+1 − 2xk + xk−1h2+ x13k−1 = 0, (5.7)and this second–order nonlinear difference equation has a structure similarto Eq. (5.4).In summary, we have constructed four finite difference schemes to numer-ically integrate the CRDEẍ + x13 = 0. (5.8)All of the NSFD schemes produced numerical solutions that were dynami-cally consistent with the actual solutions of the original differential equation.However, the SFE scheme solutions were unstable and this fits in with pre-vious experiences with this scheme [13]–[17]. The conclusion is that NSFDschemes can provide valid discretizations for the CRDE. A future researchproblem is to carry out a mathematical and discretization analysis of therelated nonlinear ODEẍ +1x13= 0. (5.9)Phase space analysis and the existence of a first–integral shows that all so-lutions to Eq. (5.9) are periodic.AcknowledgementThe research and related work of the second author were supported by theUnited States Department of Energy. We both thank Kale Oyedeji (More-house College) for his assistance with various aspects of the numerical workand their display.10References[1] R.E. Mickens, Oscillations in a x1/3 potential, J. Sound and Vibration246 (2001), 375–378.[2] R.E. Mickens, Iteration method solutions for conservative and limit–cycle x1/3 force oscillators, J. Sound and Vibration 292 (2006), 964–968.[3] R.E. Mickens, Exact solutions to a finite–difference model of a nonlinearreaction–advection equation: Implications for numerical analysis, Num.Meth. Partial Differential Eqs. 5 (1989), 313–325.[4] R.E. Mickens, Nonstandard Finite Difference Models of DifferentialEquations, World Scientific, Singapore, 1994.[5] R.E. Mickens, editor, Applications of Nonstandard Finite DifferenceSchemes, World Scientific, Singapore, 2000.[6] S.L. Ross, Differential Equations, Blaisdell, London, 1964.[7] R.E. Mickens, Mathematical Methods for the Natural and EngineeringSciences, World Scientific, Singapore, 2004.[8] I.S. Gradshteyn and I.M. Ryzhik, Tables of Integrals, Series and Prod-ucts, Academic Press, New York, 1980.[9] F.B. Hildebrand, Finite–Difference Equations and Simulations,Prentice–Hall, Englewood Cliffs, NJ, 1968.[10] J.D. Lambert, Numerical Methods for Ordinary Differential Systems,Wiley, New York, 1991.[11] R.E. Mickens, Construction of finite dfference schemes for coupled non-linear oscillators derived from a discrete energy function, J. Diff. Eqs.Appl. 2 (1996), 185–193.[12] R. Anguelov and J.M.-S. Lubuma, Contributions to the mathematics ofthe nonstandard finite dfference method and applications, Num. Meth.Partial Differential Eqs. 17 (2001), 518–543.[13] X. Wang, E.K. Blum and Q. Li, Local dynamics and bifurcation con-sistencies of continuous–time dynamical systems and their numericaldiscretizations, Proceedings of Symposia in Applied Mathematics 48(1994), 399–402.[14] R.E. Mickens, Analytical and numerical study of a non–standard finitedifference scheme for the unplugged van der Pol equation, J. Sound andVibration 245 (2001), 757–761.11[15] J.M. Sanz–Serna, An unconventional symplectic integrator of W. Kahan,Appl. Numer. Math. 16 (1994), 245–250.[16] R.E. Mickens, A nonstandard finite difference scheme for the Lotka–Volterra system, Appl. Numer. Math. 45 (2003), 309–314.[17] L.-I.W. Roeger, Nonstandard finite difference schemes for the Lotka–Volterra systems: Generalization of Mickens’ method, J. Diff. Eqs. Appl.12 (2006), 937–948.[18] J.R. Dormand and P.J. Prince, A family of embedded Runge-Kutta for-mulae, J. Comp. Appl. Math. 6 (1980), 19–26.[19] R.E. Mickens, Dynamic consistency: A fundamental principle for con-structing nonstandard finite differences for differential equations, J. Diff.Eqs. Appl. 11 (2005), 645–653.12

Discrete Models for the Cube-Root Differential Equation

https://depositonce.tu-berlin.de/bitstream/11303/16878/1/Report-40-2007.pdf

Discrete Models for the Cube-Root Differential Equation

Abstract

Similar works

Full text

Available Versions

DepositOnce