Euler method (and others) for unbounded intervals

Question

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.

Answer 1

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.

Regarding your 3rd question:

"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.

Answer 2

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,