# Euler method (and others) for unbounded intervals

In the course of my research I'm confronted with performing a numerical approximation of the solution of an initial value problem $$\begin{cases} y'=f(y,t),\\ y(t_{0})=y_{0} \end{cases}\quad\quad(1)$$ with $$y:I\rightarrow \mathbb{R}^d,\quad I\subseteq\mathbb{R}$$ and where $$f$$ satisfies standard regularity assumptions that guarantee unique solutions. One particular twist is that $$I$$ in my setup is known to be infinite, e.g., $$I=[0,\infty)$$ or $$I=\mathbb{R}$$.

I am familiar with standard approximation schemes, such as the Euler method to name a simple one. The issue is that these are only valid for bounded intervals $$J$$ of $$I$$, such as $$J=[-T,T]$$ for some $$T>0$$. Searching further on the internet, all texts and books on numerics of ODEs that I have consulted so far also solely develop various numerical approximation schemes that approximate the original solution on a predetermined, bounded interval $$J$$.

1. Can you point me to text, or mention results, where approximations schemes are developed for unbounded intervals? I am interested in obtaining an explicit formula for the recurrence equation(s) that determine the approximation (though I do not need a closed-form solution of the recurrence equation(s)).
This is because I have a strong feeling that this discretization is actually equivalent (perhaps after rewriting it in a clever way) to a model that already is well-known. Thus only qualitative information about the solution, such as whether it approaches a periodic orbit, is not sufficient.

2. Conversely, are there (collections of) results (counterexamples) about how fast one runs into trouble, when one takes a standard approximation scheme, valid on $$[-T,T]$$, but where one does not stop the iteration once $$T$$ has been hit after a finite number of steps, but rather where one simply continues iterating?

3. Lastly, if $$f$$ in $$(1)$$ does no depend on $$t$$ and solutions exist globally, it is straight-forward to see how $$(1)$$ induces a flow of a dynamical system. In this case, and supposing that solutions on unbounded intervals exist, any approximation scheme can be thought of inducing a discrete dynamical system that approximates the original. Do you know of any text that treats these objects abstractly?
E.g. by considering a "discretization operator" that maps the flow of a smooth dynamical system such as (1) to its discrete approximation?
All books on numerics of ODE I have consulted so far do not take an abstract view and I think it would be very enlightening if the relationship of a smooth dynamical system and the discrete system that arises as its approximation would by discussed abstractly.
(Then questions such as "given a discrete dynamical system, does there exist a smooth dynamical system such that the discrete one arises as its Euler approximation?" would be very interesting to investigate.)

Any information or help is much appreciated.

• One possible solution (of course not the best) is simply to transform the unbounded interval into a bounded one. May 31, 2021 at 15:14
• @DieterKadelka I guess that would be possible, though it would then not be a priori clear to me what other properties of the original solution such a transform would destroy/preserve. On the other hand, if there is a result that state that "nothing goes wrong" (I 'm myself unsure what a good formalization of this would be), I'd be happy with such an approach. May 31, 2021 at 15:24
• Can you explain what it is that you want to know about the solution behavior over the unbounded interval? Often solutions approach a steady state or periodic orbit, and there are techniques for determining that kind of thing. Also, you will probably get more answers on scicomp.SE (though there are people like me who read both that and this). Jun 1, 2021 at 5:12
• @DavidKetcheson I made an edit to my question to address your question. (And thanks for pointing me towards scicomp.SE, I will keep that in mind for future questions.) Jun 1, 2021 at 12:29
• @someuser I saw you removed your comment; I would have been interested to understand what you meant by your example. Jun 1, 2021 at 13:09

Regarding 1 and 2:

Perhaps the main reason for considering only bounded intervals is that numerical analysts are interested in provably (pointwise) convergent schemes. At least for traditional methods (Runge-Kutta, linear multistep, and many others) the step-by-step nature of the approximation precludes the possibility of uniform convergence on unbounded intervals. That is because the local truncation errors from each step accumulate -- not in a catastrophic way, if the method is stable, but still in general they accumulate and the error grows without bound as the number of steps taken goes to infinity. So uniform pointwise convergence can only be proved for finite intervals.

So it's not that you magically run into trouble when you go past the end of the chosen interval. After all, the bounds of your interval don't enter into the numerical method, and it knows nothing about them. It's simply that your solution will be less accurate at the end of the interval than at the start, and this effect will continue as you integrate further.

There are some areas where very long time integrations are needed, and special techniques are used, for instance in the study of stability of the solar system.

Perhaps a helpful example is Exercise 7.1 from this set, although you should keep in mind that it involves an intentionally poor choice of integrator for the specified problem, so the accumulation of error will not generally be quite so bad.

Of course, there's also the fact that if "one does not stop the iteration ... after a finite number of steps" then your computation will never terminate.

You didn't fully specify the problem, but based on the fact that you mentioned Euler's method, I have assumed that you're dealing with an initial value problem. If this question were about a boundary value problem, that would of course be a completely different matter.

"given a discrete dynamical system, does there exist a smooth dynamical system such that the discrete one arises as its [discrete numerical] approximation)?"

This is often known as modified equation analysis and is a form of backward error analysis. For linear PDEs, it is described for instance in Section 10.9 of LeVeque's book. For Runge-Kutta methods applied to arbitrary initial value ODEs, it is much more difficult and is described most fully in Chapter IX of Hairer, Lubich, & Wanner. I think the last book is also a good example of one that takes the abstract viewpoint you are wishing for.

• Thank you very much for you insight. You were right to assume i had an IVP in mind (I edited my question to clarfiy this as well as other things). I'm completely happy with 2. and 3.; in particular I saw that the book by Hairer et al you linked states in the intro "An important shift of view-point came about by ceasing to concentrate on the numerical approximation of a single solution trajectory and instead to consider a numerical method as a discrete dynamical system which approximates the flow of the differential equation". This was exactly what I had in mind! Jun 1, 2021 at 12:35
• Only regarding 1. I'm puzzled. You said "the step-by-step nature of the approximation precludes the possibility of uniform convergence on unbounded intervals" as " the error grows without bound as the number of steps taken goes to infinity". I'm not 100% this must be true; if, for example, the error would form a convergent series, then surely an approximation on an unbounded interval must be possible? Jun 1, 2021 at 12:49
• If one can somehow prove that for "step-by-step" methods indeed there always exists an example on which they fail --perhaps an interesting question in its own right--- I assume there must be other approximation methods for unbounded intervals (e.g. by constructing a suitable set of "basis functions" and then combining them in some way to approximate the original solution that is sought)? Jun 1, 2021 at 12:50
• >"an approximation on an unbounded interval must be possible" Yes, it's possible for instance using appropriate implicit methods on a problem with an asymptotically stable steady state. But numerical methods are designed to be convergent for general problems. I do happen to know of a recent paper that tries to establish the kind of convergence you are looking for (for a PDE discretization, in fact): arxiv.org/abs/2010.15365 Jun 1, 2021 at 16:40
• I guess your last comment nails it, regarding to my question 1: I am not looking for a numerical method that works for any $f$ on a finite time interval, but rather one where (even if I may have to make a series of additional assumptions on $f$,) I get convergence on unbounded time intervals. The paper you linked states in section 2. "To study numerical schemes for eq. (1) in terms of their long-time, or even infinite time,accuracy, a departure from established numerical analysis notions of convergence of numerical approximations is needed", which is in line with my finding [...] Jun 2, 2021 at 9:46

This is not exactly an answer; but I've done a lot of work on "Euler's Method" for unbounded intervals. There's a bit of a preamble before I can state the theorems though. I forgo all matrix methods, and this isn't to do with numerical approximation, but rather an analytic science. We're going to keep this discussion as elementary as possible; so for that we are assuming all the functions are holomorphic; and we're ignoring some of the finer nuances of Riemann Surfaces.

Begin by defining a functor like object on holomorphic functions. Let $$\mathcal{S}$$ be a domain in $$\mathbb{C}$$; and let $$\phi_j(s,z) : \mathcal{S} \times \mathbb{C} \to \mathbb{C}$$ be a sequence of holomorphic functions. Then, the inner-composition of $$\phi_j$$ is written as,

$$\Omega_{j=n}^m \phi_j(s,z)\bullet z = \phi_n(s,\phi_{n+1}(s,...\phi_m(s,z)))\\$$

These objects arise fairly naturally in the study of first order difference equations (that's difference, not differential); and we can use these objects to approximate; or equal in the limit, a first order differential equation. You can limit these things to infinity, to get what is normally called an "infinite composition"; given you have nice convergence criterion. As I'm going to keep this simple, I'll just state that,

If

$$\sum_{j=1}^\infty |\phi_j(s,z) - z| < \infty\\$$

And this sum converges "compactly normally" on $$\mathcal{S} \times \mathbb{C}$$, then the function,

$$\Phi(s,z) = \Omega_{j=1}^\infty \phi_j(s,z)\,\bullet z\\$$

Is a holomorphic function on $$\mathcal{S} \times \mathbb{C}$$. This can be intuited as the limit,

$$\Phi(s,z) = \lim_{n\to\infty} \phi_1(s,\phi_2(s,...\phi_n(s,z)))\\$$

This is an infinite composition of the first kind, or as a discrete system. This system is great for solving first order different equations. Where if we write $$\phi_j(s,z) = z + q(s-j,z)$$; then the function,

$$Q(s,z) = \Omega_{j=1}^\infty z + q(s-j,z)\,\bullet z\\ Q(s+1) - Q(s) = q(s,Q)\\$$

With this out of the way, we can move onto infinite compositions of the second kind; or the continuous case; which is a slightly more formal construction of Euler's method. I'll keep the conversation on the real line; but these functions should definitely be entire.

Choose an interval (usually I use arcs, so we can think of this as an arc if you want) $$[a,b] \subset \mathcal{S}$$. For now, let's assume this interval is arbitrarily small. Take a descending partition of this interval (arc), $$b=\gamma_0 > \gamma_{1} > \gamma_{2} > ... > \gamma_n = a$$; with sample points $$\gamma_{j} \ge \gamma_j^* \ge \gamma_{j+1}$$. Let $$\Delta \gamma_j = \gamma_{j} - \gamma_{j+1}$$ and let $$||\Delta|| = \max_{0 \le j \le n-1} |\Delta \gamma_j|$$.

Now, let's assume that $$\phi(s,z) : \mathcal{S} \times \mathbb{C} \to \mathbb{C}$$, and is holomorphic. Then we write the "partial compositions" of the integral as,

$$I = \Omega_{j=0}^{n-1} z + \phi(\gamma_j^*,z)\Delta\gamma_j \,\bullet z\\$$

Where the limit as $$||\Delta|| \to 0$$ converges in a specific manner. This is really just Euler's method in disguise, but we've spiced it up with a Riemann-Stieltjes like construction using partitions.

Now, I'm going to skip ahead a bit; because proving this thing converges is very technical; so just bear with me when I say it converges in a nice enough manner to continue this discussion.

I like to call this thing "the compositional integral" for reasons which will appear as follows. Write,

$$Y_{ba}(z) = \int_a^b \phi(s,z)\,ds\bullet z\\$$

Then the following formal laws are always satisfied.

$$\frac{d}{db} Y_{ba} = \phi(b,Y_{ba})\\ Y_{cb}(Y_{ba}) = Y_{ca}\\ Y_{aa}(z) = z\\ Y_{ba}^{-1}(z) = Y_{ab}(z)\\$$

And additionally, we get that Leibniz substitution works. So if $$u(\alpha) = a$$ and $$u(\beta) = b$$ we get,

$$\int_a^b \phi(s,z)\,ds\bullet z = \int_{\alpha}^\beta \phi(u(x),z)u'(x)\,dx\bullet z\\$$

This acts as a strict generalization of Riemann-Stieltjes integration. If I take $$\phi(s,z) = p(s)$$ (so that it is constant in the $$z$$ variable) then we're reduced to the following (just plug in the above definition),

$$\int_a^b p(s)\,ds\bullet z = z + \int_a^b p(s)\,ds\\$$

And the aforementioned group laws just become the usual additivity properties of the integral. So with all this out of the way, we can talk about compositional integration on unbounded domains similarly to integration on unbounded domains. And this will amount to a discussion of "Euler's Method" on unbounded domains (just choose Euler's partition).

So for example; I'll write the first result which brought me to this question; and what I feel is the answer you are looking for.

Let's take,

$$y = \int_{-\infty}^x e^{st}\,ds\bullet t\\$$

Which is the solution to the system of equations,

$$y' = e^{xy(x)}\\ y(-\infty) = t\\$$

So long as you keep $$t > 0$$, this thing will converge. Again, the proof is very involved, but I can sketch it.

The solution to the equation,

$$y_h(x+h) - y_h(x) = he^{xy_h(x)}\\ y(-\infty) = t\\$$

Is always solvable for $$x \in \mathbb{R}$$ and $$h, t>0$$. And its closed form expression is given as,

$$y_h(x,t) = \Omega_{j=1}^\infty t + he^{(x-jh)t}\,\bullet t\\$$

This looks a lot like our integral, but it's being defined on an unbounded domain. So, to prove that our integral converges for $$t>0$$? We just limit $$h \to 0$$. As this will satisfy,

$$\frac{y_h(x+h) - y_h(x)}{h} = e^{xy_h(x)}\\$$

Which is the equation for the derivative if we set $$h\to 0$$. Now, doing this is very very difficult. And the function I just chose here is very manufactured. But, if we view "Euler's method" on unbounded domains, as compositional integration on unbounded domains; there are many more tools at your disposal. And much of them just look like integrations. Doing this in complex scenarios is very difficult though.

But, all in all;

$$\text{`Euler's Method' on unbounded domains} = \text{compositional integrations on unbounded domains}$$

This says nothing about speed of convergence, or how well we approximate; this is, again, an analytic solution to your problem. I hope, maybe, this helps.

Regards,