Computability of Differential Equations

In this chapter, we provide a survey of results concerning the computability and computational complexity of differential equations. In particular, we study the conditions which ensure computability of the solution to an initial value problem for an ordinary differential equation (ODE) and analyze the computational complexity of a computable solution. We also present computability results concerning the asymptotic behaviors of ODEs as well as several classically important partial differential equations.


Introduction
A differential equation is an equation that relates a function with its (partial) derivatives of various orders. Differential equations are widely used in many fields and have been extensively studied since the time when the calculus was invented. Given a differential equation, the goal is to find its solutions. Unfortunately, most differential equations do not possess solutions in explicit forms. In order to obtain useful information about the solutions without having an explicit representation formula for them, various methods and tools have been developed. Classically, the following methods are of most significance: • Numerical methods: where solutions are approximated numerically.
• Qualitative methods: where the behavior of a differential equation and of its solutions is analyzed from a qualitative perspective. Numerical methods generate numerical algorithms which calculate approximations to the exact solution that cannot be found explicitly, and the approximations are expected to converge to the exact solution. The convergence is a fundamental property for a numerical algorithm to be of any use. In numerical analysis, the accuracy of approximation (the convergence rate) is usually not specified, sometimes it is impossible to be specified (see [82] for an example), as an input to the algorithm that computes the approximations, but rather it is provided externally and implicitly by finding certain analytical forms of approximation errors, which frequently depend on additional assumptions; for example, the twice differentiability of the exact solution. This requirement is often crucial, yet not applicable to many initial value problems for the first order ODEs whose solutions may not be C 2 but nevertheless are approximable by convergent numerical solutions.
The qualitative method is another key tool, often complements numerical methods, in the study of ordinary differential equations; in particular, in describing the asymptotic behaviors of solutions of ODEs. For example, it is used to describe the limit set of a typical solution as a spatial object as time tends to infinity.
In addition to the classical methods, there is also a relatively new approach based on the Turing machine model over the real field (see [67], [56], and [87] and references therein): • Computable analysis methods: where computability and computational complexity of problems related to differential equations are studied.
In computable analysis, a problem is called computable if it can in principle be solved or approximated with arbitrary precision by a computing device. The computable analysis method provides useful insights on many problems related to a given ODE from the perspective of computation in the sense that (1) the method can be used to identify whether a problem is computable; (2) if a problem is proved to be computable using the computable analysis method, often the proof would generate a (Turing) algorithm that computes approximations with arbitrary precision and, moreover, with the accuracy specified as a part of the input (computability); and (3) the method has a mechanism for assessing the resources needed for a computation (computational complexity).
The computable analysis method has been used by many authors to explore the computational aspects of differential equations; in particular, to address questions related to computing solutions of differential equations. In this chapter, we present a brief survey of some existing results.
The chapter is organized as follows. Section 0.2 addresses computability of solutions of first order ODEs, while Section 0.3 analyzes computational complexity of the computable solutions. Section 0.4 discusses computability of asymptotic behaviors of ODEs. Finally, section 0.5 presents some computability results of several classically important partial differential equations.

Computability of the Solutions of Ordinary Differential Equations
In this section, we consider ordinary differential equations (ODEs) of the form y = f (y) (0.1) and associated initial-value problems (IVPs) where f : E → R n is continuous over E ⊆ R n . Since an ODE of the more general form y = g(t, y) can be reduced to (0.1), if one replaces the independent variable t by a new variable y n+1 satisfying y n+1 = 1, y n+1 (t 0 ) = t 0 , there is no loss of generality to study only the ODEs in the form of (0.1).

Computability over Compact Sets
Assume that E is a compact subset of R n . Classically, it is well known that continuity of f ensures the existence of solution(s) to the IVP (0.2) (Peano's existence theorem, see e.g. [11, pp. 191-193]), and if f is Lipschitz continuous then the solution is unique and C 1 smooth (Picard-Lindelöf theorem, see e.g. [43, pp. 8-10]). The compactness of E is essential for the proof of Peano's existence theorem in which a sequence of functions is constructed and then the Arzelà-Ascoli theorem is employed to obtain a subsequence that converges toward a solution of (0.2). The proof is non-constructive. On the other hand, the classical proof of existence and uniqueness for Lipschitz continuous f is constructive, such as Picard's scheme of successive approximations; such constructive schemes lead naturally to some computer programs for computing the unique solution, uniformly in f . Recall that f is said to be Lipschitz continuous if it satisfies a Lipschitz condition on E; i.e. there exists some positive constant K (a Lipschitz constant) such that, for all x, y ∈ E Thus the question remaining to be investigated is whether solutions are still computable when the vector field f is continuous on E but fails to be Lipschitz continuous. In this case, whether the IVP (0.2) has a unique solution plays a pivotal role in determining if (0.2) admits a computable solution. A theorem by Osgood [10] shows implicitly that the initial-value problem y = f (t, y), y(0) = y 0 , does admit a computable solution on [0, δ ] for some δ > 0, provided that f is computable on E = [0, 1] × [−1, 1] and the solution of the IVP is unique. The idea used for the proof is to compute sequences of functions which converge to a maximal solution y max from above and a minimal solution y min from below. If the solution is unique, then the convergence is effective and thus results in a computable solution. However, this technique only works in the plane, since in higher dimensions there is no similar notion of a maximal or a minimal solution. In [71] Ruohonen generalizes this result to higher dimensions, using approximation funnels to compute the solution of (0.2), obtaining the following result.
Theorem 0.2.1. Suppose that f is computable and not identically zero on the closed ball B(y 0 , b) = {x ∈ R n : |y 0 − x| ≤ b}, where b > 0,t 0 , y 0 are computable. Then the solution of (0.2) is computable on [t 0 ,t 0 + δ ] for δ = b max x∈E f (x) > 0, provided that the solution of (0.2) is unique on that interval.
Another result that relies on Picard's construction to compute the solution of (0.2) is presented in [28]. There a domain theoretic version of Picard's operator is presented to compute the solution of the IVP (0.2), provided that the function defining the IVP is Lipschitz continuous.
It is also shown in [19] that, under certain conditions, solutions of the class of implicit ODEs in the form of A(y) · y = f (y, λ ) are computable, where A(y) is an n × n matrix formed by polynomials in y 1 , . . . , y n , with y = (y 1 , . . . , y n ), and λ is a parameter.
The scenario changes completely if (0.2) admits multiple solutions. In this case, it is possible that none of the solutions is computable. The following result was first proved by Aberth [2, p. 152] for the case where f is assumed just to be computable; later a different proof was presented by Pour-El and Richards [65]; then the result was extended in [55] to include the case where f is polynomial-time computable.
Theorem 0.2.2. There exists a polynomial-time computable function f : does not have a computable solution y on [0, δ ] for any δ > 0.
The idea for the proof is to construct a function f that embeds a non-computable problem P and show that if a computable solution for (0.3) exists, then P can be solved computationally, a contradiction. For details of the proof, we refer the reader to the original papers [65], [55] or to [56, pp. 216-219].

Computability over Non-Compact Sets
In this subsection, we assume that E is an open subset of R n . The Peano existence Theorem indicates that the IVP (0.2) has a solution existing on some time interval containing t 0 in its interior, provided that f is continuous on E. An interval is called an interval of existence if it contains t 0 and if a solution of (0.2) exists on it. The following classical result indicates how large an interval of existence may be (see e.g. [43, pp. 12-13]).
Proposition 0.2.3. (Maximal Interval of Existence) Let f : E → R n be continuous on E. Then the IVP (0.2) has a maximal interval of existence, and it is of the form (α, β ), with α ∈ [−∞, ∞) and β ∈ (−∞, ∞]. There is a solution y(t) of (0.2) on (α, β ) and y(t) leaves every compact subset of E, i.e. y(t) tends to the boundary ∂ E of E, as t → α from above and as t → β from below.
The proposition is a simple classical existence statement "the IVP (0.2) has a maximal interval of existence and a solution exists on it." It turns out that this simple statement gives rise to a variety of questions about computability. Is the maximal interval (α, β ) of existence computable from the data ( f ,t 0 , y 0 ) that defines the IVP (0.2)? Does (0.2) admit a computable solution on (α, β )? If E = R n , then a solution of (0.2) is said to blow up in finite time if α and β are two real numbers. Is the problem whether (0.2) admits a finite blowup solution decidable? We will discuss these questions in this subsection.
It is known from the subsection 2.1 that if the IVP (0.2) has multiple solutions then it may occur that none of the solutions is computable. Hence, in this subsection, we consider only those IVPs to which the solution is unique.
We begin by examining the proof of the proposition. The proposition can be proved by using the Peano existence theorem repeatedly: The theorem first produces a solution y(t) of the IVP (0.2) on an interval [a 0 , b 0 ] containing t 0 ; then it extends , v(a 0 ) = y(a 0 ); and the process is repeated infinitely many times to reach the maximality: Thus, if the IVP (0.2) has a unique solution y(t) over (α, β ), then it follows from Theorem 0.2.1 that a j , b j , and y(t) on interval [a j , b j ] are computable from j, f , t 0 , and y 0 . But the argument is insufficient to show that y(t) is computable on (α, β ) because it is possible that there does not exist a "master" algorithm computing y(t), independent of j. One possible way of dealing with this shortcoming is to strengthen f to be a smoother function, which is actually an approach frequently used in classical analysis of ODEs: the smoother the vector field f is, the better, possibly, the solution y behaves. This is the approach used in [36]. By strengthening f to be a C 1 function on E, the following result is shown in [36].
Theorem 0.2.4. Assume that E is a r.e. open subset of R n . Consider the initial value problem (0.2) where f is C 1 on E. Let (α, β ) be the maximal interval of existence of the solution y(·) of (0.2) on E. Then: The operator ( f ,t 0 , y 0 ) → (α, β ) is semi-computable (i.e. α can be computed from above and β can be computed from below).
Recall that E is a r.e. open subset of R n if it can be filled up by a computergenerated sequence {B j } ∞ j=1 of rational balls (balls having rational centers and rational radii), i.e., E = ∞ j=1 B j . The main role played by the C 1 smoothness in the proof is that C 1 smoothness implies local Lipschitz continuity (and the converse is false); that is, if f is C 1 on E, then f is Lipschitz continuous on every compact subset of E. An algorithmic version of this fact is proved in [36]. Thus, when giving a C 1 -name of f , it becomes possible to compute a Lipschitz constant K j for each E j , where E j is the closure of j i=1 B i , then compute the solution of (0.2) on E j , and thus compute the solution over the maximal interval (α, β ) of existence, because E j are getting ever larger and eventually fill up E.
Although for many important ODEs their vector fields f are indeed C 1 , the requirement for f being C 1 is nevertheless not necessary for ensuring the existence of a unique solution. It is then natural to ask whether it is possible to compute the unique solution of the IVP (0.2), whenever it exists, over the maximal interval (α, β ); in other words, whether Theorem 0.2.4 holds true over (α, β ) when E is an open subset of R n . The problem is studied in [26] and a positive answer is presented there. The proof uses a quite different approach; the idea underlying the approach is to try to cover the solution with rational boxes and to test (i) and (ii), in an algorithmic way: (i) if a given set of rational boxes is an actual covering of the solution of (0.2), and (ii) if the "diameter of the covering" is sufficiently small so that the rational boxes provide an approximation of the solution with the desired accuracy, whenever (i) is satisfied. The process of enumerating all possible families of rational boxes and applying the tests (i) and (ii) to each family generates, effectively, better and better approximations to the unique solution of (0.2), and thus proves that the solution is computable. This procedure can also be applied to study the computability of differential inclusions (see [26] for more details). The following theorem proved in [26] (see also [25]) strengthens Theorem 0.2.4. The operator ( f ,t 0 , y 0 ) → y(·) is computable; 2. The operator ( f ,t 0 , y 0 ) → (α, β ) is semi-computable (i.e. α can be computed from above and β can be computed from below).
In particular, if f is a computable function and t 0 , y 0 are computable points, then (α, β ) is a recursively enumerable open set and the solution y(·) is a computable function.
A question that remains is whether the maximal interval of existence is computable. The answer is negative, as shown in [36].
Theorem 0.2.6. There is an analytic and computable function f : R → R such that the unique solution of the problem is defined on a non-computable maximal interval of existence.
The theorem indicates that the computability of the maximal interval (α, β ) of existence is not tied up with the smoothness of f ; rather it has more to do with the "globalness" of (α, β ). Later, in section 4, we will see again that a "global property" is usually difficult to compute. The counterexample is constructed by creating a computable odd and bijective function ϕ : (−α, α) → R, where α is a noncomputable real number, such that ϕ satisfies the following conditions: (i) ϕ is the solution of (0.4) for an analytic and computable function f ; and (ii) ϕ(x) → ±∞ as x → ±α ∓ .
When E is the whole space R n and f is defined on R n , if a solution y(t) of (0.2) is defined for all t ≥ t 0 , it is called a (positively) global solution that does not blow up in finite time. In general, it is difficult to predict whether or not a solution will blow up in finite time from a given initial datum (t 0 , y 0 ), as Theorem 0.2.6 already indicates, because it often requires extra knowledge on some quantitative estimates and asymptotics of the solution over a long period of time. Actually the problem of determining whether the solution of (0.2) blows up in finite time is undecidable even if f : R n → R n has only as components polynomials [39]. Nevertheless, certain computational insight on the blow up problem is obtainable from f as the following theorem in [69] suggests.
Theorem 0.2.7. Consider the IVP (0.2), where f is locally Lipschitz. Let Z be the set of all initial values (t 0 , y 0 ) for which the corresponding unique solution of (0.2) is (positively) global. Then Z is a G δ -set and there is a computable operator determining Z from f .

Computational Complexity of the Solutions of Ordinary Differential Equations
In this section we analyze the computational complexity of solving an initial value problem (0.2). Note that some of the topics discussed in this section overlap with the topics explored in the chapter "Uniform Computational Complexity in Analysis" of this handbook, which we strongly recommend to the reader. We begin by noting that, as mentioned before, any C 1 function defined on a compact set E ⊆ R n satisfies a Lipschitz condition. The analysis of any standard algorithm (e.g. Euler's method, or the algorithm implicit in Picard's iteration method) to solve an initial value problem (0.2) applied to a continuous function which satisfies a Lipschitz condition over a compact set, with a Lipschitz constant K > 0, shows that the value of the Lipschitz constant affects the number of iterations needed to compute the solution of (0.2) with a given accuracy. However, since the Lipschitz constant is constant in a compact, it can be hidden in the big O notation. When we solve (0.2) over a non-compact set, we typically have infinitely many distinct local Lipschitz constants, and the complexity results need to take into account this effect. This is the reason both cases are analyzed in separate sections, like for the case of computability.

Results for Compact Sets
We have just discussed above the importance that a Lipschitz condition seems to have when computing the solution of an initial value problem (0.2). First it makes sense to ask how the presence of a Lipschitz constant affects the computational complexity of the solution of (0.2). We already know from Theorem 0.2.1 that if (0.2) has a unique solution, then it must be computable. So what is the complexity of the solution of (0.2) if we assume that f is "easy to compute" (polynomial-time computable)? The following result from [56, p. 219], based on another result of Miller [58, p. 469] shows that in that case the solution of (0.2), while computable, can have arbitrarily high complexity. Note that in (0.5) we have used the independent variable t explicitly, i.e. we considered a non-autonomous ODE instead of an autonomous ODE as in (0.2). Technically speaking, that is not needed as we have seen in Section 0.2, since a nonautonomous IVP can be converted into an autonomous IVP. However, in the remainder of this section we will explicitly consider the independent variable t since, as we will see later in this section, this will allow us to do a more refined analysis on what happens at a computational complexity level, by analyzing how the smoothness of f relative to t or to y can affect the complexity of the solution of (0.5). In this section we also suppose that the initial condition for (0.5) is y(0) = 0, although WLOG any initial condition y(t 0 ) = y 0 , where t 0 , y 0 are polynomial-time computable, could be used to obtain the same results. We note that for non-autonomous IVPs (0.5) the Lipschitz condition is not needed for the variable t. Therefore, in that case, we say that f : Lipschitz, then when we convert the ODE of the initial-value problem (0.5) into an autonomous ODE (0.1), then f in (0.1) is Lipschitz continuous. Therefore using non-autonomous ODEs with (right-)Lipschitz ODEs or autonomous Lipschitz continuous ODEs are equivalent approaches, although the non-autonomous case (0.5) has the advantage mentioned above.
In (0.5) we also assume that the function f is defined on [0, 1] × B(0, 1). However those results can easily be extended to functions defined over [0, T ] × B(0, M) for T, M > 0 via rescaling. More concretely let f : [0, T ] × B(0, M) → R n and let y : [0, T ] → B(0, M) be the solution of (0.5). Note that we suppose that y(t) is always defined for all t ∈ [0, T ] (and hence is always inside B(0, M)). Define the function y * : [0, 1] → R n as y * (t) = y(tT ). In other words y * (1) = y(T ). It is not difficult to see that y * is the solution of the IVP y * = T f (tT, y * ) and y * (0) = 0.
Similarly one can define the function y * * : [0, 1] → R n as y * * (t) = y * (t)/M = y(tT )/M. In other words y * * (t) is just like y * , but scaled down by a factor of M so that y * * (t) ∈ B(0, 1) for all t ∈ [0, 1]. It is also not difficult to check that y * * is the solution of the IVP In [23, p. 450] it is shown that if f of (0.5) satisfies a Lipschitz condition on As mentioned by Ko, the proof of this result follows from a careful analysis of Euler's method. For example, it is shown in [7, pp. 349-350] that when n = 1 and we use Euler's method to compute the solution y of (0.4) over the time interval [0, 1] with time step bounded by h > 0 and rounding error bounded by ρ > 0, the error e n in step n is bounded by where A, B,C > 0 are some constants (the original formula (6.2.32) [7] explicitly tells how A, B,C depend on the Lipschitz constant, etc., but that level of detail is not needed for our analysis). So to compute the solution of (0.4) over the time interval [0, 1] with accuracy 2 −n we could use Euler's method with h = a2 −n and ρ = b(2 −n ) 2 for some appropriate constants a, b > 0 which are independent of n.
Euler's method would then run in space polynomial (indeed quadratic) in n which shows the result. The case n > 1 can be obtained by considering the Euler method for dimensions n > 1 as explained in [11,Section 7.2]. The analysis of the error (done in [11,Section 7.3], except that the rounding error is not considered, but this can also be taken into account using a procedure similar to that of [7, p. 350]) gives a formula (0.6), but with constants A, B,C > 0 which have different values than from the previous case, and therefore the conclusions remain valid for the case where n > 1.
In [56,Section 7.4] (also implicitly in [55, p. 159, Question D],) Ko asks whether polynomial space is also a lower bound for the complexity of the solution y of (0.5), assuming the conditions of the previous theorem for n = 1. Kawamura provides an affirmative answer in [49] (see also [48]) using an appropriate notion of reduction (see [49,Section 2.2]) which allows us to define C -hard and C -complete functions for a complexity class C . Namely Kawamura shows the following theorem.
The previous results are non-uniform since we fix the complexity of the input to be polynomial time. But one can naturally ask what is the (uniform) complexity of the operator LipIV P which maps a function f satisfying a Lipschitz condition to the solution y of (0.5). A major problem to tackle is that until recently there was no general theory to measure the complexity of operators which map functions into functions. This problem was solved in [50] (see also the chapter of this handbook mentioned in the beginning of this section) using appropriate representations of functions and appropriate notions of reduction, hardness, and completeness. Using this setting one can define C -hard and C -complete functions for a complexity class C . Since we are talking of functions, it makes sense to take C = FP (functions computable in polynomial time) or C = FPSPACE (functions computable in polynomial space). By using appropriate representations δ for continuous functions in C([0, 1]) and δ L for functions satisfying a Lipschitz condition (which are continuous), and an appropriate reduction ≤ 2 mF (see [50] for more details) it was also shown in [50,Theorem 4.10] that the operator LipIV P is FPSPACE-complete.
An interesting question is whether smoothness of f helps reduce the computational complexity of solving an IVP (0.5). We say that f : if the partial derivative ∂ n+m f (t, y)/∂t n ∂ y m exists and is continuous for all n ≤ i and m ≤ j; it is said to be of class C (∞, j) if it is of class C (i, j) for all i ∈ N. In [51] this question is analyzed and Theorem 0.3.3 is generalized to C (∞.1) functions.
Theorem 0.3.6. There exists a polynomial-time computable function f : In that paper it is also shown [51,Theorem 2] that if f is more than once differentiable, then the unique solution can be CH-hard, where CH ⊆ PSPACE is the counting hierarchy.
Theorem 0.3.7. Let k be a positive integer. There exists a polynomial-time computable function f : In the most extreme case, when f is analytic, the solution of (0.5) is polynomialtime computable on a compact set. This follows from results of Müller [60] and Ko and Friedman [57], which show that polynomial-time computability of an analytic function on a compact interval is equivalent to polynomial-time computability of the sequence of its Taylor coefficients at a rational point. Given a sequence of Taylor coefficients for f one can compute the Taylor coefficients for the solution y of (0.5) [61, Theorem 2.1]. Note that it is known that if f is analytic, then so is the solution y of (0.5) [5,Section 32.4]. By analytic continuation [60] polynomial time computability of y follows. As remarked in [51, last paragraph of Section 5.2], the previous argument also shows uniform polynomial computability of the operator LipIV P, if we represent analytic functions by their sequence of Taylor coefficients, obtaining a representation δ Taylor , since the Taylor sequence of y is easy to compute from the Taylor sequence of f . In the next theorem LipIV P D represents the operator LipIV P restricted to the class of (real) functions D and C ω represents the class of analytic functions.
Theorem 0.3.6 can also be extended to an uniform version by showing that the operator LipIV P C (∞,1) [51, last paragraph of Section 5.2] is FPSPACE complete, using the polynomial-time Weihrauch reduction ≤ W .
Computability in polynomial time of the solution of (0.5) is also obtained in [29] by using polynomial enclosures in Picard's method, provided that some assumptions are made to the function f in (0.5), such as that f is Lipschitz and that its representation via polynomial enclosures satisfies certain properties (e.g. this representation should not take "too much space").

Results for Non-Compact Sets
In the previous section we have analyzed the computational complexity of solutions of an IVP (0.5) (or, equivalently, (0.2)) where f is considered over a compact set But what happens when f : R n+1 → R n and the solution y of (0.5) is considered over its maximal interval of definition? That is, what is the computational complexity of y : R → R n , or more generally of the operator mapping f to y?
As a first attempt we could consider the approach of, given t ∈ [0, +∞), to rescale the IVP (0.5) given by f : [0,t] × B(0, M) → R n , where M is big enough so that y(t) ∈ B(0, M) for all t ∈ [0,t], to an IVP (0.5) which uses a function f : [0, 1] × B(0, 1) instead of f , as explained in the previous section. Since, for a fixed M (and fixed T , with t ≤ T ), the complexity of solving (0.5) is the same as the complexity of solving (0.5) wheref is used instead of f , we might be led to conclude that the complexity of obtaining the solution of (0.5) for the noncompact case where f : R n+1 → R n is the same as when f is taken over a compact set [0, 1] × B(0, 1). However this reasoning is incorrect, since M (and T ) is not fixed in this case and must be uniformly computed from y, i.e. M = M(y(t)) = M(t). This is not a trivial task since to compute y(t), we need to know M(t). On the other side, to know M(t), we need to know max 0≤u≤t y(u) and therefore we get into a circular argument where M(t) is needed to compute y(t) and vice versa. To better illustrate this problem, consider the following system taken from [64] . . . y d (t)= y 1 (t) · · · y n (t) The results of Section 0.3.1 show that for any fixed, compact I = [0, T ], y is polynomial-time computable since the ODE is analytic. On the other hand, this system can be solved explicitly and yields: One immediately sees that y d being a tower of exponentials prevents y from being polynomial-time computable over R, since one might need (supra-)exponential time just to write down an integer approximation of y d with precision 1/2 = 2 −1 for d ≥ 2. This example shows that the solution of an analytic IVP (or even of a polynomial IVP) can be polynomial time computable on any fixed compact, while it may not necessarily be polynomial time computable over R.
A possible way to solve this problem is to analyze the complexity of (0.5) using a bound on the growth of the solution y as a parameter on the function used to measure the complexity (this is an example of parametrized complexity), since the problem of knowing how quickly y can grow is not generally well-understood, even when f is constituted by polynomials. This approach is taken in [13] for the case of polynomial IVP (solutions of polynomials IVPs (0.5) where f is formed by polynomials are sometimes called PIVP functions) and in [12] for the case of analytic functions. The motivation for studying PIVPs is that the class of PIVP functions is well-behaved and includes many interesting functions. In particular it contains all of the usual functions of analysis (exponential, trigonometric functions, polynomials, their inverses) and is closed under the usual arithmetic operations, composition and ODE solving [21], [34]. It can also be shown that this class is exactly the class of functions generated by Shannon's General Purpose Analog Computer (GPAC) [77], [35], which is the mathematical model for analog computers (differential analyzers) used before the advent of digital computers [20] (see also the chapter "Analogue Models of Computation on the Real Numbers" of this handbook). The following result is from [13] and shows that the solution of a polynomial IVP can be computed in polynomial time with respect to a bound of Y (t) = max 0≤u≤t y(u) .
Theorem 0.3.11. There exists an algorithm A such that for any vector of polynomials p with polynomial time computable coefficients, µ ∈ N, t ∈ Q with t ≥ 0, and Y ∈ Q such that Y max 0≤u≤t y(u) , satisfies where y is the solution of (0.5). Furthermore A (p, µ,t,Y ) is computed in time polynomial in the value of µ,t and Y .
This result can be extended to IVPs where f is not necessarily a polynomial, assuming polynomial-time computability of the higher derivatives of f and an appropriate (polynomial) bound on the growth of those derivatives (see [13,Section 7]).
We note that Theorem 0.3.11 can be extended to the case where the initial condition is of the form y(t 0 ) = y 0 , where t 0 ∈ Q and y 0 is polynomial-time computable (and it is actually stated in this way in the original paper [13]). This result is proved by using a variable order method in which the solution of (0.5) is computed over a succession of subintervals [a i , a i+1 ], with i = 0, 1, 2, . . . and a 0 = 0, whose union gives the maximal interval of definition of the solution y of (0.5). On each subinterval a Taylor approximation of the solution is computed, but using a variable order method: instead of using an approximation of fixed order for each subinterval [a i , a i+1 ], the order of the approximation is allowed to change on each interval. Note that the usual methods for numerical integrations (including basic Euler's method, etc.) fall in the general theory of n-order methods for some n. It was already suggested in [79, Section 3] that fixed order methods might only run in exponential time when solving (0.5) for certain polynomial-time computable functions, but that variable order methods might solve the same problem in polynomial time. Variable order methods have also been used in certain contexts to solve IVPs, but usually without complexity results or with complexity results valid only for compact sets, see e.g. [27], [45], [8], [1].
The previous result has the drawback that we need to know a bound Y max 0≤u≤t y(u) (which is used as input) to be able to compute y(u) for u ∈ [0,t].
This is improved in [64], where T and Y are replaced by a single parameter -the length of the solution curve y, with the added benefit that this parameter is not needed as an input to the algorithm. We recall that the length of the curve defined by the graph of a function f between x = a and x = b is In the case of the solution of (0.5) where f = p is a vector of polynomials (i.e., each component of the vector p is a polynomial), we note that the derivative of the solution y is given by p(y). If the degree of p is k (in the case where p has more than one component, the degree of p is the maximum of the degrees of its components), the length of the solution has a value which has an order of magnitude similar to the following quantity where Σ p denote the sum of the absolute values of all the coefficients of p (or the maximum of such sums made over each component of p, if p has more than one component). The following theorem is from [64]. where Len(t 0 ,t) = t t 0 Σ p max(1, y(u) ) k du, k is the degree of p, n is the number of components of p, and poly(a 1 , . . . , a j ) means polynomial in a 1 , . . . , a j .
More recently [85], [52] a similar result was established for IVPs with the form (0.2), where f is analytic, using techniques similar to those of [13] (i.e. using power series to compute the solution of the IVP). The authors also characterize the complexity of solving (0.2) using parametrized complexity. In particular, if f is polynomial-time computable in several input quantities (such as a bound for a complex extension of f ; the inverse of the distance to a singularity; a bound on the norm of the point where f is evaluated; the precision up to which f must be calculated), then the solution of the IVP (0.2) is also polynomial-time computable on those inputs. See [85,Theorem 4.4.1] for more details.
In [85], [54] the above result is applied to some classes of volume preserving systems such as Hamiltonian systems. In particular, for some classes of such systems, it can be shown that their average time complexity is polynomial-time computable, provided that hard instances are relatively rare. The authors also apply this result to show that the planar circular restricted three-body problem is polynomial-time computable on average under certain assumptions.

Computability of Qualitative Behaviors of Ordinary Differential Equations
Qualitative analysis of asymptotic (long term) behaviors is a key tool for exploring ordinary differential equations because most ODEs do not have representation formulas for their solutions and numerically computed solutions are usually limited for short terms only. In some circumstances, in particular for hyperbolic systems, qualitative analysis may also provide algorithmic descriptions, to various degrees, of certain asymptotic behaviors. There are only few results of such nature though, which are mainly about computation of trajectory behaviors near a hyperbolic equilibrium point. We note that a statistical approach can also be used to obtain information about the asymptotic behavior of a dynamical system. However, since this approach involves measure theory, which is discussed in its own dedicated chapter of this book, results of this nature are not discussed in this chapter.
The asymptotic behavior of a differential equation is captured by its limit sets, which are the states the solutions reach after an infinite amount of time has passed. A limit set can be a point, a finite set of points, a curve, a manifold, or even a complicated set with a fractal structure known as a strange attractor. Limit sets are well understood in dimensions one and two: in dimension one and for continuous f in (0.1), the only possible limit sets are equilibrium points; and in dimension two and for C 1 smooth f , a closed and bounded limit set other than a periodic orbit or an equilibrium point consists of equilibria and solutions connecting them according to the Poincaré-Bendixson Theorem, a celebrated result in the field of dynamical systems. Structures of limit sets in dimensions greater than two become much richer and more complex; for instance, strange attractors abound in R 3 .
An equilibrium point is the simplest type of limit sets. has zero real part, where D f (ỹ) is the Jacobian of f atỹ. If the real part of every eigenvalue of D f (ỹ) is negative (positive, respectively), thenỹ is called a sink (source, respectively).
The structure of solutions, also known as trajectories or orbits, near a hyperbolic equilibrium is characterized by the Stable Manifold Theorem, which is one of the most important results concerning the local qualitative theory of differential equations. From the perspective of computation, hyperbolicity ensures that D f does not have a computational singularity nearỹ, an indication of possible existence of an algorithmic characterization for the orbit structure nearỹ. Note that D f presents the directions of the trajectories. We recall that the condition x = 0 for computable real numbers is considered a computational singularity because it cannot be decided effectively. The following theorem is an effective version of the Stable Manifold Theorem, which demonstrates that the trajectory structure near a hyperbolic equilibrium is indeed computable.
Moreover, if k < n (k > 0), then a rational number η and a ball D centered atỹ can be computed from f such that for any solution φ (t, y 0 ) to the equation (0.1) with , respectively) no matter how close the initial value y 0 is toỹ.
The classical proof of the stable manifold theorem relies on the Jordan canonical form of A = D f (ỹ). To reduce A to its Jordan form, one needs to find a basis of generalized eigenvectors. Since the process of finding eigenvectors from corresponding eigenvalues is not continuous in general, it is a non-computable process. Thus if one wishes to construct an algorithm that computes some S and U of (0.1) at y 0 , a different method is needed. In [37] an analytic, rather than algebraic, approach to the eigenvalue problem is used to allow the computation of S and U without calling for eigenvectors. The analytic approach is based on function-theoretical treatment of resolvents (see, e.g., [84], [46], [47], and [70]).
The local stable and unstable manifolds can be extended to the global stable and unstable manifolds by taking respectively. Is the global stable or unstable manifold computable for a computable hyperbolic equilibrium point? Since the global stable manifold of a sink s coincides with the basin of attraction of s, a less ambitious question is whether the basin of attraction of a computable sink is computable? The answer depends on the function f in (0.1). It is known that for hyperbolic rational functions, there are (polynomialtime) algorithms for computing basins of attraction and their complements (Julia sets) with arbitrary precision [9]; in other words, basins of attraction and Julia sets of hyperbolic rational functions are (polynomial-time) computable. On the other hand, for more general but still well behaved functions f , the answer can be negative, even for computable and analytic functions (see [59], [100] for examples of noncomputable basins of attraction occurring in C k -and C ∞ -systems). The following theorem (taken from [33]) shows that the basin of attraction of a computable sink in a computable and analytic system can be non-computable.
Theorem 0.4.3. There exists an ODE (0.1) defined by an analytic and computable function f , which admits a computable hyperbolic sink s such that the global stable manifold (basin of attraction) of s is not computable.
Although the basin of attraction constructed in the theorem above is not computable, it is an r.e. open subset of R n (see [100]); in other words, it can be filled up by a computable sequence of open balls. In general, it can be shown (see [37,Section 5]) that the degree of unsolvability of computing a global stable manifold is essentially Σ 0 2 .
The proof of Theorem 0.4.3 is fairly complicated. It starts by encoding a wellknown non-decidable problem into the basin of a computable hyperbolic sink in a discrete-time system; then this discrete-time system is embedded into a continuoustime system. The standard suspension method (see Smale [78], Arnold and Avez [6]) for embedding a discrete-time system into a continuous-time system is not sufficient for the construction; instead, a new method is developed for the embedding.
Another problem that is close in spirit to the stable manifold theorem concerns linearization of the flow of a nonlinear ODE (0.1) near a hyperbolic equilibrium point. Poincaré originated the study and, later, a number of researchers contributed to progress on related problems which resulted in the Hartman-Grobman theorem (see e.g. [63, Section 2.8]); the theorem asserts that near a hyperbolic equilibrium pointỹ of a nonlinear ODE (0.1), there is a homeomorphism H such that H • φ = L • H, where φ is the solution to the ODE (0.1) and L is the solution to its linearizationẏ = D f (ỹ) y. In other words, nearỹ, the ODE (0.1) has the same qualitative structure as its conjugated linear system. This theorem remains important, since it shows the structural stability of hyperbolic equilibria in sufficiently smooth dynamical systems. The computability study of this linearization is motivated by one of the seven open problems in the addendum to the book Computability in Analysis and Physics by Pour-El and Richards [67], who asked, "What is the connection between the computability of the original nonlinear operator and the linear operator which results from it?" It is proved in [38] that the homeomorphism H that performs the linearization is computable. The precise statement is given below. [Effective Hartman-Grobman theorem [38]] There is a computable map Θ : F → O × O × C(R n ; R) × C(R n ; R n ) such that for any f ∈ F , f → (U,V, µ, H), where (a)H : U → V is a homeomorphism ; (b)the unique solution y(t, y 0 ) = y(y 0 )(t) to the initial value problemẏ = f (y) and y(0) = y 0 is defined on (−µ(y 0 ), µ(y 0 )) for all y 0 ∈ U; moreover, y(t, y 0 ) ∈ U for all y 0 ∈ U and −µ(y 0 ) < t < µ(y 0 ) ; (c)H(y(t, y 0 )) = e D f (0)t H(y 0 ) for all y 0 ∈ U and −µ(y 0 ) < t < µ(y 0 ) .
Recall that for any y 0 ∈ R n , e D f (0)t y 0 is the solution to the linear probleṁ y = D f (0)y, y(0) = y 0 . So the theorem shows that the homeomorphism H, computable from f , maps trajectories near 0, a hyperbolic equilibrium point, of the nonlinear problemẏ = f (y) onto trajectories near the origin of the linear probleṁ y = D f (0)y. In other words, H is a conjugacy between the linear and nonlinear trajectories near the origin. Since the classical proofs of the Hartman-Grobman linearization theorem are not constructive, the effective version of the theorem cannot be obtained from a classical proof.
The results discussed so far are centered at hyperbolic equilibrium points. Evidently, hyperbolicity and local (in the spatial sense) are two key conditions for obtaining algorithmic descriptions about asymptotic behaviors of trajectories of the ODE (0.1) near an equilibrium; on the other hand, the role played by the structure f , near an equilibrium, of the ODE (0.1) turns out not to be as crucial, since the computations for a local stable manifold, a local unstable manifold, and the conjugacy H at a hyperbolic equilibrium are all uniform in f . Equilibrium points are the simplest limit sets, and there are many more limit sets other than equilibrium points. What about computability of other limit sets; in particular, those with complicated fractal structures such as strange attractors? The question is largely open. There is however one type of strange attractors -geometric Lorenz attractors -that has been proved to be computable recently. The Lorenz attractor is perhaps the most famous strange attractor. It was introduced in 1963 by E. N. Lorenz as one of the first examples of strange attractors (see e.g. [41], [44] for a treatment of the Lorenz attractor and [4] for a more in-depth analysis of Lorenz-like attractors). However, Lorenz' research was mainly based on (nonrigourous) numerical simulations and, until recently, the proof of the existence of the Lorenz attractor remained elusive. To address that problem, Afraimovich, Bykov, and Shil'nikov [3], and Guckenheimer and Williams [42] originated the study of flows satisfying a certain list of geometric properties based on the behavior observed in the numerical simulations of the Lorenz equation. In particular, they proved that any such flow must contain a strange attractor. These examples came to be known as geometric Lorenz models, and the strange attractor contained in a geometric Lorenz flow is called the geometric Lorenz attractor. In 2002 Tucker [86] provided, using a combination of normal form theory and rigorous numerics, a formal proof on the existence of the Lorenz attractor by showing that the geometric Lorenz models do indeed correspond to the Lorenz system for certain parameters. Since a geometric Lorenz model supports a strange attractor, so does the Lorenz system. Recently, the following result concerning computability of geometric Lorenz attractors is proved in [40].
Theorem 0.4.5. For any geometric Lorenz flow, if the data defining the flow are computable, then its attractor is a computable subset of R 3 . Moreover, the physical measure supported on this attractor is a computable probability measure.
We should also note that some related problems can be found in verification theory. For example, in [24] the reachability problem is analyzed from a computability perspective. The reachability problem consists of finding the set of all points which can be reached in finite time from trajectories starting in an initial set. It is shown in [24] that the reachable set is lower-computable for ODEs (or, more generally, differential inclusions). However computability only holds in certain conditions.

Computability of Partial Differential Equations
The area of partial differential equations is vast and there is no general theory known concerning the solvability of all PDEs. In fact, most PDEs do not even have classical solutions (a solution of a PDE of order k is classical if it is at least k times continuously differentiable in the uniform norm). Instead, they may only be solvable in the sense of weak or generalized solutions defined in some function spaces equipped with an energy-type norm such as Sobolev spaces or in the space of generalized functions. Computability theories of these spaces are established in [98], [101]. Many modern approaches to the subject -energy methods within Sobolev spaces, the calculus of variation, conservation laws, semigroups, etc -are also useful tools in the recursion theoretic study of the subject; in particular, energy methods augmented by a fix point argument.
An equation is exactly solvable if there is an explicit formula for its (classic or weak) solution. There are only few solvable partial differential equations which are exactly solvable. Pour-El and Richards have studied three important second-order linear exactly solvable PDEs from the viewpoint of computability [67]: the heat equation, Laplace's equation, and the wave equation. They showed, by exploiting the solution formulas, that the heat equation and Laplace's equation admit directly computable classical solutions, at least in the case where good domains such as rectangular regions with computable corners are considered. For the wave equation, they demonstrated that a computable initial datum may generate a non-computable wave propagation measured in the uniform norm (also see [66], [68]. A discussion of those results can also be found in [31]); but in a more physically relevant setting, the wave equation does admit computable solutions in the energy norm (see, for example, [92], [88], [91]). This study indicates that computability of the solution operator of a PDE, if the PDE is uniquely solvable, depends not only on the structure of the equation but also on the physical measurements used for the initial/boundary function(s) and for the solution. Another notable fact revealed by the study is that the uniqueness of solution of a computable initial value problem for a PDE doesn't guarantee the solution being computable -a significant distinction between computability theories of ODEs and PDEs. A related result can be found in [22], where the difficulty of solving (ordinary and partial) differential equations is studied using index sets. Several cases of differential equations are considered (e.g. equations with locally computable solutions) and are shown to correspond to certain levels in the arithmetical hierarchy, which depend on the particular problem being considered. This type of results provides more precise characterizations of non-computability results for differential equations presented in [65] and [66].
Since there is no general theory known concerning the solvability of all partial differential equations, research in the field focuses on various particular PDEs which are important for applications within and outside of mathematics; for example, the defines a nonlinear map (the solution operator) K R from the Sobolev space H s (R) to the space C(R; H s (R)) for real numbers s ≥ 0, where u represents the amplitude of the wave. It is shown in [94], [93] that for any integer s ≥ 3, the map K R : H s (R) → C(R; H s (R)) is Turing computable, which means that the solution K R (ϕ) can be computed with arbitrary precision on Turing machines when sufficiently good approximations to the initial function ϕ are available (in [30] a preliminary computability result was presented for the periodic case). Note that the solution operator K R is defined for real numbers s ≥ 0 while the computability of K R is established only for integers s ≥ 3 (the proof can be extended to real numbers s ≥ 3 with some modifications). The stronger smoothness requirement is used for ensuring that several derivatives in the construction of approximations remain to be computable. It would be very interesting to know whether K R remains computable for real numbers 0 ≤ s < 3.
We also note that computability results exist for similar PDEs. For example, in [99] the initial and boundary value problem (IBVP) is considered; it is shown that the operator which maps the initial and boundary data to the solution of the IBVP is computable.

The (incompressible) Navier-Stokes Equation
describes the motion of a viscous incompressible fluid filling a rigid box Ω . The vector field u = u(x,t) = u (1) , u (2) , . . . , u (d) represents the velocity of the fluid, P = P(x,t) is the scalar pressure, and f is a given external force. The question of global existence and smoothness of the solutions of Equation (0.8), even in the homogeneous case f ≡ 0, is one of the seven Millennium Prize Problems posted by the Clay Mathematics Institute at the beginning of the 21st century. Local existence has been established, though, in various L p settings [32]; and uniqueness of weak solutions in dimension 2, but not in dimension 3 [80], [14]. Nevertheless, numerical solution methods have been devised in abundance, often based on pointwise (or even uniform, rather than L p ) approximation and struggling with unphysical artifacts [62]. Thus, it becomes useful to know whether it is possible to compute a local solution with arbitrary precision in a rigorous mathematical/physical setting. A recent result provides a positive answer: it is shown in [81] that the solution operator is indeed Turing computable. More precisely, it is shown that there exist a computable map T , T : L σ 2,0 (Ω ) ×C([0, ∞), L σ 2,0 (Ω )) → (0, ∞), (a, f) → T (a, f) and a computable map S , such that for every a ∈ L σ 2,0 (Ω ) and f ∈ C([0, ∞), L σ 2,0 (Ω )), the restriction S (a, f) | [0,T (a,f)] = (u, P) constitutes a (strong local) solution to Equation (0.8), where Ω = (−1, 1) 2 and L σ 2,0 (Ω ) is the set of all square-integrable, divergence-free and boundary-free functions defined on Ω .
The proofs of both theorems make use of energy methods augmented by a fixed point argument by converting the PDE into an integral equation, assembling an iterate scheme, and then proving that the sequence generated by the iteration is computable as well as effectively convergent and the unique limit is the computable solu-tion sought after. There are several typical difficulties in such constructions. Firstly, there is often a need of some custom-designed representations for the spaces of input and output functions, which are different from known canonical ones. For example, the weak solution of the Navier-Stokes equation (0.8) defined on a 2-dimensional good region is a locally square-integrable, solenoidal (i.e. divergence-free), and boundary-free vector field. The canonical representation for an L p space using polynomials with rational coefficients is not rich enough to encode a boundary-free solenoid. Yet a strengthened representation for the boundary-free solenoids must be sufficiently general and robust so that the information carried by the representation won't be destroyed by elementary analytical operations such as (distributive) differentiation and (primitive) integration. The second challenging issue is that the iteration usually involves several operations, which may move back and forth among different spaces. Additional care for representations may become inevitable in order to keep the computations flowing in and out from one space to another. Thirdly, the iterate sequence may (effectively) converge only in a special subspace of the space where the sequence is defined, which may further complicate matters. For instance, Bourgain's space is needed in showing that the iterate sequence is contracting and therefore converges to the solution of the KdV equation.
In addition to particular equations, there are studies on computability aspects for several classes of PDEs. One such class is linear parabolic equations. An initial value problem of a linear parabolic equation can be written in the form: where f : [0, T ] → X is a continuous map and A is a strong elliptic operator, bounded or unbounded, on a Banach space X. When X = L 2 (R n ), the solution to the initial value problem (0.9) is given by the formula: where W (t) is a C 0 semigroup generated by the infinitesimal generator A . Thus if the semigroup W is computable from A , then the solution operator K : (A , f , x) → u is computable uniformly from the data (A , f , x) defining the problem. The study concerning computability of semigroups is carried out in [96], [97]. It is shown in [96] how to compute on Turing machines a uniformly continuous or strongly continuous semigroup uniformly from its infinitesimal generator and vice versa.
Another class for which a master algorithm exists for computing (local) solutions of any of its members is the class of general (nonlinear) first-order PDEs for scalar functions with several independent variables. A problem in this class can be written in the form: where x ∈ U, U is an open subset of R n , F : R n × R ×V → R and g : ∂U → R are given, V is an open subset of R n containing U, the closure of U, and u : U → R is the unknown and Du is the gradient of u. A function u is called a (strong) solution of the boundary value problem if it solves (0.10), (0.11). A special feature of equation (0.10) is that it can, locally at least, be transformed into a system of first-order ordinary differential equations, called the characteristic equations, as follows: to compute the solution u(x) at a point x ∈ U, one finds some curve x(s) lying within U, connecting x with a point x 0 ∈ ∂U, such that along the curve x(s) the boundary value problem for the PDE is reduced to an initial value problem for the characteristic equations. As discussed in Section 2, a system of ODEs is in principle algorithmically uniformly solvable, sometimes quite explicitly. It is shown in [83] that the method of characteristics can be used to compute local solutions of (0.10), (0.11), at least for feasible instances and when the boundary ∂U is effectively locally C 1 . In the same article, it is also shown that although the solution is locally computable the maximal region of existence of a computable local solution may not be computable. Moreover, the problem whether a boundary value problem has a global solution is not algorithmically decidable. The negative results retain even within the class of quasilinear equations defined by analytic computable functions over particularly simple domains (quasilinear equations are among the simplest firstorder nonlinear partial differential equations). This fact shows that the algorithmic unsolvability is intrinsic.
In [53]  In [95] (see also [89], [90]) the computability of the initial value problems for the linear Schrödinger equation u t = i∆ u + φ and the nonlinear Schrödinger equa-tion i u t = −∆ u + mu + |u| 2 u is studied and it is shown that the solution operator is computable if the initial data are Sobolev functions but noncomputable in the linear case if the initial data are L p -functions and p = 2.
In [73] the authors study symmetric hyperbolic systems of PDEs of the form where A, B 1 , . . . , B m are computable constant symmetric n × n matrices, t ≥ 0, and x ∈ [0, 1] m . There it is shown that if the matrices are fixed and if the first and secondorder derivatives of the initial function ϕ are bounded, the operator which maps the initial condition ϕ to the solution of the PDE is computable (see also [74] for a correction and [75,Remark 5.5.3] for the statement given here, which removes some unneeded assumptions used in the original formulation of [73]). Later it was shown in [75,Theorem 5.3] that the operator which maps an initial condition and the matrices A, B 1 , . . . , B m to the solution of the PDE is uniformly computable, provided the matrices satisfy some algebraic conditions. More recently, the computational complexity of solving PDEs with the format (0.12) was studied in [76], [72]. There the (uniform) complexity of this problem is shown to be in EXPTIME in general and in PTIME under certain conditions.
The Dirichlet Problem was also analyzed in [16], [17], [18] using Bishop's framework of constructive mathematics. It is shown that weak solutions of the Dirichlet Problem exist if certain conditions hold.

Some Open Problems
In this chapter, we have presented several results concerning the computability and computational complexity of differential equations. However, the field of differential equations is vast, and many more problems in the field await to be explored from the perspectives of computable analysis. Here we merely give two open problems.
1. The first problem is to give general principles from which the computability or noncomputability of limit sets of ODEs follows as corollaries. For example, it has been shown that a basin of attraction of an ODE defined by a hyperbolic rational function f is computable, uniformly in f ( [9]); however, on the other hand, if f is analytic, computable but non-hyperbolic, a basin of attraction is no longer necessarily computable ( [33]). It would be interesting to know whether a computable, analytic and hyperbolic f ensures the computability of a basin of attraction. 2. The second problem is to characterize the computational complexity of wellknown computable PDEs such as the KdV equation, the Schrödinger equation, and the Navier-Stokes equation.