When constructing an algorithm for the numerical integration of a differential equation, one should first convert the known ordinary differential equation (ODE) into an ordinary difference equation. Given this difference equation, one can develop an appropriate numerical algorithm. This technical note describes the derivation of two such ordinary difference equations applicable to a first order ODE. The implicit ordinary difference equation has the same asymptotic expansion as the ODE itself, whereas the explicit ordinary difference equation has an asymptotic that is similar in structure but different in value when compared with that of the ODE

Freed, Alan D.

Walker, Kevin P.

English

NASA Technical Reports Server

/i!_i-',_
NASA Technical Memorandum 104530
"L
//
,' ¢"
/< /
From Differential to Difference
Equations for First Order ODEs
Alan D. Freed
Lewis Research Center
Cleveland, Ohio
and
Kevin P. Walker
Engineering Science Software, Inc.
Smithfield, Rhode Island
k
July 1991
(,,j/'"A-[','_-_L_,',.-_':'L_} rr_,M '_'I::;--CRrNTiAL 1_,,
i,"!A':;A ) _ _, C(;Ct
7::
,'-s/o-_
% :_j.-._ 7 : " :.
https://ntrs.nasa.gov/search.jsp?R=19910018572 2020-03-19T16:45:22+00:00Z

From Differential to Difference Equations
for First Order ODEs*
Alan D. Freed
National Aeronautics and Space Administration
Lewis Research Center
Cleveland, Ohio 44135
Kevin P. Walker
Engineering Science Software, Inc.
Smithfield, Rhode Island 02917
Abstract
When constructing an algorithm for the numerical integration of a differential equation,
one should first convert the known ordinary differential equation (ODE) into an ordinary
difference equation (OAE). Given this difference equation, one can then develop an
appropriate numerical algorithm. This technical note describes the derivation of two
such OAEs applicable to a first order ODE. The implicit OAE has the same asymptotic
expansion as the ODE itself; whereas, the explicit OAE has an asymptotic expansion
that is similar in structure but different in value when compared with that of the ODE.
Many physical processes can be represented by systems of first order ODEs of the form
or equivalently as
3ro
JC_ + P_X_ = Q_ (for a = 1,2,...,N), (1)
=
- Q_[X_,t]-Pa[X_,t]Xa (fora,/3= 1,2,...,N), (2)
where the X_ are the N independent variables to be solved for, which may be scalar, vector
or tensor valued in applications. The parameters P_ (a scalar), and F_ and Qo (of the
same rank/type as X_), are functions of the variables XZ and t, in general. If neither
P_ or Q_ depends on X_ for all a, fl, then the system of equations is said to be linear;
otherwise, it is nonlinear. The dot '" is used to denote differentiation with respect to time,
t. We choose time to be the dependent variable for illustrative purposes, as it so often is in
physical applications; however, this is not a necessary restriction on the theory presented
herein.
*A technical note prepared for the Journal o.[ Computational Physics.
The ODE of equations (1) and (2) is a point function in time. Numerical integration
algorithms, however, are based on functions evaluated over an interval in time. We therefore
introduce the first order OAE
AX_
At - J:o[X_,At,r], (3)
where
and therefore
AX. = X.[t+At]- Xo[t], (4)
obtain the recursive integral equation [1]
X_,[t+At] = exp a_=t
+ J(=t exp -
P_ [X_, (] d(] X_ [t]
e [xe,(ld( a(,
J(=(
(6)
which is an exact solution to this first order ODE. John Bernoulli [2] developed a nonre-
cursive solution similar to (6) in 1697 for the equation )( = aX + bX n, whose solution
was sought by Jacob Bernoulli [3] in 1695. John's solution was expressed as a quadrature
because the integral of dz/z in the form of a logarithm was not generally known until later
that same year [4].
The authors have obtained a variety of approximate and exact solutions to (6) by rep-
resenting the parameters P_ and Q_ with series expansions, and integrating them term
by term. The resulting linear approximation, which was acquired by using Taylor series
expansions, is given by [1]
This approximation becomes exact whenever P_ and Q_ are both constants, i.e. whenever
the ODE is linear.
Rearranging the recursive solution (7) into a difference equation, one obtains two new
relationships; they are, for the implicit case,
(1--e-P_[X_'t+At]'At_ (P_[Xz,t+Atl.AX_ F_[X_,t+:_[Xz,At, t+At] = -p--[_Z,t+At_.At .] + At])
o[ o q [x ,t+At]+ O(t_At) _',_ ' (8)
x.[t + At] = x.[t] + Tl./,,t. (5)
The OAE represents the ODE in numerical integration algorithms. The formulation is
implicit when v = t + At; whereas, it is explicit when r = t. The development of existing
algorithms begins with the tacit assumption that ._'_ = F_. We shall now show that this
assumption is in error. This is accomplished by analytically integrating the ODE of (1) in
a recursive manner, which permits the difference function 9:"a in (3) to be determined for
the OAE.
One can introduce an integrating factor into the differential equation (1) and thereby
andfor the explicit case,
&Ex ,At, tl= (1_ e-PotX.J , 
-fi__ / F_[XI_'t]
+0 o q.[xz, t]t] (9)
Notice the presence of an extra AXa term in the implicit relation (8), which acts as a
correction to the derivative F_, and which is not present in the explicit relation (9).
The coefficient (1 - e-P_[X_'@/'t)/P,_[X_, r].At is a correction factor of the first order
that results from taking the differential function F_ at time v and converting it into a
difference function ._'_ at time v taken over the interval (t, t+At). (Higher order solutions
for the difference function .T'_ will be given in a future paper.) This coefficient goes to 1 in
the limit as At goes to 0, as it must so as to recover the differential; in other words,
lim .%'_[Xz, At, r] = F_[Xz, T]
At ---*0
(lo)
for both the implicit and explicit cases. A note of caution when writing computer code. In
the neighborhood of P_[Xz, r].At _ 0, one needs to expand (1- e-P°[X_'_]'At)/P_,[X_, "*'].At
into a power series to secure a sound computational algorithm.
The presence of the coefficient (1- e-P_[Xt_'r]'At)/Po[Xz, r].At also introduces desirable
asymptotic characteristics into our numerical approximations. In particular, for the implicit
case where r = t + At, the asymptotic expansion for the OAE (3 with 8) is given by
Q_[X_,t+At]
lim X.[t+At]× (11)
_t_l_g¢ P_ [X_, t + At] '
P_>O
which is also the asymptotic expansion of the ODE (1 & 2). In contrast, for the explicit
case where r = t, the asymptotic expansion for the OAE (3 with 9) is given by
lim X_[t+At] _ Q_[Xe't]
P_>O
(12)
These two asymptotic expansions differ only in when their parameters P_ and Qa are
evaluated. For large time steps, the implicit case is asymptotically accurate and stable for
exponentially decaying solutions, i.e. when P_[X_, t] > 0 V t. Stability becomes an issue
only when P_ < 0. In constrast, the explicit case will oscillate around the true solution for
large time steps, but with much less potential of becoming unbounded when compared with
equivalent algorithms constructed without our correction coefficient. These oscillations can
be mitigated only by choosing smaller time steps.
The second integral in (6) is a Laplace [5] integral where the integrand has its largest
value at the upper limit, t+At, and possesses an evanescent memory of the forcing function,
Q_[Xz,{], provided that P_[Xz, _] > 0 over the interval (t, t+At). This fading memory
means that the solution will depend mainly on the recent values of the forcing function,
and that by concentrating the accuracy on the recent past we obtain accurate asymptotic
representations of the solution.
In the implicit solution, the integrands in (6) were expanded in Taylor series about their
upper limits [1], where each integrand has its largest value and contributes the most to the
3
integral. By retainingbut a singleterm in theTaylor seriesexpansionsthe integrandsare
accuratelyapproximatedwheretheyarelargest,andthe neglectof the higherorderterms
is only felt nearthe lowerlimits whereeachintegrandcontributesonly a smallamountto
the integralbecauseof its exponentialdecayfromtheupperlimit. Theneglectof thehigher
ordertermsin the Taylorseriesthusresultsin analgorithm that is asymptoticallycorrect
at the upper limit. Normally,whentreatingasymptoticexpansions,the exponentialdecay
of the integrandallowsthelowerlimit to bereplacedwith zeroor minusinfinity to easethe
integration.This wasnot donein the presentcase,however,sothat by retainingthelower
limit as t, we obtain a uniformly valid asymptotic algorithm in the implicit approximation
provided that P_[X_, (] > 0 over the interval (t, t+At).
In the explicit solution, the integrands in (6) were expanded in Taylor series about their
lower limits [1], where the neglect of the higher order terms in the Taylor series results in
integrands that become progressively more inaccurate as they approach their upper limits
where the contribution from each integrand is most important. The explicit approximation
is not, therefore, a valid asymptotic representation of the integral when the Taylor series
is truncated at a finite number of terms. However, when P_[X_, _] < 0 over the interval
(t, t+At), the reverse situation occurs. In this case the asymptotic solution is now obtained
by expanding the integrands about their lower limits, where this region of each integrand
now contributes the most to the integral. Here the implicit method, obtained by expanding
about the upper limit, does not give a valid asymptotic representation of the integral.
In conclusion, one only needs to replace the differential function F_[X/3, r] with the ap-
propriate difference function :T'_[X/_, At, r] in many existing numerical integration methods
(e.g. Euler and Runge-Kutta) to construct an appropriate OAE for the numerical integra-
tion of a given ODE, and thereby obtain substantial improvements in their performance.
We have demonstrated this in references [1, 6].
References
[1] K.P. WALKER and A.D. FREED, "Asymptotic Integration Algorithms for Nonhomoge-
neous, Nonlinear, First Order, Ordinary Differential Equations," (NASA TM-103793,
1991). Submitted to: J. Comp. Phys.
[2] JOHN BERNOULLI, "DE CONOIDIBUS ET SPH/EROIDIBUS, quedam. Solutio analyt-
ica 2Equationis," Acta Eruditorum publicata Lipsice, pp. 113-118 (1697).
[3] JACOB BERNOULLI, "EXPLICATIONES, ANNOTATIONES ET ADDITIONES, Ad
ea, qumin Actis fup. anni de Curva Elaftica, Ifochrona Paracentrica, O Velaria, binc
inde memorata, O partita controverfa leguntur; ubi de Linea mediarum directionum,
aliisque novis," Acta Eruditorum publicata Lipsiae, pp. 179-193 (1695).
[4] E.L. INCE, Ordinary Differential Equations, (Dover, New York, 1956) p. 531.
[5] P.S. LAPLACE, Thdorie Analytique des Probabilitds, Vol. 1, (Paris, 1820).
[6] A.D. FREED and K.P. WALKER, "Exponential Integration Algorithms Applied to Vis-
coplasticity," (NASA TM-104461, 1991). To appear in: Cornplas III, edited by E. On£te
et al., Barcelona, 6-10 April, 1992.
rU/ A
National Aeronautics and
Space Administration
1. Report No.
NASA TM- 104530
4. Title and Subtitle
Report Documentation Page
2. Government Accession No.
3. Recipient's Catalog No.
5. Report Date
From Differential to Difference Equations for First Order ODEs
7. Author(s)
Alan D. Freed and Kevin P. Walker
9. Performing Organization Name and Address
National Aeronautics and Space Administration
Lewis Research Center
Cleveland, Ohio 44135- 3191
12. Sponsoring Agency Name and Address
National Aeronautics and Space Administration
Washington, D.C. 20546-0001
July 1991
6. Performing Organization Code
8. Performing Organization Reporl No.
E - 6406
10. Work Unit No.
553-13-00
11. Contract or Grant No.
13. Type of Report and Period Covered
Technical Memorandum
14. Sponsoring Agency Code
15. Supplementary Notes
Alan D. Freed, NASA Lewis Research Center; Kevin P. Walker, Engineering Science Software, Inc., Smithfield,
Rhode Island 02917. Responsible person, Alan D. Freed, (216) 433- 3262.
16. Abstract
When constructing an algorithm for the numerical integration of a differential equation, one should first convert the
known ordinary differential equation (ODE) into an ordinary difference equation (OzkE). Given this difference equa-
tion, one can the develop an appropriate numerical algorithm. This technical note describes the derivation of two such
Oz_s applicable to a first order ODE. The implicit OAF. has the same asymptotic expansion as the ODE itself;
whereas, the explicit OZkE has an asymptotic expansion that is similar in structure but different in value when compared
with that of the ODE.
17. Key Words (Suggested by Author(s))
Numerical methods
Integration
18. Distribution Statement
Unclassified - Unlimited
Subject Category 64
19. Security Classif. (of the report) 20. Security Classif. (of this page) 21. No. of pages 22. Price*
Unclassified Unclassified 6 A02
NASAFORM1626OCT86 *Forsale by theNationalTechnicalInformation Service, Springfield,Virginia 22161
PRECEDING PA_E BLA;,_K NOT FILMED
_llliBIW(lllll m',,---



From differential to difference equations for first order ODEs

Abstract

Similar works

Full text

Available Versions

NASA Technical Reports Server