CC BY
2URAL MATHEMATICAL JOURNAL, Vol. 10, No. 1, 2024, pp. 18-27
DOI: 10.15826/umj.2024.1.002
AN EXPLICIT ESTIMATE FOR APPROXIMATE SOLUTIONS OF ODES BASED ON THE TAYLOR FORMULA
Abdulla Kh. Abdullayev
V.I. Romanovskiy Institute of Mathematics of Uzbekistan Academy of Sciences, 9 University Str., 100174 Tashkent, Uzbekistan
abduganiax@mail.ru
Abdulla A. Azamov
V.I. Romanovskiy Institute of Mathematics of Uzbekistan Academy of Sciences, 9 University Str., 100174 Tashkent, Uzbekistan
abdulla.azamov@gmail.com
Marks B. Ruziboev
Kimyo International University, 156 Shota Rustaveli Str., 100121 Tashkent, Uzbekistan
marx.ruziboev@gmail.com
Abstract: In this paper, we consider a third-order explicit scheme based on Taylor's formula to obtain an approximate solution for the Cauchy problem of systems of ODEs. We prove an estimate for the accuracy of the approximate solution with an explicit constant that depends only on the right-hand side of the equation and the domain of the solution.
Keywords: Dynamical systems, Cauchy problem, Approximate solution, Taylor formula, Accuracy of approximate solution, Level of accuracy, Error term.
1. Introduction
It is needless to note the importance of estimating accuracy for approximate solutions of ODEs. Here we consider the problem
X(t) = f (x(t)), x(0) = X0, (1.1)
where x € D c Rd and D is a convex domain. In what follows, we assume that the function f : D ^ Rd is three times differentiable with continuous derivatives in D. In practice, establishing the highest possible accuracy of an approximate solution is one of the key problems. Thus, the efficiency of an approximate solution is determined by its accuracy. Let x*(t) be a solution of (1.1) in the interval 0 < t < T for some T > 0, and let X(t) be its approximate solution (obtained by some scheme) on the same interval. The accuracy of the scheme is expressed by an inequality of the form
sup |x*(t) - X(t)|< CeLThs, (1.2)
0<t<T
where L is the Lipschitz constant of the function f, h is the mesh size, and s is the order of accuracy of the method. Approximate solution schemes can be implicit and explicit. In this paper, we are
concerned with explicit methods. Frequently used one-step approximate solution schemes can be divided into two groups:
(1) Schemes based on Taylor's expansion of the solution;
(2) Runge-Kutta-type methods.
Schemes based on Taylor's expansion are easy to implement, but experts prefer Runge-Kutta-type methods. This preference is caused by the fact that the error estimates for schemes based on Taylor's formula contain derivatives of the function f, which can be challenging to estimate. However, with the development of computer algebra, computations of derivatives of a rather wide range of functions can be automated [7, 12, 13]. Therefore, schemes based on Taylor's formula can be implemented without any extra hurdles, due to the simplicity of implementation.
On the other hand, since the 1960s, computer simulations have been used extensively to study dynamical systems described by nonlinear systems of ordinary differential equations (ODEs). Recently, with the rise of computational power and computers being widely available, computerassisted proofs have come into play. These proofs, however, require verification of the accuracy of the scheme, i.e., proof of inequality (1.2) with explicit constants.
Estimates for the accuracy of approximate solutions to ODEs are studied extensively in literature [2, 4, 6]. For example, the monographs [10, 14, 16] contain estimates for Runge-Kutta methods. Part III of the well-known monograph [3] is devoted to approximate solutions of ODEs. Schemes based on Taylor's formula are briefly discussed in the first section from a methodological viewpoint. But no estimates are provided. It is surprising that, among the vast amount of literature devoted to the approximate solutions of ODEs, we did not find estimates with explicit constants. Most estimates give the order of approximation, which is insufficient if we want to use approximate solutions in the proofs. Except [10], which, citing [5], gives an inequality for the Runge-Kutta method. In [15, Sect. 2, Part II], the authors mention a scheme based on Taylor's formula, but do not consider the problem of accuracy estimates. Some authors claim that if one takes the first n terms of Taylor's expansion, then the error term will be of the form
But they neither provide proof nor speak about constants involved in O(hn).
In [7], which is one of the most comprehensive monographs on approximate solutions of ODEs, the authors claim that the difference between exact and approximate solutions is estimated by the remainder term of Taylor's expansion and in just one step it will be O(hp+1) [8, Sect. 318, p. 180]; no further details are given. In [17, 18], the problem of estimating the error is investigated for methods of approximation of the integral
Paper [8] provides explicit estimates for the approximate computation of this integral. This is equivalent to considering the first term of Taylor's expansion, which provides the first-order Taylor approximation scheme. The author considers multi-step approximate solution schemes but does not give inequalities of the form (1.2).
It turns out that schemes based on Taylor's formula are more convenient than schemes based on Runge-Kutta methods for obtaining explicit constants in (1.2). For example, in [1], an estimate of type (1.2) is obtained for a second-order scheme based on Taylor's formula. Also, they are easier to implement in numerical approximations than Runge-Kutta-type methods. Keeping in mind an application of the approximate solution schemes in computer-aided proofs, in the present paper, we consider a third-order scheme based on Taylor's formula for Cauchy problem (1.1) and give an explicit constant C, for which inequality (1.2) holds with s = 3.
max |x
0 <n<N
n
y(xn) | = O(hn).
t+h
t
2. Approximate solution schemes based on Taylor's formula
For a continuously differentiable function f (x), the initial value problem (1.1) has a unique solution; however, it is challenging to estimate the interval where the solution exists [11]. On the other hand, for approximation schemes, we need to know the existence of the solution. Therefore, our first standing assumption is the following.
Assumption A. Fix T > 0. The solution x*(t) to the Cauchy problem (1.1) is defined on the interval [0,T].
Usually, when considering approximate solutions, one tries to find values of the approximate solution x*(t) on a mesh 0 = t0 < ti < t2 < ■ ■ ■ < tn = T. Here, to simplify the exposition, without loss of generality, we consider the uniform mesh tn = nh, where h = T/N, N € N+.
We start with the Taylor expansion of the exact solution with accuracy O(h3):
h,2 h3
x*(t + nh) = x*(t) + hx*(t) + —x*(t) + —x*(t) + E4(i, h).
26
By definition, we have
X *(t) = f [x*(t)j, X*(t) = f '[x*(t)]f [x*(t)j, £.*(t) = f"Mt)]f [x*(t)]f [x*(t)j + f '[x*(t)]f '[x*(t)]f [x*(t)j.
We also need the fourth derivative xIV(t). To shorten the notation, we interpret derivatives of the function f as operators acting on x*(t) and write
x *(t) = f [x*(t)],
x*(t) = (f 'f )[x*(t)],
x.*(t) = (f ''ff + f'f'f )[x*(t)]. In particular, on the uniform mesh,
h2
+ 1 )h] = x*(nh) + hf[x*(nh)] + — (f' f)[x*(nh)\
+ Ku"ff + /77)k*(«/^)] + n = 0,1,2,..., TV - 1.
6
By neglecting the remainder term, we obtain the following recurrent formula for the approximate solution:
= hf(xn) + y (/7)(-^) + ]-jU"ff + f'f'f)(xn). (2.1)
It is expected that elements of the sequence xn (n = 1,2,...,N) defined by (2.1) are close to the values of the exact solution at the points x*(h),x*(2h),... ,x*(Nh). An intuitive way to measure this closeness is to compute the value
p = max |x*(nh) — xn|. (2.2)
1<n<N
However, the quantity (2.2) does not provide any information about the behavior of the solution on the interval ((n — 1)h, nh). The main aim of this paper is to derive such estimates. If we want to apply numerical solutions in computer-aided proofs, then we cannot ignore the behavior of the
system on the interval (nh, (n + 1)h). Certainly, for t e (nh, (n + 1)h), we can estimate x*(t) as follows. For t e [0,T], let n(t) = [t/h]. Then |t - n(t)h| < h and
h
|x*(t) - x*(n(t)h)| < J |f [x(s)]|ds < Moh,
n(t)h
where M0 is the maximum value of |f (x)| on some compact subset of D. The above inequality together with (2.2) imply that
|x*(t) - xn| < Moh + |x* [n(t)h] - xn| < Moh + p,
which shows that the difference between the exact and approximate solutions is, at best, of order h. However, this estimate is very rough and insufficient for our purposes. To get a better estimate, we use generalized Euler polygons.
From now on, we fix N, h, and T > 0 such that Nh = T and define
at = nh if t e [nh, (n + 1)h) for n = 0, ...,N - 1.
We start with the definition of an approximate solution.
Definition 1. A continuous function x(t) : R ^ [0 : T] satisfying the equation t
m =X0 + J (f[x(<Ts)] + (s- as)(f'f)[x(as)} + ~°s) (f'ff + f f f)[x(a8)])ds (2.3) o
is called an approximate solution of (1.1).
Although (2.3) looks like an integral equation, it is a recurrent formula, and we can construct a?(t) explicitly step by step. Therefore, (2.3) defines a function x(t) on [0,T] as an approximate solution; i.e., our notion of approximate solution is well defined.
Lemma 1. The equality i(nh) = xn holds for all n > 0.
The lemma is proved easily by induction on the intervals [0,nh]. We use this lemma to reduce the problem to deriving estimates of the form (1.2) for the difference |x*(i) — i(t)|.
The remainder term of the Taylor expansion is given by the formula (see for example, [9])
1
R4(t, h) = J (1~S) xiIV)(t + sh)h4ds, 0
where
xiIV )(s) = (f '''fff + 3f f 'ff + f 'f ''ff + f 'f 'f 7) [x*(i)].
By the formula for the error term is estimated by the maximum values of the derivatives of f. To estimate the derivatives effectively along the solution and approximate solution, we need them to stay in some compact set. Therefore, we assume the following.
Assumption B. Let K be a convex and compact domain in Rd. We assume that values of the exact solution x*(i) and the approximate solution ic(i) remain in K for all t € [0,T].
Since, K is compact by definition, Assumptions A and B allow us to define
Mo = max |f (x)|, Mi = max ||f'(x)||,
xeK xeK (2 4)
M2 = max ||f''(x)||, M3 = max ||f'''(x)||. (.)
x€K x€K
Using these quantities, we obtain the following estimate for the remainder term of the Taylor formula:
h4 h4
\Ri(nh, h)| < — max^ |4V (s)|ds < — (M3Mq + 4M2MiM^ + A/,*" A/..) .
In the literature, when estimating the accuracy of approximate solutions, many authors claim that the error is bounded by |R4(nh, h)|. However, this is not true, since R4(nh, h) is the difference between x*(t) and its Taylor expansion, which does not directly imply any conclusions for the difference x*(t) — x(t). In the present paper, we show that it is possible to obtain an explicit estimate for the latter. The main result of the this paper is the following theorem.
Theorem 1. Under Assumptions A and B, the following inequality holds:
e 1 - > - > 2\L3
Mi) - x(t)I < ——-(Lo + L\h + L2hz)h
oMi
where
L0 = 5M02M1M2 + M0M3 + M03M3, Li = \(MoM2 + 4Mq M1M3 + 9MQM'IM2), l2 = \(M*M2M3 + MqM^MS + 2MQM\M2 + 2M$MfM2).
Note that our estimate for the accuracy of the method is explicit and can be computed effectively in terms of the right-hand side of the initial value problem (1.1).
3. Proof of Theorem 1
In this section, we prove the main theorem. The key ingredient of the proof is a discretization of the time t using a piecewise constant function
In what follows, we repeatedly use the following formula for the derivative of the approximate solution. By (2.3), for t = nh (n = 1,2,..., N — 1),
m = f[x(<rt)] + (t- at)(f'f)[x(at)\ + i^lV/2 + (f')2f)[x(at)]. (3.1)
Further, taking into account that 0 < t — at < h, we have
I Hat) I < M0 (l + hMi + y (M2M0 + Mf)). (3.2)
Proof. For the exact solution, we have the equality
t
x*(t) = xo + | f [x*(s)]ds.
Using this equality and equation (2.3) and adding and subtracting the term
/ f[x(s)]ds, Jo
we obtain
t t
|x* (t) - £(t)| < JI (s)ds + J \f [x*(s)] - f [£(s)]|ds,
0 0
where
I(s) = f[x(s)\ - f[x(as)\ -(s- as)(f'f)[x(as)\ - {$ ~°s)\f"f2 + f'2f)[x(as)}.
We are going to derive an upper estimate of the form Ch3 for I(s). By the fundamental rule of the calculus, we have
s s
f[x(s)} - f[x(as)} = f = j f'[x(r)]B(r)dr.
o.
Substituting the derivative of the approximate solution x(r) given in (3.2) into the right-hand side, we obtain
s
f[x(s)]-f[x(as)}= J f'[^r)}{f[^ar)}Hr-ar)(f'fmar)} + (r U" f+f'2 f)[x^r)] }dr.
0s
Note that ar = as in this equation since r € [as,s] and a is piecewise constant. Denote the latter term on the right-hand side by C1:
s
Ci= f (r -as)2 (/„/2 + f2fmar)]mr)]dr_
o,
Using (2.4), we obtain the estimate
|Ci| < (M2MiMq + MfMo) [ ^ ^ dr < Dxh3, (3.3)
where D1 = (M2M1M02 + M3M0)/6. Therefore, we obtain
s
I(s) = / if'[^(r)] — f'[^(as)]}f [®(as))]dr
0s
s
+1 f'[x(r)](ff)[x(as)}dr - {S ~2as)2 (f'f + fr2f)[x(as)} + Ci-Denote the first term of this expression by J(s):
s
J (s)= / f' [X(r)] — f' [£(as)]}f [x(as)]dr.
(3.4)
o
s
o
s
o
m
We estimate J(s), using the fundamental rule of calculus and the derivative of the approximate solution given by (3.2):
Define
s r s r
J(s) = J J ^^duf[x(as)]dr = 11 f"[x(u)]{f[x(au)}
as CTs CTs CTs
+ (u - au)(ff)[x(au)\ + {U~2au)\f"f + (f')2f)[x(au)]}f[x(a8)]dudr.
s r
C = J J {(u — au)f//[J(u)](f/f)[£K)]f [a^)]
as as
+ {ll~2au)2 f"[x(uW"f2 + (f')2f)[x(<ru)]f[x(<Ts)]}dudr. Thus, taking into account that au = as and using (2.4), we obtain the inequality
(3.5)
(3.6)
s r
|c2|<//
M2M\Mq(u - a8) + ( Л/7Л/.? + М2М?М^и
dudr < + D3h4, (3.7)
where
D2 = iMiM2M02, £>3 = T^i.UjU,'' + M2M2M2). o 24
Consequently, substituting (3.6) into (3.5) and then (3.5) into (3.4) and denoting C1 + C2 by C3, we obtain
s r s
I (s)=/ J f"[ir(u)]f [£(u)]f [i(as)]dudr + J f [x(r)](r — ar )(ff )[i(ar )]dr
as
(f'f2 + f '2f )[^(^s)] + C3 = A(s) + B (s) + C3,
(3.8)
(s — rll ¿2 , f/2
2
where we used the following notation:
s r
A(s) = J j {(f"f)[i(u)] — (ff)[®(as)]}f [i(as)]dudr,
as as s
B (s) =J(r — as){f [®(r)] — f [ir(as)]}(ff )[i(^s)]dr.
as
Combining (3.3) and (3.7), we obtain
|C3| < D1h3 + D2h3 + D3h4. (3.9)
It remains to estimate A(s) and B(s). We have
s r u
a a
s
which implies
s r u
A(s) = //(/ f^M^f[X(v)] + f"M)]f [x(v)]£(v)}dv)f[£(*s)]dudr. (3.10)
Vs Vs Vs
Therefore, taking the absolute value of the expression under the outer integral (integration with respect to r), using estimates (2.4) and (3.2), and taking into account that t — at < h and au = as = ar, we obtain
s r u ^
|A(s)| < j j j {Mo + hMiMo + yi.\/2.\/,f + .\/f.\/„))(M3Mq + M0M1M2)dvdudr
J s Vs Vs
h2 , h3
< i.\/,.\/,f + M0MiM2) ( M0 + + — (M2MQ + MfM0)J — .
Similarly,
r
B(s) = J(r-<r8){J -^f'[x(u)]duyf'f)[x(a8)]dr
s
Vs Vs
s r
= I (r - ffsW f//[£(u)];?(u)du}(f/f)[x(CTs)]dr.
s Vs
Again, using (2.4) and (3.2), and taking into account that t — at < h and as = ar in the above equation, we obtain
(r_n" )3 / h2
|5(s)| < --—^-MoMiM2(Mo + // A /1 A A. + — (M2Mq + .\/f.\/„)
3 2
u2 k 3
< M0MiM2 (Mo + /; /1 /,, + — i.\/2.\/,T + .\/f.\/,,)N 1
(3.11)
2 v ' 0 1 3
Finally, substituting (3.9), (3.10), and (3.11) into (3.8), we obtain the inequality
|/(s)| < i(M2MiM02 + .\/i;.\/,,))/; '' + ^MiM2Mq/?3 + —i-( \/7 \/.? + M2M2Mq )/?4 6 6 24
/ h2 \ h3
+(M3M2 + M0M1M2) f Mo + + — (M2Mq + M2M0)J —
/ h2 \ h3
+M0MiM2 f M0 + + — (M2Mq + M:2M0)J y
= ^(L0 + £i/?, + L2/?2)/?3, 6
with the required constants L0, L1, and L2.
Now, we use the compactness of the domain K and smoothness of f to obtain
|f [x*(s)] — f [x(as)]| < M1|x,(s) — x(s)|. (3.12)
Using inequalities (3.1), (3.2), and (3.12), we get the inequality
t t
|x*(t) _ x(t)| < J |I(s)|ds + J |f Ms)] _ f [x(as)]|ds oo
t
<^{L0 + L\h + L2h2)hH + J Mi Ms) - £(«)!•
(3.13)
s
V
s
In (3.13), we apply Gronwall's inequality. For our purposes, the following version is the most convenient. Let u : R ^ R be a continuous function such that u(t) > 0 for t > 0 and
t
u(t) < Ct + M J u(s)ds
0
for some C, M > 0. Then, the following inequality holds:
Mt — i
Applying Grownall's inequality to estimate (3.13), we obtain
1 eMit — 1
Ix*{t) - £(i)| < -(L0 + L\h + L2/?2)/?3 —--.
6 M1
This completes the proof. □
4. Conclusion
1. If a nonautonomous system is considered in a d-dimensional space, then we can interpret it as an autonomous system in the (d + 1)-dimensional space. In particular, we can consider the following Cauchy problem:
i = f x(0)=xo,
x = f (t, x), x(0) = xo d
If-1- «°> = °-
Therefore, the difference between the exact and approximate solutions can be estimated by the same expression with the constant sjMq + 1 insted of Mo.
2. In the proof of the estimate for the difference between the exact and approximate solutions, we obtained
eMi t — 1
max \xJt) - x(t)\ < C——-h3
o<t<T ' Wl " M1
with the coefficient that is considerably larger than expected, where the constant C is a fifth-order polynomial of the constants M0, M1, M2, and M3. If these constants are not very large, i.e., on the order of 1, then the coefficient C does not affect the choice of the mesh size. In this case, the mesh size would mostly depend on
(eMlt — 1)/M1. On the other hand, if the constants M0, M1, M2, and M3 are on about 10, then the coefficient at h3 is on the order of 105, which would affect the choice of h essentially, in certain cases, it may even invalidate the approximate solution scheme.
3. Another interesting question is whether it is possible to simplify the proof of the main theorem. The authors think that the proofs cannot be simplified considerably.
4. It is possible to prove a similar theorem with the order of accuracy h4; i.e., we can consider the fourth-order scheme
h2 h3
= Zn + hf(xn) + -U'f){xn) + - (/7/ + /77 )(Xn) +
(/'"/// + 3/77/ + /77/ + /777) M
and prove an analogous theorem.
REFERENCES
1. Azamov A. A., Abdullaev A. X., Tilavov A.M. On the derivation of an inequality for estimating the accuracy of an approximate solution to the initial value problem. Bull. Inst. Math2022. Vol. 5, No. 5. P. 105-111.
2. Babenko K.I. Foundations of Numerical Analysis [Osnovy chislennogo analiza]. 2nd edition. Moscow, Izhevck: Regulyarnaya i haoticheskaya dinamika, 2002. 848 p. (in Russian)
3. Bakhvalov N. S. Numerical Methods (Analysis, Algebra, Ordinary Differential Equations) [CHislennye metody (analiz, algebra, obyknovennye differencial'nye uravneniya)]. Moscow: Nauka, 1975. 632 p.
4. Berezin I.S., Zhidkov P.N. Computing Methods. Vol. 2. Pergamon Press, 1965. 464 p.
5. Bieberbach L. On the remainder of the Runge-Kutta formula in the theory of ordinary differential equations. J. Appl. Math. Phys. (ZAMP), 1951. Vol. 2. P. 233-248. DOI: 10.1007/BF02579687 (in German)
6. Brenan K. E., Campbell S. L., Petzold L. R. Numerical Solution of InAtia,l-Va,lue Problems in Differential-Algebraic Equations. Philadelphia: SIAM, 1995. XII+251 p. DOI: 10.1137/1.9781611971224
7. Buchberger B., Collins G.E., Loos R., Albrecht R. Computer Algebra: Symbolic and Algebraic Computation. Vienna: Springer, 1983. VII+283 p. DOI: 10.1007/978-3-7091-7551-4
8. Butcher J.C. Numerical Methods for Ordinary Differential Equations. 3rd edition. United Kingdom: John Wiley & Sons Ltd., 2016. 544 p.
9. Cartan H. Formes Differentielles. Hermann Paris, 1967. 190 p. (in French)
10. Hairer E., Norsett S., Wanner G. Solving Ordinary Differential Equations I. Nonstiff Problems. 2nd edition. Berlin, Heidelberg: Springer, 2008. XV+528 p. DOI: 10.1007/978-3-540-78862-1
11. Hartman P. Ordinary Differential Equation. SIAM, 2002. XVIII+624 p. DOI: 10.1137/1.9780898719222
12. Gathen J., Gerhard J. Modern Computer Algebra. 3rd edition. Cambridge University Press, 2013. 796 p. DOI: 10.1017/CBO9781139856065
13. Geddes K.O., Czapor S.R., Labahn G. Algorithms for Computer Algebra. NY: Springer, 2007. XXII+586 p. DOI: 10.1007/b102438
14. Kendall A. E. An Introduction to Numerical Analysis. 2nd edition. John Wiley & Sons, Inc., 1989. 693 p.
15. Krilov V. I., Bobkov V. V., Monastirniy P.I. Computational Methods [Vychislitel'nye metody]. Vol. II. Moscow: Nauka, 1977. 399 p. (in Russian)
16. Lambert J. D. Numerical Methods for Ordinary Differential Systems: The Initial Value Problem. Wiley, 1992. 304 p.
17. Milne W. E. Numerical Calculus. Princeton University Press, 1949. 404 p.
18. Milne W. E. The Remainder in linear methods of approximation. J. Res. Nat. Bur. Stand., 1949. Vol. 43. Art. no. RP2401. P. 501-511.
