Bit complexity of computing solutions for symmetric hyperbolic systems of PDEs with guaranteed precision

Abstract

We establish upper bounds on bit complexity of computing solution operators for symmetric hyperbolic systems of PDEs, combining symbolic and approximate algorithms to obtain the solutions with guaranteed prescribed precision. Restricting to algebraic real inputs allows us to use the classical (“discrete”) bit complexity concept.

Keywords

Symmetric hyperbolic system bit complexity solution operator guaranteed precision symbolic computations algebraic real approximation symmetric matrix spectral decomposition difference scheme

1. Introduction

In this paper, which extends the conference paper [35], we establish upper bounds on (classical “discrete”) bit complexity of computing solution operators for symmetric hyperbolic systems of PDEs (restricting the consideration to algebraic real inputs).

Algorithms used in mathematics-oriented software can be divided into two big classes: symbolic algorithms (e.g. [6,10]) which aim to find precise solutions, and approximate algorithms (e.g. [4,17,19,39]) which aim to find “good enough” approximations to precise solutions. The symbolic algorithms are implemented e.g. in computer algebra systems while the approximate algorithms are included into numerical mathematics packages. Both classes of algorithms are widely used in applications and in mathematical research.

Symbolic algorithms correspond well to computations on discrete structures (with mathematical foundations in the classical computability and complexity theory) while the approximate algorithms help to carry out computations on continuous structures (with mathematical foundations in the field of computability and complexity in analysis which however is not yet completely bridged with numerical methods used in practical applications). An important idea relating both classes of algorithms is to look for approximations to the precise solutions with “guaranteed precision”. It is one of the approaches in the fast developing branch of reliable computations which are crucially important for safety-critical applications.

The statement of a guaranteed-precision version of some problem on a continuous structure may reduce it to a problem on a discrete structure which enables to apply the classical computability and bit complexity theory. This is the approach we adopt here to measure the complexity of algorithms for differential equation problems: searching for an approximate solution with guaranteed prescribed precision and using exact symbolic calculations on each step of the algorithm. Classical bit complexity is fundamental because it estimates the amount of computational resources needed to solve a problem on a computing device. We briefly recall the definitions in Section 2.4 (also there are comments on other complexity approaches) and explain how to encode algebraic reals: using encodings of the corresponding minimal (rational) polynomial and the number of the root to which the given algebraic real corresponds, see Section 2.3.

In this paper, we consider initial-value (IVP), or Cauchy, and well posed boundary-value (BVP) problems for symmetric hyperbolic systems of PDEs $A \frac{\partial u}{\partial t} + \sum_{i = 1}^{m} B_{i} \frac{\partial u}{\partial x_{i}} = f (t, x_{1}, \dots, x_{m})$ where $A = A^{*} > 0$ and $B_{i} = B_{i}^{*}$ are symmetric $n \times n$ -matrices, $t ⩾ 0$ , $x = (x_{1}, \dots, x_{m}) \in Q = {[0, 1]}^{m}$ , $f : [0, + \infty) \times Q ⇀ R^{n}$ . Such systems can be used to describe a wide variety of physical processes like those considered in the theories of elasticity, acoustics, electromagnetism etc., see e.g. [13,15,16,23]. Using a well known grid-based method (difference scheme), for algebraic real coefficients and rational polynomials as initial data, we establish an EXPTIME complexity upper bound, and indicate restrictions when it boils down to PTIME: the precision needs to be fixed and certain upper bounds on inputs need to be imposed.

We consider both the IVP and BVP cases as very general and illustrative: they are different in terms of domains of existence and uniqueness as well as properties of the solutions, but almost the same w.r.t. the difference approach. The methods of the paper are applicable to IVP and BVP not only for symmetric hyperbolic systems, but for broader classes of PDE systems.

We stick here to symmetric hyperbolic systems as a quite broad and practically important class capturing one of the following common problems with (explicit) difference scheme for evolutionary systems of PDEs. To construct a stable difference scheme it is necessary to know the number of positive/ negative/ zero eigenvalues and corresponding eigenvectors of matrix pencils related to $A, B_{i}$ , which are known to be non-computable unless e.g. the number of different eigenvalues is given as input [48]. The computational instabilities occuring because of that are well known in practice; here we are carefully avoiding them by restricting to algebraic real inputs, while in [31,34] we add extra restrictions for rigorously proving computability of solution operators in the real setting, see Section 2.1.

Moreover, for the considered case of symmetric hyperbolic systems we provide an explicit PTIME algorithm of computing, from the given precision and input data, the (space and time) grid steps, using the difference scheme with which (or any smaller) provides the solution with this given precision, see Proposition 4 in Section 4.1. A similar approach can be used for broader classes of systems of PDEs, as well as for the real setting [31], which helps to adapt the difference scheme approach to exact real computation (ERC) packages, see [32]. In fact, while the present paper was being reviewed, the first author has published, with two other coauthors, the paper [22] with real complexity estimates for linear evolutionary systems of PDEs (including symmetric hyperbolic systems as a particular case) in the exact real computation approach. The therein obtained estimates help to improve the ones in the present paper: PSPACE instead of EXPTIME, due to replacing the step-by-step difference scheme iterations by a more efficient procedure (recursive matrix powering). However, in [22] we consider fixed real-valued initial functions and matrix coefficients while in the present paper we treat them as inputs, so these papers complement each other; in the difference scheme part [22] partially uses the methods of the present paper; cf. also [27] where the coefficients are fixed as well.

In the next section we formally state the problems to be investigated, and recall some notions and background to be used. In Section 3 we establish upper complexity bounds for several auxiliary linear algebra problems closely related to finding the solutions of symmetric hyperbolic systems of PDEs. In Section 4 we prove the bit cost estimates for the solutions of PDEs, and we conclude in Section 5 with a short discussion on possible applications of this work.

2. Preliminaries and formulations of main results

2.1. Cauchy and boundary-value problems

The Cauchy problem for a symmetric hyperbolic system is stated as follows: $\begin{matrix} (1) & \{\begin{matrix} A \frac{\partial u}{\partial t} + \sum_{i = 1}^{m} B_{i} \frac{\partial u}{\partial x_{i}} = f (t, x_{1}, \dots, x_{m}), & t ⩾ 0, \\ u |_{t = 0} = φ (x_{1}, \dots, x_{m}), \end{matrix} \end{matrix}$ where $A = A^{*} > 0$ and $B_{i} = B_{i}^{*}$ are symmetric $n \times n$ -matrices, $t ⩾ 0$ , $x = (x_{1}, \dots, x_{m}) \in Q = {[0, 1]}^{m}$ , $φ : Q \to R^{n}$ , $f : [0, + \infty) \times Q ⇀ R^{n}$ and $u : [0, + \infty) \times Q ⇀ R^{n}$ is a partial function acting on the domain H of existence and uniqueness of the Cauchy problem (1). The set H is known to be (see e.g. [15]) the intersection of semi-spaces $\begin{matrix} t ⩾ 0, x_{i} - μ_{\max}^{(i)} t ⩾ 0, x_{i} - 1 - μ_{\min}^{(i)} t ⩽ 0 (i = 1, \dots, m) \end{matrix}$ of $R^{m + 1}$ where $μ_{\min}^{(i)}, μ_{\max}^{(i)}$ are respectively the minimum and maximum of the eigenvalues of the matrix $A^{- 1} B_{i}$ .

The boundary-value problem is stated as follows: $\begin{matrix} (2) & \{\begin{matrix} A \frac{\partial u}{\partial t} + \sum_{i = 1}^{m} B_{i} \frac{\partial u}{\partial x_{i}} = f (t, x_{1}, \dots, x_{m}), \\ u |_{t = 0} = φ (x_{1}, \dots, x_{m}), \\ Φ_{i}^{(1)} u (t, x_{1}, \dots, x_{i - 1}, 0, x_{i + 1}, \dots, x_{m}) = 0, \\ Φ_{i}^{(2)} u (t, x_{1}, \dots, x_{i - 1}, 1, x_{i + 1}, \dots, x_{m}) = 0, \\ i = 1, 2, \dots, m, \end{matrix} \end{matrix}$ where $A, B_{i}, Q, φ, f, u$ are as above, and $Φ_{i}^{(1)}$ , $Φ_{i}^{(2)}$ are rectangular matrices (the boundary coefficients) that meet the following conditions:

1) The number of rows of $Φ_{i}^{(1)}$ (respectively, $Φ_{i}^{(2)}$ ) is equal to the number of positive (respectively, negative) eigenvalues of the matrices $A^{- 1} B_{i}$ ; coincidence constraints of the initial and boundary conditions hold (such constraints depend on the particular problem and on the smoothness which we want to obtain).

2) The boundary conditions are assumed to be dissipative which means that $\begin{matrix} (3) & (B_{i} u, u) ⩽ 0 for x_{i} = 0, (B_{i} u, u) ⩾ 0 for x_{i} = 1, i = 1, 2, \dots, m . \end{matrix}$

Condition 1) guarantees the existence of solution of the boundary-value problem (2) in the cylinder $[0, \infty) \times Q$ , while 2) implies its uniqueness.

For both problems (1) and (2), theorems on continuous dependence of solutions on the input data hold, i.e. the problems are correctly posed. Both problems are practically important, since many physical processes (including linear elasticity, acoustic, Maxwell equations) are described by such kind of systems.

For the considered problems, there exist different numerical methods, from which we use those developed in [15,16], see also [23]. Their convergence relies on the following well-known theorem (see e.g. [18,23,37]): if a difference scheme approximates the given differential problem and is stable (which is an intrinsic property of the difference scheme), then the (discrete) solution of the corresponding difference equations converges to the exact solution of the differential problem in an appropriate grid norm; the speed of convergence corresponds to the order of approximation.

The Godunov scheme, which we use, is of first order of approximation and is stable with a Courant number (relating the time and space steps of the grid) depending on spectral characteristics of the matrix coefficients A, $B_{i}$ (see Section 4.1 below). As it is known, in order to construct a stable difference scheme for a symmetric hyperbolic system, one needs to compute eigenvectors of symmetric matrices, which is actually a discontinuous operation [29], hence not computable. However, due to [47,48], eigenvectors are computable, provided that the spectrum cardinality (i.e., the number of distinct eigenvalues of the given matrix) is given as an input.

In [31 ,33 ,34] we developed an approach to study computability properties of PDEs based on the Godunov difference scheme and established the computability, in the sense of the TTE approach [7,42], of the solution operator of the Cauchy and boundary-value problems. In this section we briefly recall main properties of such systems and the related computability results.

For computability of the solution of (1) from the initial data and the matrix coefficients, the following result has been established.

Theorem 1 ([31]).

Let $M_{φ} > 0, M_{A} > 0, p ⩾ 2$ be integers, let $i = 1, \dots, m$ , and let $n_{A}, n_{1}, \dots, n_{m}$ be cardinalities of spectra of A and of the matrix pencils $λ A - B_{1}, \dots, λ A - B_{m}$ , respectively (i.e., $n_{i}$ is the number of distinct roots of the characteristic polynomial $det (λ A - B_{i})$ ). Then the operator $\begin{matrix} R : (A, B_{1}, \dots, B_{m}, n_{A}, n_{1}, \dots, n_{m}, φ) \mapsto u \end{matrix}$ sending any sequence $A, B_{1}, \dots, B_{m}$ of symmetric real matrices with $A > 0$ such that the matrix pencils $λ A - B_{i}$ have no zero eigenvalues, $\begin{matrix} (4) & ‖ A ‖_{2}, {‖ A^{- 1} ‖}_{2}, ‖ B_{i} ‖_{2} ⩽ M_{A}, λ_{min}^{(i)} < 0 < λ_{max}^{(i)}, i = 1, 2, \dots, m, \end{matrix}$ the sequence $n_{A}, n_{1}, \dots, n_{m}$ of the corresponding cardinalities, and any function $φ \in C^{p + 1} (Q, R^{n})$ satisfying the conditions $\begin{matrix} (5) & {‖ \frac{\partial φ}{\partial x_{i}} ‖}_{s} ⩽ M_{φ}, {‖ \frac{\partial^{2} φ}{\partial x_{i} \partial x_{j}} ‖}_{s} ⩽ M_{φ}, i, j = 1, 2, \dots, m, \end{matrix}$ to the unique solution $u \in C^{p} (H, R^{n})$ of ( 1 ) (with $f = 0$ ) is a computable partial function from the space $S_{+} \times S^{m} \times N^{m + 1} \times C_{s}^{p + 1} (Q, R^{n})$ to $C_{s L_{2}}^{p} (H, R^{n})$ .

Here $C_{s}^{p + 1} (Q, R^{n})$ , $C_{s L_{2}}^{p} (H, R^{n})$ are the spaces of continuously differentiable functions, being in $C^{p + 1}$ or $C^{p}$ on the corresponding sets, such that all of their first and second derivatives are uniformly bounded by $M_{φ}$ , with the sup-norm $\begin{matrix} ‖ φ ‖_{s} = sup_{x \in Q} | φ (x) | \end{matrix}$ on $C^{p + 1} (Q, R^{n})$ and $s L_{2}$ -norm $\begin{matrix} ‖ u ‖_{s L_{2}} = sup_{{t : (t, x) \in H}} \sqrt{\int_{Q} {| u (x, t) |}^{2} d x} \end{matrix}$ on $C^{p} (H, R^{n})$ , respectively; $| u (x, t) |^{2} = ⟨ u, u ⟩$ is the standard scalar product.

By S and $S^{+}$ we denote respectively the spaces of all symmetric and symmetric positively definite matrices with euclidean norms uniformly bounded by a constant $M_{A}$ .

In [31], an analogue of Theorem 1 is established for the boundary-value problem (2) with fixed computable real matrices $Φ_{i}^{(1)}$ , $Φ_{i}^{(2)}$ $(i = 1, 2, \dots, m)$ , provided that the strong dissipativity property (with strict inequalities in (3)) holds. The solution for the boundary-value problem is defined on a cylinder $[0, T] \times Q$ , for any computable $T > 0$ , i.e. $u \in C_{s L_{2}}^{p} ([0, T] \times Q, R^{n})$ .

Remark 1.
The computability of the operator $R : (A, B_{1}, \dots, B_{m}, φ) \mapsto u$ without the additional information on spectra cardinalities is by now an open question in general, though in particular cases like Cauchy or periodic boundary-value problems computability can be established by different methods like e.g. Fourier transformation. In case of non-computability, it is instructive to study the topological complexity of computing the operator. For the case when the number of different eigenvalues of a matrix is given as input, a PTIME algorithm of computing eigenvectors has been recently devised by Sewon Park (the paper soon to appear); also eigenvalues are known to be PTIME computable [44], see also [28]. Combined with proofs of our Theorems 1, 5 these facts give an analog of Theorem 5 in the real complexity approach [21]; those estimates can then be further improved as in [22]. For different approaches to define complexity see Remark 2.

Interestingly, if we work in decidable (in the Russian terminology, strongly constructivizable [11]) fields (for example, the field of algebraic reals $A$ used in the present paper), we do not need to add cardinalities as inputs.
Theorem 2 ([31]).

Let $(B, β)$ be a decidable (i.e., strongly constructive) real closed ordered subfield of $R$ . Then the solution operator $R : (A, B_{1}, \dots, B_{m}, φ) \mapsto u$ is uniformly computable (w.r.t. the numbering β) on matrices $A, B_{1}, \dots, B_{m}$ with coefficients in $B$ .

For the boundary-value problem (2) an analogue of this result is proved, with uniformity on the matrices $Φ_{i}^{(1)}$ , $Φ_{i}^{(2)}$ having coefficients in $B$ .

In [34], dependence on the right-hand part f in (1) is added, and also the time $T > 0$ is included into the arguments of the solution operator R, i. e. the solution belongs to $C_{s L_{2}}^{p} (H \cap [0, T] \times Q, R^{n})$ . In [34] the technique of using decidable (i.e., strongly constructive) fields for the computability results is explained in detail.

Also, Theorem 2 implies computability for fixed computable real matrices (which we earlier established using the result of [48]).

Theorem 3 ([31,33]).

Let $Q = {[0, 1]}^{m}$ ; $T > 0$ be a computable real and $M_{φ} > 0$ , $p ⩾ 2$ be integers. Let $A, B_{1}, \dots, B_{m}$ be fixed computable symmetric matrices, such that $A = A^{*} > 0$ , $B_{i} = B_{i}^{*}$ .

If $φ \in C^{p + 1} (Q)$ satisfies the conditions ( 5 ), then the operator $R : φ \mapsto u$ mapping the initial function φ to the unique solution $u \in C^{p} (H, R^{n})$ of the Cauchy problem ( 1 ) is a computable partial function from $C_{s}^{p + 1} (Q, R^{n})$ to $C_{s L_{2}}^{p} (H, R^{n})$ .

An analogue of Theorem 3 is established in [31] for the boundary-value problem (2) with fixed computable real matrices $Φ_{i}^{(1)}$ , $Φ_{i}^{(2)}$ $(i = 1, 2, \dots, m)$ .

The next natural step is to study computational complexity of the mentioned problems, which we start to do in the present paper.

2.2. Discretization of the problems

To investigate the complexity of computing solutions of the systems (1) and (2) (and even to formulate the results), we need discrete approximations of the given and unknown functions, as well as their interpolations. Therefore we first recall some discretization details.

Consider, for any positive integer N, the uniform rectangular grid $G_{N}$ on $Q = {[0, 1]}^{m}$ defined by the points $\begin{matrix} (\frac{i_{1} - \frac{1}{2}}{2^{N}}, \frac{i_{2} - \frac{1}{2}}{2^{N}}, \dots, \frac{i_{m} - \frac{1}{2}}{2^{N}}) \end{matrix}$ where $1 ⩽ i_{1}, i_{2}, \dots, i_{m} ⩽ 2^{N}$ . Let $h = 1 / 2^{N}$ be the corresponding spatial grid step and τ be a time step. Denote $G_{N}^{τ} = G_{N} \times {l τ}_{l = 1}^{L}$ , where L is the number of the time steps. The choice of steps h and τ, which guarantee good properties of the difference scheme, will be specified below in Section 4.1.

Note that the number of points in the grid $G_{N}$ is $2^{N m}$ , so the set $Q^{G_{N}}$ (or $A^{G_{N}}$ ) of grid functions $g^{(h)} : G_{N} \to Q^{n}$ (resp. $g^{(h)} : G_{N} \to A^{n}$ ) may be identified with $Q^{n \cdot 2^{N m}}$ (resp. $A^{n \cdot 2^{N m}}$ ). We will consider the following grid norms $\begin{matrix} {‖ g^{(h)} ‖}_{s} = {max}_{x \in G_{N}} | g^{(h)} (x) |, {‖ g^{(h)} ‖}_{L_{2}}^{2} = h^{m} \sum_{x \in G_{N}} ⟨ g^{(h)} (x), g^{(h)} (x) ⟩ . \end{matrix}$

We will consider the $s L_{2}$ -norm on the vector spaces $Q^{G_{N}^{τ}}$ (or $A^{G_{N}^{τ}}$ ) of grid functions $v^{(h)} (t, x)$ on such grids: $\begin{matrix} {‖ v^{(h)} ‖}_{s L_{2}} = {max}_{t \in {l τ}_{l = 1}^{M}} h^{m} \sum_{x \in G_{N}} ⟨ v^{(h)} (t, x), v^{(h)} (t, x) ⟩ . \end{matrix}$

Recall that multilinear interpolations $\tilde{u}$ (linear on each coordinate and coinciding with u at the grid points) have the following properties: $u \mapsto \tilde{u}$ and $u^{(h)} \mapsto \tilde{u^{(h)}}$ are linear, and the following estimate holds [38,45]: $\begin{matrix} (6) & ‖ u - \tilde{u} ‖_{s} ⩽ max_{i, j} {{‖ \frac{\partial^{2} u}{\partial x_{i} \partial x_{j}} ‖}_{s L_{2}}, {‖ \frac{\partial^{2} u}{\partial x_{i} \partial t} ‖}_{s L_{2}}} \cdot h^{2} . \end{matrix}$

Further in Section 4 we will construct, by means of a stable difference scheme approximating the differential system (1), a grid function υ on $G_{N}^{τ}$ such that $\begin{matrix} (7) & ‖ u - \tilde{υ ∣_{H}} ‖_{s L_{2}} < \frac{1}{a}, \end{matrix}$ where $a > 1$ is a given integer (which determines the precision $\frac{1}{a}$ of computation) and u is the solution of (1).

Recall that, for an abstract boundary-value problem $\begin{matrix} (8) & \{\begin{matrix} L u (y) = f (y) \in C^{p} (Ω, R^{n}), & y \in Ω \subset R^{k} \\ L u (y) |_{Γ} = φ (y ∣_{Γ}) \in C^{q} (Γ, R^{n}), & Γ \subseteq \partial Ω . \end{matrix} \end{matrix}$ (where L and $L$ are differential operators with the differential order of $L$ less than that of L, Γ is a part of the boundary $\partial Ω$ of some area Ω), a difference scheme is a system of algebraic equations $\begin{matrix} (9) & L_{h} u^{(h)} = f^{(h)}, L_{h} u^{(h)} = φ^{(h)} . \end{matrix}$

Here $L_{h}, L_{h}$ are difference operators (which are in our case linear); $u^{(h)}$ and $φ^{(h)}$ are grid functions.

The scheme (9) approximates the differential equations (8) with order of accuracy l (where l is a positive integer) on a solution $u (t, x)$ of (8) if $\begin{array}{l} {‖ (L u) |_{G_{k}} - L_{h} u^{(h)} ‖}_{F_{h}} ⩽ M_{1} h^{l}, {‖ f |_{G_{k}} - f^{(h)} ‖}_{F_{h}} ⩽ M_{2} h^{l}, \\ {‖ (L u) |_{G_{k}} - L_{h} u^{(h)} ‖}_{Φ_{h}} ⩽ M_{3} h^{l} and {‖ φ |_{G_{k}} - φ^{(h)} ‖}_{Φ_{h}} ⩽ M_{4} h^{l} \end{array}$ for some constants $M_{1}, M_{2}, M_{3}$ and $M_{4}$ not depending on h and τ.

The difference scheme (9) is called stable if its solution $u^{(h)}$ satisfies $\begin{matrix} {‖ u^{(h)} ‖}_{U_{h}} ⩽ N_{1} {‖ f^{(h)} ‖}_{F_{h}} + N_{2} {‖ φ^{(h)} ‖}_{Φ_{h}} \end{matrix}$ for some constants $N_{1}$ and $N_{2}$ not depending on h, τ, $f^{(h)}$ and $φ^{(h)}$ .

Fact 1 ([18]).

Let the difference scheme be stable and approximate (1) on the solution u with order l. Then the solution $u^{(h)}$ uniformly converges to the solution u in the sense that $‖ u |_{G_{k}^{τ}} - u^{(h)} ‖_{U_{h}} ⩽ N h^{l}$ for some constant N not depending on h and τ.

The difference scheme which we use for approximating (1) will be described in detail in Section 4.

2.3. Algebraic preliminaries

In the study of computability of solution operators we considered rather general classes of initial data $A, B_{i}, φ_{i}, \dots$ (matrices with real coefficients, broad classes of smooth functions and so on). In contrast, the study of complexity suggests to consider first more restricted classes of objects admitting fast enough computations. Here we briefly (and not very systematically) recall some relevant algebraic notions and facts. For details see e.g. [40].

Coefficients of matrices and polynomials will usually be taken from a fixed (ordered) field $F \subseteq R$ of reals. For a field $F$ , let $F [x_{1}, x_{2}, \dots]$ denote the ring of polynomials with coefficients in $F$ and variables $x_{1}, x_{2}, \dots$ , and let $F (x_{1}, x_{2}, \dots)$ be the corresponding field of fractions (i.e., the field of rational functions with variables $x_{1}, x_{2}, \dots$ ). A polynomial $p \in F [x]$ is written as $p = a_{0} + a_{1} x + \dots + a_{n} x^{n}$ where $a_{n} \neq 0$ for $p \neq 0$ ; n is the degree of p denoted as $deg (p)$ ; if $a_{n} = 1$ , p is called unitary. A polynomial $p \in F [x]$ is irreducible (over $F$ ) if it is not the product of two polynomials from $F [x]$ of lesser degrees.

The arithmetics of polynomials is very similar to the arithmetics of the integers $Z$ , in particular any non-zero polynomial p, $deg (p) ⩾ 2$ , has a canonical (i.e., unique up to permutation of factors) factorisation $p = a p_{1}^{m_{1}} \dots p_{k}^{m_{k}}$ where $a \in F$ , $k, m_{1}, \dots, m_{k} ⩾ 1$ and $p_{1}, \dots, p_{k}$ , $deg (p_{i}) ⩾ 1$ , are unitary irreducible polynomials such that $(p_{i}, p_{j}) = 1$ for $i \neq j$ ( $(p_{i}, p_{j})$ is the greatest common divisor of $p_{i}, p_{j}$ ). It is known that if $0 \neq α \in C$ is algebraic over $F$ (i.e. α is a root of some polynomial over $F$ ) then the smallest subfield $F (α)$ of $C$ containing $F \cup {α}$ is isomorphic to the quotient field $F [x] / (p_{α})$ where $p_{α}$ is the unique unitary polynomial of minimal degree $p \in F [x]$ with $p (α) = 0$ and $(p_{α})$ is the ideal of $F [x]$ generated by $p_{α}$ ; $deg (α) = deg (p_{α})$ is called the degree of α over $F$ .

In this paper we most often work with ordered fields $F \in {Q, Q (α), A ∣ 0 \neq α \in A}$ where $Q$ is the ordered field of rationals, $A$ is the ordered field of algebraic reals (which consists of the reals algebraic over $Q$ ). One can also consider the smallest subfield $Q (α_{1}, \dots, α_{n})$ of $R$ containing given $α_{1}, \dots, α_{n} \in A$ . For any such field there is a “primitive element” $α \in A$ with $Q (α) = Q (α_{1}, \dots, α_{n})$ . Note that $A$ is the smallest real closed ordered field. With any non-zero $α \in A$ we associate the unique pair $(p_{α}, k)$ such that k satisfies $α = α_{k}$ where $α_{1} < \dots < α_{m}$ is the increasing sequence of all real roots of the polynomial $p_{α}$ .

Let $M_{n} (R)$ be the set of $n \times n$ -matrices over a (commutative associative with a unit element 1) ring R, and $M (R)$ be the union of all $M_{n} (R)$ , $n ⩾ 1$ . We use standard terminology and notation from linear algebra. In particular, $det (A)$ is the determinant of $A = (a_{i j}) \in M_{n} (R)$ , $diag (a_{1}, \dots, a_{n})$ is the diagonal matrix with the diagonal elements $a_{1}, \dots, a_{n} \in R$ , so in particular $I = I_{n} = diag (1, \dots, 1)$ is the unit matrix. The roots of the polynomial $c h_{A} = det (λ I - A)$ are called eigenvalues of $A \in M_{n} (C)$ . In general, the eigenvalues of a real matrix are complex numbers. The eigenvalues of a symmetric real matrix are always real.

2.4. Encodings and bit complexity

Computations on existing computers (as well as on theoretical computing devices like Turing machines) don’t work with abstract mathematical objects (like integers, rationals or polynomials) but rather with words over a finite alphabet. Here we briefly recall some relevant notions and facts (for more details see e.g. [2,3,5,9,25,30]).

Let Σ be a finite alphabet, $Σ^{*}$ the set of words over Σ, $S \subseteq {(Σ^{*})}^{n}$ , and $t : S \to ω$ . A function $f : S : \to Σ^{*}$ is computable in time t if there exists a k-tape Turing machine, $k ⩾ n + 1$ , such that starting to work with words $x_{1}, \dots, x_{n} = \bar{x} \in S$ written on the first n tapes, finishes within $t (\bar{x})$ steps with $f (\bar{x})$ written on the $(n + 1)$ -st tape. The set S is computable in time t if so is its characteristic function $χ_{S} : {(Σ^{*})}^{n} \to {0, 1}$ (assuming that $0, 1 \in Σ$ ). A structure $S = (S; \dots)$ of a finite signature is computable in time t if so are its universe S and all the signature functions and relations. An abstract structure $A$ is t-time-presentable if it is isomorphic to a structure $B$ computable in time t. Any isomorphism from $A$ onto $B$ is a t-time-presentation of $A$ . Note that usually people work not directly with Turing machines but rather with informal algorithms on words; the algorithms should use only elementary enough steps to make it clear how to translate them to the syntax of Turing machines.

If the function $t (\bar{x})$ is bounded by some polynomial on $| \bar{x} | = max {| x_{1} |, \dots, | x_{k} |}$ (where $| x |$ is the length of a word x) we say that the function f (resp., the set S) is computable in polynomial time (p-computable for short). An abstract structure $A$ is p-presentable if it is isomorphic to a p-computable structure $B$ . Any isomorphism from $A$ onto $B$ will be called a p-presentation of $A$ . Similarly one can define computability in linear, quadratic or exponential time.

In order to enable Turing machines work with abstract objects, we have to encode the objects by words in a finite alphabet. First note (see e.g. [5]) that for any finite alphabet Σ there is a natural encoding (injective function) $c : Σ^{*} \to {0, 1}^{*}$ of words over Σ by binary words such that the function c, its range $r n g (c)$ , and the inverse function $c^{- 1}$ are computable in linear time, so in all interesting cases we can without loss of generality stick to binary encodings (this is the reason why the classical computational complexity is often called bit complexity). We give several examples of such encodings for different sets A of abstract objects. In describing such an encoding $b : A \to {0, 1}^{*}$ we often use the trick of first describing an auxiliary encoding $e : A \to Σ^{*}$ for some bigger alphabet Σ and then setting $b = c \circ e$ . Below we often take $Σ = {0, 1, *}$ where ∗ is a new symbol used for separating binary words.

A positive integer n is usually encoded by its binary notation $b (n)$ , so $| b (n) | = log (n)$ . Adding an additional bit for the sign, we obtain a binary encoding b of the integers. Identifying rationals with fractions $\frac{p}{q}$ where $p, q \in Z$ , $q ⩾ 1$ and $(p, q) = 1$ , we can define an encoding $c : Q \to {0, 1, *}^{*}$ by $c (\frac{p}{q}) = b (p) * b (q)$ (and, by the mentioned trick with alphabets, we can modify c to obtain a binary coding b of the rationals). The defined encodings give p-presentations of the ordered ring $Z$ and the ordered field $Q$ .

The binary encoding b of $Q$ induces the encoding $e : Q [x] \to {0, 1, *}^{*}$ which associates with a non-zero polynomial $p = a_{0} + a_{1} x + \dots + a_{n} x^{n}$ , $a_{n} \neq 0$ , the code $b (a_{0}) * \dots * b (a_{n})$ . In a similar way one can define natural induced encodings of $Q [x_{1}, \dots, x_{n}]$ and of $Q (x_{1}, \dots, k_{n})$ for each $n ⩾ 1$ which provide p-presentations of the corresponding rings and fields (see e.g. [3] for additionl details). Moreover, in the field $Q (x_{1}, \dots, x_{n})$ also the subtraction and division, as well as the evaluation function $Q [x_{1}, \dots, x_{n}] \times Q^{n} \to Q$ are p-computable. Furthermore, these encodings induce encodings of $Q [x_{1}, \dots] = ⋃_{n} Q [x_{1}, \dots, x_{n}]$ and of $Q (x_{1}, \dots)$ modulo which any of $Q [x_{1}, \dots, x_{n}], Q (x_{1}, \dots, x_{n})$ is p-computable; moreover, the evaluation partial function $Q [x_{1}, \dots] \times Q^{*} \to Q$ (and the similar function for $Q (x_{1}, \dots)$ ), where $Q^{*}$ is the set of finite strings over $Q$ , is p-computable.

If $p \in Q [x]$ is unitary irreducible then the quotient-field $Q [x] / (p)$ (formed by the polynomials of degree $< deg (p)$ ) also has a natural p-presentation (see e.g. [1,3]). Therefore, the ordered field $Q (α)$ , for each non-zero $α \in A$ , has a natural p-presentation. For the induced binary presentation of $F = Q (α)$ we can repeat the constructions of the previous paragraph and obtain natural p-presentations of $F [x_{1}, \dots]$ and of $F (x_{1}, \dots)$ . Again, the polynomial evaluation $F [x_{1}, \dots] \times F^{*} \to F$ will be p-computable w.r.t. the natural induced encodings.

We also define a natural encoding of $A$ into ${0, 1, *}^{*}$ by associating with any non-zero $α \in A$ the word $b (p_{α}) * b (k)$ where $(p_{α}, k)$ is the pair from the previous subsection and $b : Q [x] \to {0, 1}^{*}$ is the natural binary coding specified above. From deep results of computer algebra (including the polynomial time algorithm for factoring rational polynomials [24]) it follows that in this way we obtain a natural p-presentation of the ordered field $A$ in which also subtraction and division are p-computable (see e.g. [3,25]).

Simplifying notation, we always denote the natural binary p-presentations of any of the ordered fields $F \in {Q, Q (α), A ∣ 0 \neq α \in A}$ by b. Note however that some important computational properties of the presentation of $A$ discussed above differ from those for the presentations of $Q$ and $Q (α)$ . In particular, the evaluation function $A [x_{1}, \dots] \times A^{*} \to A$ is now computable in exponential time but not in polynomial time. In fact, already the “long sum” operation $(α_{1} * \dots * α_{n}) \mapsto α_{1} + \dots + α_{n}$ is not PTIME-computable uniformly on n (even not computable in PSPACE) w.r.t. the presentation of $A$ . This follows from the results in [46] (for a detailed explanations see comments after the proof of Theorem 4 in [3]).

We will use some results from [3] about the complexity of root-finding in the field $C_{alg} = (C_{alg}; +, \times, 0, 1)$ of complex algebraic numbers, i.e. of finding all roots of an equation $α_{e} x^{e} + \dots + α_{1} x + α_{0} = 0$ where $α_{i} \in C_{alg}$ for $i ⩽ e$ . More precisely, the authors of [3] consider equations of the form $\begin{matrix} (10) & t_{e} (α_{1}, \dots, α_{k}) x^{e} + \dots + t_{1} (α_{1}, \dots, α_{k}) x + t_{0} (α_{1}, \dots, α_{k}) = 0, \end{matrix}$ where $α_{1}, \dots, α_{k} \in C_{alg}$ and $t_{j} (\bar{x}) \in Q [x_{1}, \dots, x_{k}]$ . The problem is to find a list of (codes of) all roots of (10) from given $b (α_{1}) * \dots * b (α_{k})$ and $b (t_{0} (\bar{x})) * \dots * b (t_{e} (\bar{x}))$ where b is a natural binary encoding of $C_{alg}$ induced by the presentation of $A$ described above and by Gauss representation of complex numbers as pairs of reals. As shown in Theorem 8 [3], the problem is solvable in polynomial time for any fixed k. Moreover, the same estimate holds for the version of this problem when one computes the list of all distinct real roots of (10) in increasing order.

The introduced binary presentations of fields, polynomials and rational functions are natural in the sense that they are p-equivalent to some presentations really used in computer algebra systems. See e.g. [3,25] for additional details.

Associate with any matrix $A = (a_{i j}) \in M_{n} (F)$ its code $c (A) \in {0, 1, *}^{*}$ by $c (A) = b (a_{11}) * \dots * b (a_{1 n}) * b (a_{21}) * \dots * b (a_{2 n}) * \dots * b (a_{n 1}) * \dots * b (a_{n n})$ . These encodings induce the binary encoding of the set $M (F) = ⋃_{n} M_{n} (F)$ of all square matrices over $F$ in which any set $M_{n} (F)$ is p-computable. Furthermore, many matrix properties like symmetricity are also p-computable.

It is well known (see e.g. [30]) that our encodings give p-presentations of the rings $M_{n} (Q)$ uniformly on n (uniformity means that there is a polynomial bound working for all n). Moreover, it is easy to check that evaluation of some “long” terms in these rings are also p-computable uniformly on n w.r.t. these presentations (in particular the function $A_{1} * \dots * A_{n} \mapsto A_{1} \times \dots \times A_{n}$ is p-computable). Even more involved matrix algorithms like computing of the determinant, computing of the inverse of a non-degenerate matrix, and Gauss method also work in polynomial time (see Chapter 3 [30] for additional details).

From results in [3,25] (using also the arguments in [30]) it follows that the polynomial time estimates of the previous paragraph remain true for the rings $M_{n} (Q (α))$ for each non-zero $α \in A$ . Moreover, our presentation of the ring $M_{n} (A)$ is a p-presentation for any fixed $n ⩾ 1$ , but not uniformly on n. The “long terms” are not p-computable w.r.t. our presentation for $M_{n} (A)$ , even for a fixed n.

The introduced encodings of matrices are natural in the sense they are closely related to standard encodings of matrices in numerical analysis. The only difference is that here we use a precise symbolic encoding of matrix coefficients while in numerical analysis the floating-point approximations of coefficients are usually used.

Remark 2.
Note that, along with bit complexity, there are other approaches to measuring complexity of computations on, say, Euclidean spaces, including algebraic complexity (see e.g. [8]), topological complexity (see e.g. [36,41]), real complexity, i.e. bit complexity adapted to reals by means of computable analysis (see e.g. [7,20,42,43]). Algebraic complexity counts the arithmetical operations needed for finding a symbolic solution (the notion we use is close to that, but also taking into account the bit costs of performing the operations). Topological complexity counts the number of equality and order tests needed to find a solution of a discontinuous (hence, non-computable in the sense of [7,42]) problem. Real complexity counts steps of a Turing machine working on the representations of reals (when restricted to algebraic reals, it is close to our approach).

2.5. Formulations of main results

Now we have enough notions and terminology to state the guaranteed-precision problems in a rigorous form.

First we consider the task of computing the domain H of existence and uniqueness of the Cauchy problem. As mentioned in Section 2.1, the set H is the intersection of semi-spaces $\begin{matrix} t ⩾ 0, x_{i} - μ_{\max}^{(i)} t ⩾ 0, x_{i} - 1 - μ_{\min}^{(i)} t ⩽ 0, (i = 1, \dots, m) \end{matrix}$ of $R^{m + 1}$ where $μ_{\min}^{(i)}, μ_{\max}^{(i)}$ are respectively the minimum and maximum of the eigenvalues of the matrix $A^{- 1} B_{i}$ . Therefore, the computation of H reduces to the computation of the eigenvalues of the matrices $A^{- 1} B_{i}$ .

Our algorithms for solving the Cauchy problem will be for technical reasons presented only for the case when H satisfies the condition $μ_{\min}^{(i)} < 0 < μ_{\max}^{(i)}$ for all $i = 1, \dots, m$ ; this condition often holds for natural physical systems. Note that this condition implies that H is a compact subset of $[0, + \infty) \times Q$ .

In [33] we observed that the domain H for the problem (1) is computable from matrices $A, B_{1}, \dots, B_{m}$ (more exactly, the vector $(μ_{max}^{(1)}, \dots, μ_{max}^{(m)}, μ_{min}^{(1)}, \dots, μ_{min}^{(m)})$ is computable from $A, B_{1}, \dots, B_{m}$ ; this implies the computability of H in the sense of computable analysis [42]).

The next result establishes the complexity of computing H satisfying the mentioned condition.

Theorem 4.
Let $m, n ⩾ 2$ be any fixed integers. There is a polynomial time algorithm which for any given $A, B_{1} \dots, B_{m} \in M_{n} (A)$ finds the vector $(μ_{max}^{(1)}, \dots, μ_{max}^{(m)}, μ_{min}^{(1)}, \dots, μ_{min}^{(m)})$ and checks the condition $μ_{\min}^{(i)} < 0 < μ_{\max}^{(i)}$ for all $i = 1, \dots, m$ . Thus, the algorithm finds the domain H satisfying the condition above, or reports on the absence of such a domain.

Now we state two guaranteed-precision versions of a restricted (to algebraic matrix coefficients and rational polynomials as initial data) Cauchy problem (1). Let $m, n$ be fixed positive integers.

Task 1. $CP (m, n, A, Q)$ is the following computational task:

INPUT: Integer $a ⩾ 1$ , polynomials $φ_{1} \dots, φ_{n} \in Q [x_{1} \dots, x_{m}], f_{1} \dots, f_{n} \in Q [t, x_{1} \dots, x_{m}]$ and matrices $A, B_{1} \dots, B_{m} \in M_{n} (A)$ .

OUTPUT: (Codes of) a rational $T > 0$ with $H \subseteq [0, T] \times Q$ , a spatial rational grid step h dividing 1, a time grid step τ dividing T and a rational $h, τ$ -grid function $v : G_{N}^{τ} \to Q$ such that $\begin{matrix} (11) & ‖ u - \tilde{υ ∣_{H}} ‖_{s L_{2}} < \frac{1}{a} . \end{matrix}$

This problem is a standard computational task on a discrete structure which asks for an algorithm (and its complexity estimation) which, for any given input computes a suitable output. The second version of the Cauchy problem is formulated as follows.

Task 2. $CP (m, n, a, M, A, Q)$ is the following computational task:

INPUT: (Codes of) polynomials $\begin{matrix} φ_{1} \dots, φ_{n} \in Q [x_{1} \dots, x_{m}], f_{1} \dots, f_{n} \in Q [t, x_{1} \dots, x_{m}] \end{matrix}$ and matrices $A, B_{1} \dots, B_{m} \in M_{n} (A)$ such that the quantities $‖ A ‖_{2}$ , $\frac{λ_{\max} (A)}{λ_{\min} (A)}$ , $\begin{matrix} max_{i} {‖ B_{i} ‖_{2}, {‖ {(A^{- 1} B_{i})}^{2} ‖}_{2}, max_{k} {| μ_{k} | : det (μ_{k} A - B_{i}) = 0}, sup_{t, x} {‖ \frac{\partial^{2} f}{\partial x_{i} \partial t} (t, x) ‖}_{2}}, \end{matrix}$ and $\begin{matrix} max_{i, j} {{‖ A^{- 1} B_{i} A^{- 1} B_{j} - A^{- 1} B_{j} A^{- 1} B_{i} ‖}_{2}, sup_{t, x} {‖ \frac{\partial^{2} f}{\partial x_{i} \partial x_{j}} (t, x) ‖}_{2}, sup_{x} {‖ \frac{\partial^{2} φ}{\partial x_{i} \partial x_{j}} (x) ‖}_{2}} \end{matrix}$ are bounded by M.

OUTPUT: (Codes of)a rational $T > 0$ with $H \subseteq [0, T] \times Q$ , a spatial rational grid step h dividing 1, a time grid step τ dividing T and a rational $h, τ$ -grid function $v : G_{N}^{τ} \to Q$ such that $‖ u - \tilde{υ ∣_{H}} ‖_{s L_{2}} < \frac{1}{a}$ .

The guaranteed-precision versions $BV (m, n, a, M, A, Q)$ and $BV (m, n, A, Q)$ of boundary-value problem (2) are stated in a similar way, with an additional assumption ${max}_{i} {‖ Φ_{i}^{(1)} ‖_{2}, ‖ Φ_{i}^{(2)} ‖_{2}} ⩽ M$ in the second version.

Our basic result, which concerns Tasks 1 and 2, is formulated as follows (the definitions of the well known complexity classes PTIME and EXPTIME may be found e.g. in [5]).
Theorem 5.

For any $m, n ⩾ 1$ , the problems $CP (m, n, A, Q)$ and $BVP (m, n, A, Q)$ are solvable in EXPTIME.

For any $m, n, a, M ⩾ 1$ , the problems $CP (m, n, a, M, A, Q)$ and $BVP (m, n, a, M, A, Q)$ are solvable in PTIME.

One could consider several variations of the stated problems. For instance, one could vary some of the remaining fixed parameters or take rational functions $φ_{1} \dots, φ_{n} \in Q (x_{1} \dots, x_{m}), f_{1} \dots, f_{n} \in Q (t, x_{1} \dots, x_{m})$ instead of polynomials. We will make some comments on the first option below. For the second option, analogues of Theorem 5 may be proved using the methods of this paper. In our approach, we cannot take arbitrary algebraic polynomials as the initial and right-hand side functions because computing their values, by the above-mentioned result in [46], requires exponential time which yields worse upper bounds than in Theorem 5 (see Section 4.3). Remark 3.
From the proof in the last section of this paper it will follow that similar upper bounds hold for the algebraic complexity.

3. Auxiliary algorithms

In this section we present some linear algebra algorithms which are used in Section 4 to prove the main results.

3.1. Computing spectral decomposition

Here we give upper bounds for the complexity of symbolic computations of eigenvalues and eigenvectors for some classes of matrices and matrix pencils. The computations are w.r.t. the encodings specified in Section 2.4.

By spectral decomposition of a symmetric real matrix $A \in M_{n} (R)$ we mean a pair $((λ_{1}, \dots, λ_{n}), (v_{1}, \dots, v_{n}))$ where $λ_{1} ⩽ \dots ⩽ λ_{n}$ is the non-decreasing sequence of all eigenvalues of A (each eigenvalue occurs in the sequence several times, according to its multiplicity) and $v_{1}, \dots, v_{n}$ is a corresponding orthonormal basis of eigenvectors.

Proposition 1.

For any fixed $n ⩾ 1$ , there is a polynomial time algorithm which, given a symmetric matrix $A \in M_{n} (A)$ , computes a spectral decomposition of A.

There is an algorithm which, given a symmetric matrix $A \in M (Q)$ , computes (uniformly on n) the spectrum $(λ_{1}, \dots, λ_{n})$ of A in PTIME and a spectral decomposition of A in EXPTIME. 1
¹
Formulation of this item is weaker than in [35] where the estimate for orthonormalized eigenvectors is incorrect because the Gram–Schmidt process does not work in PTIME uniformly on n.

Proof.
1. For fixed n, the coefficients of $c h_{A} = det (λ I - A) \in A [λ]$ are fixed integer polynomials from the elements of A, so we can compute them in polynomial time using Theorem 4 in [3]. By results in [25] and Theorem 8 in [3] cited in Section 2.3, we can compute in polynomial time the increasing sequence $μ_{1} < \dots < μ_{m}$ of all roots of $c h_{A}$ and the corresponding multiplicities $r_{1}, \dots, r_{m}$ (hence we can also compute the sequence $λ_{1}, \dots, λ_{n}$ ). Compute in polynomial time a primitive element $β \in A$ with $Q (β) = Q (α, μ_{1}, \dots, μ_{m})$ . Working in $Q (β)$ , we find in polynomial time (in the usual way, solving corresponding linear systems by the Gauss method), for each $j = 1, \dots, m$ , a basis $w_{1}^{i}, \dots, w_{r_{j}}^{i}$ for the eigenspace of $μ_{j}$ . Applying the Gram–Schmidt orthogonalisation process (which also works in polynomial time for fixed n) and normalizing the obtained orthogonal basis, we obtain an orthonormal basis for this eigenspace. Putting together the orthonormal bases for all j, we obtain a desired orthornormal basis $(v_{1}, \dots, v_{n})$ for the whole space.

2. Given $A \in M (Q)$ , compute first the order n of A. Next we compute the characteristic polynomial $c h_{A} = det (λ I - A) \in Q [λ]$ , written in the form $c h_{A} = λ^{n} - p_{1} λ^{n - 1} - p_{2} λ^{n - 2} - \dots - p_{n}$ . By remarks at the end of Section 2.4, we can compute in polynomial time the traces $s_{1}, \dots, s_{n}$ of matrices $A^{1}, \dots, A^{n}$ respectively. It is known (see e.g. Section 4.4 of [14]) that the equalities $p_{1} = s_{1}$ , and $(k + 1) p_{k + 1} = s_{k + 1} - p_{k} s_{k}$ for $k = 1, \dots, n - 1$ , hold. From these we subsequently compute $p_{1}, p_{2}, \dots, p_{n}$ in polynomial time.

Since factorisation in $Q [λ]$ is computable in polynomial time [24], we can compute the canonical decomposition $c h_{A} = q_{1}^{k_{1}} \dots q_{l}^{k_{l}}$ where $l, k_{1}, \dots, k_{l}$ are positive integers. Let $λ_{j, 1} < \dots < λ_{j, d_{j}}$ be the sequence of all roots of $q_{j}$ where $d_{j} = deg (q_{j})$ . Then $c (q_{j}, i)$ is the code of $λ_{j, i}$ for each $i = 1, \dots, d_{j}$ . Since ${λ_{j, i} ∣ j = 1, \dots, l, i = 1, \dots, d_{j}} = {λ_{1}, \dots, λ_{n}}$ and all $λ_{j, i}$ are pairwise distinct, we have computed the (codes of the) eigenvalues $λ_{1}, \dots, λ_{n}$ . Since the multiplicity of $λ_{j, i}$ is $k_{j}$ , we have also computed the multiplicity $r_{i}$ of any eigenvalue $λ_{i}$ (the multiplicity of $λ_{j, i}$ is $k_{j}$ ).

Since the eigenvalues $λ_{j, i}$ are pairwise distinct, it remains to find an orthonormal basis of the eigenspace ${x ∣ (λ_{j}^{i} \cdot I - A) x = 0)}$ corresponding to any fixed $λ_{j, i}$ (the dimension of this space is $k_{j}$ ). Applying the Gauss method (in the field of coefficients $Q (λ_{j, i})$ ) we find a desired basis $w_{j, 1}, \dots, w_{j, k_{j}}$ (note that, since the field $Q (λ_{j, i})$ is isomorphic to $Q [λ] / (p_{j})$ for each $i = 1, \dots, d_{j}$ , the systems may be solved uniformly on i in the quotient field, with subsequent evaluation of the computed polynomials at $λ_{j, i}$ ). By remarks in Section 2.4, this computation runs in polynomial time. Applying the Gram–Schmidt orthogonalization process and normalizing the obtained orthogonal bases for each $i = 1, \dots, d_{j}$ , we put them together to get a resulted orthonormal basis for the whole space in EXPTIME uniformly on n, using Theorems 4 and 8 in [3]. □
Remark 4.
In fact, Theorem 8 in [3] implies that a spectral decomposition of any given symmetric algebraic matrix may be computed in exponential time uniformly on n, though in this case we cannot compute in polynomial time even coefficients of $c h_{A}$ , again because of the result in [46].

By matrix pencil we mean a pair $(A, B)$ (often written in the form $μ A - B$ ) of real non-degenerate symmetric matrices such that A is positive definite (i.e., all of its eigenvalues are positive). By spectral decomposition of such a pencil we mean a tuple $\begin{matrix} ((λ_{1}, \dots, λ_{n}), (v_{1}, \dots, v_{n}), (μ_{1}, \dots, μ_{n}), (w_{1}, \dots, w_{n})) \end{matrix}$ such that $((λ_{1}, \dots, λ_{n}), (v_{1}, \dots, v_{n}))$ and $((μ_{1}, \dots, μ_{n}), (w_{1}, \dots, w_{n}))$ are spectral decompositions of the symmetric matrices A and $D^{} L^{} B L D$ respectively, where L is the matrix formed by vectors $v_{1}, \dots, v_{n}$ written as columns and $D = diag {\frac{1}{\sqrt{λ_{1}}}, \frac{1}{\sqrt{λ_{2}}}, \dots, \frac{1}{\sqrt{λ_{n}}}}$ . Proposition 2.
For any fixed $n ⩾ 1$ , there is a polynomial time algorithm which, given a matrix pencil $(A, B)$ with $A, B \in M_{n} (A)$ , computes a spectral decomposition of $(A, B)$ .
Proof.
By item 1 of the previous theorem, we can find in polynomial time a spectral decomposition $((λ_{1}, \dots, λ_{n}), (v_{1}, \dots, v_{n}))$ of A. Since we can solve polynomial equations in $Q (λ_{1}, \dots, λ_{n})$ in polynomial time (see e.g. [3,25] for details), we can compute in polynomial time the matrix $D^{} L^{} B L D$ . Applying item 1 of the previous theorem to this matrix (and working now in the field $Q (λ_{1}, \dots, λ_{n})$ which is also computable in polynomial time [25]), we compute the remaining items $(μ_{1}, \dots, μ_{n}), (w_{1}, \dots, w_{n})$ . Note that $(μ_{1}, \dots, μ_{n})$ coincides with the sequence of eigenvalues of (in general, non-symmetric) matrix $A^{- 1} B$ . □

3.2. Computing data for the difference scheme

Here we explain how to compute data needed for computations with the difference schemes in Section 4.

Let $A, B_{1}, \dots, B_{m} \in M_{n} (A)$ be matrices satisfying the conditions in Cauchy problem. We can compute the spectral decomposition $((λ_{1}, \dots, λ_{n}), (v_{1}, \dots, v_{n}))$ of A as in the proof of Proposition 1. Let $λ_{\max}$ , $λ_{\min}$ be respectively the maximum and minimum of $λ_{1}, \dots, λ_{n}$ . Let L be the orthonormal matrix formed by vectors $v_{1}, \dots, v_{n}$ written in columns, so $L^{*} A L = Λ = diag {λ_{1}, λ_{2}, \dots, λ_{n}}$ , and let $D = Λ^{- \frac{1}{2}}$ .

For each $i = 1, \dots, m$ , let $((μ_{1}^{(i)}, \dots, μ_{n}^{(i)}), (w_{1}^{i}, \dots, w_{n}^{i}))$ be the spectral decomposition of the symmetric matrix $D^{*} L^{*} B_{i} L D$ computed as in the proof of Proposition 2. Let $μ_{\max}^{(i)}$ , $μ_{\min}^{(i)}$ be respectively the maximum and minimum of $μ_{1}^{(i)}, \dots, μ_{n}^{(i)}$ . Let $M_{i} = diag {μ_{1}^{(i)}, \dots, μ_{n}^{(i)})}$ and $K_{i}$ be the orthonormal matrix formed by vectors $w_{1}^{i}, \dots, w_{n}^{i}$ written in columns, so $K_{i}^{*} D^{*} L^{*} B_{i} L D K_{i} = M_{i}$ . Let $T_{i} = L D K_{i}$ for each $i = 1, \dots, m$ .

From Propositions 1, 2 and remarks in Section 2.4 we easily obtain:

Proposition 3.
For any fixed $m, n ⩾ 1$ , there is a polynomial time algorithm which, given matrices $A, B_{1}, \dots, B_{m}, \in M_{n} (Q)$ satisfying the conditions of symmetric hyperbolic systems, computes the objects $A^{- 1}$ , $T_{i}$ , $T_{i}^{- 1}, λ_{\max}$ , $λ_{\min}, μ_{\max}^{(i)}, μ_{\min}^{(i)}$ , $μ_{k}^{(i)}$ ( $i = 1, \dots, m, k = 1, \dots, n$ ) specified above.

Note that in the proof of Theorem 5 in Section 4 we will stick (for notational simplicity) to the typical particular case $m = 2$ where the notations $T_{x} = T_{1}, T_{y} = T_{2}$ and $K_{x} = K_{1}, K_{y} = K_{2}$ are more appealing, e.g. for considering the linear transformations of variables.
4. Proof of the main results

Theorem 4 straightforwardly follows from Proposition 2. We give a proof of Theorem 5 for the Cauchy problem and the boundary-value problem simultaneously, with a minor modification of the numerical algorithm used in our proof.

4.1. Computing the grid steps

In the Proposition below we assume that all the assumptions of Section 2.5 hold. This proposition is the main technical tool of the present paper.

Proposition 4.
Let $m, n ⩾ 1$ be fixed integers. The time and space steps τ and h guaranteeing ( 11 ), when using the difference scheme ( 19 ) on the corresponding grid, are PTIME computable from $a, φ, f, A, B_{i}$ ( $i = 1, 2, \dots, m$ ).
Proof.
W.l.o.g., stick to the case $m = 2$ , denoting $B_{1} = B$ , $B_{2} = C$ . First note [16] that the Godunov scheme (see its detailed description below in the next subsection) is stable if and only if $\begin{matrix} (12) & τ ⩽ h \cdot {(\frac{1}{{max}_{i} {| μ_{i} | : det (μ_{i} A - B) = 0}} + \frac{1}{{max}_{i} {| μ_{i} | : det (μ_{i} A - C) = 0}})}^{- 1} \end{matrix}$ (see also a short summary of the proof of this fact in [33]). Recall that stability is an intrinsic property of a difference scheme, implying, together with approximation, its convergence to the corresponding differential equation in grid norms.

By Proposition 2, a rational τ satisfying (12) can be found in polynomial time, if h is computed in polynomial time. It is also obvious that τ can be chosen so that $L = \frac{T}{τ}$ is integer. Thus it suffices to estimate complexity of finding h.

In [31] we established, applying well known theorems about difference schemes and interpolations, that $\begin{matrix} (13) & ‖ u - \tilde{υ ∣_{H}} ‖_{s L_{2}} ⩽ ‖ u - \tilde{u ∣_{G_{N}^{τ}}} ‖_{s L_{2}} + ‖ \tilde{u ∣_{G_{N}^{τ}}} - \tilde{υ ∣_{H}} ‖_{s L_{2}} ⩽ c_{int} h + c_{diff} h ⩽ C h ⩽ \frac{1}{a}, \end{matrix}$ where the “constants” $c_{int}$ and $c_{diff}$ depend only on $A, B_{1}, \dots, B_{m}$ and the first and second partial derivatives of φ. So we can take $h = \frac{1}{2^{N}}$ such that $\begin{matrix} (14) & h ⩽ \frac{1}{a C}; \end{matrix}$ it remains to estimate the complexity of finding an upper bound $C$ on $(c_{int} + c_{diff})$ from the input data.

We claim that $\begin{matrix} (15) & max {c_{int}, c_{diff}} ⩽ P (A, B, C, φ), \end{matrix}$ where $\begin{array}{l} P (A, B, C, φ) = & \frac{λ_{\max} (A)}{λ_{\min} (A)} \cdot max {{‖ \frac{\partial^{2} φ}{\partial x_{i} \partial x_{j}} ‖}_{s}} \\ \times max {‖ A ‖_{2}, ‖ B ‖_{2}, ‖ C ‖_{2}, {‖ {(A^{- 1} B)}^{2} ‖}_{2}, {‖ {(A^{- 1} C)}^{2} ‖}_{2}, {‖ A^{- 1} B A^{- 1} C - A^{- 1} C A^{- 1} B ‖}_{2}} . \end{array}$

As noted in [18] (see chapter 5), $c_{diff} = c_{1} \cdot c_{2}$ where $c_{2} = \sqrt{\frac{λ_{\max} (A)}{λ_{\min} (A)}}$ comes from the stability property [16] and $c_{1}$ comes from the approximation property $‖ L_{h} u_{h} - (L u) |_{G_{k}^{τ}} ‖_{s L_{2}} ⩽ c_{1} h$ . Since our scheme has the first order of approximation, it follows from the Taylor decomposition of $L u$ that $c_{1}$ depends only on the norms of $A, B_{i}$ and $\begin{matrix} {‖ \frac{\partial u}{\partial x_{i}} ‖}_{s L_{2}}, {‖ \frac{\partial u}{\partial t} ‖}_{s L_{2}}, {‖ \frac{\partial^{2} u}{\partial x_{i} \partial x_{j}} ‖}_{s L_{2}}, {‖ \frac{\partial^{2} u}{\partial x_{i} \partial t} ‖}_{s L_{2}} . \end{matrix}$

As it is known, by the proof of the uniqueness theorem for (1) [15] (p. 155 for the Cauchy problem and p. 194 for the boundary-value problem, respectively), see also [12,26], we have $\begin{matrix} (16) & ‖ u ‖_{A, s L_{2}} ⩽ ‖ φ ‖_{A, L_{2}} \end{matrix}$ Applying an analogue of (16) to the systems for the second derivatives of u and using the equivalence of norms in $R^{n}$ $\begin{matrix} λ_{\min} (v, v) ⩽ ‖ v ‖_{A} = (A v, v) ⩽ λ_{\max} (v, v), \end{matrix}$ we obtain the desired estimate. More precisely, the estimates for the derivatives of the solution can be obtained as follows. Considering the Cauchy problem, due to the smoothness assumptions, we can construct auxiliary Cauchy problems for partial derivatives of u (we write down a couple of them, as examples): $\begin{array}{l} (17) & \begin{aligned} \{\begin{matrix} A {(u_{x})}_{t} + B {(u_{x})}_{x} + C {(u_{x})}_{y} = 0, \\ u_{x} |_{t = 0} = φ_{x}, \end{matrix} \{\begin{matrix} A {(u_{t})}_{t} + B {(u_{t})}_{x} + C {(u_{t})}_{y} = 0, \\ u_{t} |_{t = 0} = - A^{- 1} (B φ_{x} + C φ_{y}), \end{matrix} \\ \{\begin{matrix} A {(u_{t t})}_{t} + B {(u_{t t})}_{x} + C {(u_{t t})}_{y} = 0, \\ u_{t t} |_{t = 0} = - A^{- 1} (B {(u_{t} |_{t = 0})}_{x} + C {(u_{t} |_{t = 0})}_{y}) \\ = {(A^{- 1} B)}^{2} φ_{x x} + {(A^{- 1} C)}^{2} φ_{y y} + (A^{- 1} B A^{- 1} C - A^{- 1} C A^{- 1} B) φ_{x y} . \end{matrix} \end{aligned} \end{array}$ From (6) $c_{int}$ is easily estimated in a similar way as above (by bounding the derivatives of u).

From these considerations, $h ⩽ \frac{1}{a C} ⩽ \frac{1}{2 a P (A, B, C, φ)}$ can be computed in PTIME, since all the expressions in $P (A, B, C, φ)$ (eigenvalues, matrix multiplication, taking an inverse matrix, calculating upper bounds on the norms, differentiating rational polynomials) are PTIME computable. □

Note that for not fixed n, the PTIME bounds in Proposition 4 fail to hold, since there are no PTIME bounds on finding eigenvalues, see Remark 4.
4.2. Description of the algorithm

The difference scheme for the boundary-value problem (2) and the Cauchy problem (1) may be chosen in various ways. We use the Godunov scheme [16] (see also e.g. [23]), which can also be applied to a broader class of systems, including some systems of nonlinear equations. We describe it in few stages, letting for simplicity the righthand part of (1) to be zero: $f = 0$ . The scheme approximates the system (1) or (2) with the first order of accuracy (the proof of the approximation property is done by means of the Taylor decomposition). To simplify the analysis of complexity, we write the scheme in an algorithmic form.

The algorithm describes, in several steps, computation of the values $u^{i - \frac{1}{2}, j - \frac{1}{2}}$ (approximating the solution of (1) or (2) at the point $(t, \frac{i - \frac{1}{2}}{2^{N}}, \frac{j - \frac{1}{2}}{2^{N}})$ ) on the time level $t = (l + 1) τ$ from the values $u_{i - \frac{1}{2}, j - \frac{1}{2}}$ on the time level $t = l τ$ . We distinguish the values by upper and lower indices, in order to avoid using of a third index l.

0. Find the matrices $T_{x}$ , $T_{y}$ , $T_{x}^{- 1}$ , $T_{y}^{- 1}$ , $A^{- 1}$ , the eigenvalues $μ_{k} (B)$ , $μ_{k} (C)$ of the matrix pencils $μ A - B$ , $μ A - C$ ( $k = 1, 2, \dots, n$ ) and the steps $h, τ$ from given $m, n, a$ , algebraic (resp. rational) matrices $A, B, C$ and the initial function φ.

1. Let (for i from 1 to $2^{N}$ , j from 1 to $2^{N}$ ) $u_{i - \frac{1}{2}, j - \frac{1}{2}} = φ ∣_{G_{N}}$ , i.e. on the level t=0 we just take the values of initial conditions $φ (\frac{i - \frac{1}{2}}{2^{N}}, \frac{j - \frac{1}{2}}{2^{N}})$ .

The next steps 2-5 are carried out for all l from 1 to $L = \frac{T}{τ}$ .

2. For i from 1 to $2^{N}$ , j from 1 to $2^{N}$ , find the auxiliary vectors $\begin{array}{l} v_{i - \frac{1}{2}, j - \frac{1}{2}}^{(x)} = T_{x}^{- 1} u_{i - \frac{1}{2}, j - \frac{1}{2}}, v_{i + \frac{1}{2}, j - \frac{1}{2}}^{(x)} = T_{x}^{- 1} u_{i + \frac{1}{2}, j - \frac{1}{2}} \\ v_{i - \frac{1}{2}, j - \frac{1}{2}}^{(y)} = T_{y}^{- 1} u_{i - \frac{1}{2}, j - \frac{1}{2}}, v_{i - \frac{1}{2}, j + \frac{1}{2}}^{(y)} = T_{y}^{- 1} u_{i - \frac{1}{2}, j + \frac{1}{2}} . \end{array}$

3. Find the auxiliary “large values” $\begin{matrix} V_{i, j - \frac{1}{2}}^{(x)} = (\begin{matrix} W_{i, j - \frac{1}{2}}^{1} \\ W_{i, j - \frac{1}{2}}^{2} \\ \dots \\ W_{i, j - \frac{1}{2}}^{n} \end{matrix}); V_{i - \frac{1}{2}, j}^{(y)} = (\begin{matrix} W_{i - \frac{1}{2}, j}^{1} \\ W_{i - \frac{1}{2}, j}^{2} \\ \dots \\ W_{i - \frac{1}{2}, j}^{n} \end{matrix}) \end{matrix}$ for $i = 0, 1, 2, \dots, 2^{N}$ , $j = 1, 2, \dots, 2^{N}$ ; $i = 1, 2, \dots, 2^{N}$ , $j = 0, 1, 2, \dots, 2^{N}$ .

For auxiliary “interior” points i from 1 to $2^{N} - 1$ , j from 1 to $2^{N} - 1$ let $\begin{matrix} (18) & W_{i, j - \frac{1}{2}}^{k} = \{\begin{matrix} {(v_{i - \frac{1}{2}, j - \frac{1}{2}}^{(x)})}_{k}, & if μ_{k} (B) ⩾ 0, \\ {(v_{i + \frac{1}{2}, j - \frac{1}{2}}^{(x)})}_{k}, & if μ_{k} (B) < 0; \end{matrix} W_{i - \frac{1}{2}, j}^{k} = \{\begin{matrix} {(v_{i - \frac{1}{2}, j - \frac{1}{2}}^{(y)})}_{k}, & if μ_{k} (C) ⩾ 0, \\ {(v_{i - \frac{1}{2}, j + \frac{1}{2}}^{(y)})}_{k}, & if μ_{k} (C) < 0, \end{matrix} \end{matrix}$ where $μ_{k} (B), μ_{k} (C)$ are respectively the k-th eigenvalues of the matrix pencils $μ A - B$ , $μ A - C$ .

3a). In the case of the Cauchy problem (1), for the auxiliary “boundary” values $W_{0, j - \frac{1}{2}}$ , $W_{2^{N}, j - \frac{1}{2}}$ (for j from 1 to $2^{N}$ ) and $W_{i - \frac{1}{2}, 0}$ , $W_{i - \frac{1}{2}, 2^{N}}$ (for i from 1 to $2^{N}$ ) we use the same formula (18) taking $\begin{matrix} v_{- \frac{1}{2}, j - \frac{1}{2}}^{(x)} = v_{2^{N} - \frac{1}{2}, j - \frac{1}{2}}^{(x)}, v_{2^{N} + \frac{1}{2}, j - \frac{1}{2}}^{(x)} = v_{\frac{1}{2}, j - \frac{1}{2}}^{(x)} \end{matrix}$ and $\begin{matrix} v_{i - \frac{1}{2}, - \frac{1}{2}}^{(y)} = v_{i - \frac{1}{2}, 2^{N} - \frac{1}{2}}^{(y)}, v_{i - \frac{1}{2}, 2^{N} + \frac{1}{2}}^{(y)} = v_{i - \frac{1}{2}, \frac{1}{2}}^{(y)} . \end{matrix}$

3b). In the case of the boundary-value problem (2), we compute the boundary values $V_{0}, V_{2^{N}}$ with the help of the boundary conditions. On the left boundary $x = 0$ we calculate $m_{+}$ components of $V_{0}$ , corresponding to the positive eigenvalues of the matrix $A^{- 1} B$ , from the system of linear equations $Φ_{1}^{(1)} V_{0} = 0$ ; for $m_{-}$ components of $V_{0}$ , corresponding to the negative eigenvalues, we let $V_{0} : = v_{\frac{1}{2}}$ . The components corresponding to the zero eigenvalues of $A^{- 1} B$ can be chosen arbitrarily since they are multiplied by zero in the scheme. The values on the right boundary and on both boundaries by the y-coordinate are calculated in a similar way.

4. Calculate $U_{i, j - \frac{1}{2}} = T_{x} V_{i, j - \frac{1}{2}}^{(x)}$ (for i from 0 to $2^{N}$ , for j from 1 to $2^{N}$ ) and $U_{i - \frac{1}{2}, j} = T_{y} V_{i - \frac{1}{2}, j}^{(y)}$ (for i from 1 to $2^{N}$ , for j from 0 to $2^{N}$ ).

5. Find values on the next grid step: for i from 1 to $2^{N}$ , j from 1 to $2^{N}$ let $\begin{matrix} (19) & u^{i - \frac{1}{2}, j - \frac{1}{2}} = u_{i - \frac{1}{2}, j - \frac{1}{2}} - \frac{τ}{h} A^{- 1} (B (U_{i, j - \frac{1}{2}} - U_{i - 1, j - \frac{1}{2}}) + C (U_{i - \frac{1}{2}, j} - U_{i - \frac{1}{2}, j - 1})) . \end{matrix}$

Remember the calculated values as $u_{i - \frac{1}{2}, j - \frac{1}{2}}^{(l)} = u^{i - \frac{1}{2}, j - \frac{1}{2}}$ , then let $u_{i - \frac{1}{2}, j - \frac{1}{2}} = u^{i - \frac{1}{2}, j - \frac{1}{2}}$ .

6. Finally, $υ = {u_{i - \frac{1}{2}, j - \frac{1}{2}}^{(l)}}_{l = 1}^{L} ∣_{H}$ is the approximation of the solution u of the system (1), and $υ = {u_{i - \frac{1}{2}, j - \frac{1}{2}}^{(l)}}_{l = 1}^{L} ∣_{{[0, 1]}^{m} \times [0, T]}$ is the approximation of the solution u of the system (2).

Remark 5.
Stage 3 of the algorithm can be carried out without using the branching operator, by letting $\begin{matrix} W_{i, j - \frac{1}{2}} = S_{-}^{B} v_{i + \frac{1}{2}, j - \frac{1}{2}}^{(x)} + S_{+}^{B} v_{i - \frac{1}{2}, j - \frac{1}{2}}^{(x)}; W_{i - \frac{1}{2}, j} = S_{-}^{C} v_{i - \frac{1}{2}, j + \frac{1}{2}}^{(y)} + S_{+}^{C} v_{i - \frac{1}{2}, j - \frac{1}{2}}^{(y)} . \end{matrix}$ Here $S_{-}^{B}$ (resp. $S_{-}^{C}$ ) is the matrix $diag {1, 1, \dots, 1, 0, 0, \dots, 0}$ , with the number of 1s equal to the number of negative eigenvalues of $A^{- 1} B$ (resp. $A^{- 1} C$ ); similarly, $S_{+}^{B}$ (resp. $S_{+}^{C}$ ) is the matrix $diag {0, 0, \dots, 0, 1, 1, \dots, 1}$ with the number of 1s equal to the number of nonnegative eigenvalues of $A^{- 1} B$ (resp. $A^{- 1} C$ ). Note that the matrices $S_{-}^{B}$ , $S_{-}^{C}$ , $S_{+}^{B}$ , $S_{+}^{C}$ can be computed on Stage 0 before the cycles, and that these matrices depend only on the signs of the eigenvalues $μ_{k} (B)$ , $μ_{k} (C)$ .

4.3. Counting steps

With all this at hand, it is not hard to count the computation steps in the Godunov scheme (all computations are w.r.t. the p-presentation of $A$ in Section 2.3). By Proposition 4, the number of grid points in the scheme (see Section 2.2) is bounded by a polynomial. The computations in Godunov’s scheme proceed bottom-up by layers, along the time axis. At the bottom level, we just evaluate the initial functions in the grid points which requires polynomial time according to remarks in Section 2.3. To go one level up requires, for each grid point on the next level, the values at the previous levels and a fixed number of matrix multiplications by matrices, computed in advance using Propositions 1 and 2. Therefore, climbing one level up also requires polynomial time. Let $p_{i}$ , $i = 1, \dots, L$ (where $L = \frac{T}{τ}$ , computed in polynomial time), be a polynomial bounding the computation time for level i. Since the computation at level i uses only the values of υ at grid points of level $i - 1$ (note that computation of the value at any point from the i-th level requires only finite number of points at the ( $i - 1$ )th level, in our case five $2 m + 1$ points) and some matrices computed in advance, the whole computation time is (essentially) bounded by the composition $p_{L} \circ \dots \circ p_{1}$ of polynomials which lays down to EXPTIME and yields the complexity bound given in item 1 of Theorem 5.

For item 2 the argument is the same, except that one has to take into account that for a fixed n the algorithm of spectral decomposition works in polynomial time and that the number of time steps L is just constant (instead of PTIME computable) in this case. It follows from the estimates (12), (14) and (15) of Proposition 4, in particular the expression $P (A, B, C, φ)$ used to calculate h can be taken just as $M^{3}$ . □

Remark 6.
1) Note that the estimate in Theorem 5 (1) is exponential even for fixed n because we take arbitrary rational polynomials as initial functions. Taking reasonably restricted classes of initial functions and matrices (with restrictions like those in [31,33]) yields a polynomial estimate in item (2), though the degree of polynomial is high and one needs to take a fine grid with small steps h and τ. In this way, the exponential algorithm of item (1) might work out better than the polynomial one of item (2) for concrete problems. Also note that in item (2) the precision is fixed.

2) For the case of simultaneously diagonalizable (or, equivalently, mutually commuting) matrices $A^{- 1} B_{j}$ , the solution of the IVP (1) even in the sense of Task 1 is in PTIME (without using the difference scheme method). Indeed, for $n = 1$ , the scalar transport equation $u_{t} = \sum_{j = 1}^{m} b_{j} u_{x_{j}}$ can be easily solved in PTIME: $u (t, x_{1}, \dots, x_{m}) = φ (x_{1} - b_{1} t, \dots, x_{m} - b_{m} t)$ . For a fixed $n > 1$ , the system (1) can be linearly transformed to n independent transport equations via spectral decomposition of the matrix pencils $(A, B_{j})$ , which is computable in PTIME according to Proposition 2.

5. Conclusion

In this paper we obtained apparently first bit complexity upper bounds for computing solutions of the Cauchy and dissipative boundary-value problems for symmetric hyperbolic systems of PDEs (to which also many higher-order hyperbolic PDEs can be reduced) with guaranteed precision.

Although our methods do not always yield (for instance, for large n) practically feasible algorithms for guaranteed-precision problems for PDEs, we hope that investigations in this direction are fruitful for both theoretical research and applications. In particular, on the implementation level it seems useful and rewarding to enhance the existing systems of “exact real computations” by packages based on highly developed algorithms of computer algebra. We are not aware of the existence of such “hybrid” systems built under the slogan of “guaranteed-precision numerical computations”.

Footnotes

Acknowledgements

We are grateful to Sergey Goncharov for stimulating discussions, and to Pavel Alaev, Gennadiy Demidenko and Martin Ziegler for valuable comments.

This work was partially supported by RFBR-JSPS Grant 20-51-50001. The first author was partially supported by the National Research Foundation of Korea (grant 2017R1E1A1A03071032), by the NRF Brain Pool program (grant 2019H1D3A2A02102240), by the International Research & Development Program of the Korean Ministry of Science and ICT (grant 2016K1A3A7A03950702) and the Russian Foundation of Fundamental Research (grant 17-01-00801). The second author was partially supported by the Regional Mathematical Center of Kazan Federal University.

References

A.G.

Akritas, Elements of Computer Algebra with Applications, Wiley Interscience, New York, 1989.

P.E.

Alaev, Existence and uniqueness of structures computable in polynomial time, Algebra and Logic 55(1) (2016), 106–112. doi:10.1007/s10469-016-9377-6.

P.E.

Alaev and

V.L.

Selivanov, Fields of algebraic numbers computable in polynomial time. I, Algebra and Logic 58(6) (2019), 673–705. doi:10.33048/alglog.2019.58.601.

K.I.

Babenko, in: Foundations of Numerical Analysis, Moscow, Nauka, 1986, (in Russian).

J.L.

Balcázar,

Díaz and

Gabarró, Structural complexity I, in: EATCS Monographs on Theoretical Computer Science, Vol. 11, Springer-Verlag, 1988.

Basu,

Pollack and

Roy, Algorithms in Real Algebraic Geometry, Springer, Heidelberg, 2006.

Brattka,

Hertling and

Weihrauch, A tutorial on computable analysis, in: New Computational Paradigms,

S.B.

Cooper,

Löwe and

Sorbi, eds, 2008, pp. 425–491.

Bürgisser,

Clausen and

Shokrollahi, Algebraic Complexity Theory, Springer, Berlin-Heidelberg, 1997.

Cenzer and

Remmel, Polynomial time versus recursive models, Annals of Pure and Applied Logic 54 (1991), 17–58. doi:10.1016/0168-0072(91)90008-A.

10.

G.E.

Collins and

Loos, Real zeros of polynomials, in: Computer Algebra: Symbolic and Algebraic Computations, Springer-Verlag, 1982, pp. 83–94. doi:10.1007/978-3-7091-3406-1_7.

11.

Y.L.

Ershov and

S.S.

Goncharov, Constructive Models, Scientific Book, Novosibirsk, 1999, (in Russian, there is an English Translation).

12.

L.C.

Evans, Partial Differential Equations, Graduate Studies in Mathematics, Vol. 19, American Mathematical Society, 1998.

13.

K.O.

Friedrichs, Symmetric hyperbolic linear differential equations, Communication on Pure and Applied Mathematics 7 (1954), 345–392. doi:10.1002/cpa.3160070206.

14.

F.R.

Gantmacher, Matrix Theory, Nauka, Moscow, 1967, (in Russian).

15.

S.K.

Godunov, Equations of Mathematical Physics, Nauka, Moscow, 1971, (in Russian).

16.

S.K.

Godunov (ed.), Numerical Solution of Higher-Dimensional Problems of Gas Dynamics, Nauka, Moscow, 1976, (in Russian).

17.

S.K.

Godunov

et al., in: Guaranteed Precision of Solving Systems of Linear Equations in Euclidean Spaces, Novosibirsk, Nauka, 1988, (in Russian).

18.

S.K.

Godunov and

V.S.

Ryaben’kii, Introduction to the Theory of Difference Schemes, Fizmatgiz, Moscow, 1962, (in Russian). English translation: Difference Schemes: An Introduction to the Underlying Theory (Studies in Mathematics and Its Applications), Elsevier Science Ltd, 1987.

19.

John, Lectures on Advanced Numerical Analysis, Gordon and Breach, Science Publishers, Inc., 1966.

20.

Kawamura and

Ziegler, Invitation to Real Complexity Theory: Algorithmic Foundations to Reliable Numerics with Bit-Costs, arXiv:1801.07108.

21.

Ker-I, Complexity Theory of Real Functions, Birkhäuser, Boston, 1991.

22.

Koswara,

Selivanova and

Ziegler, Computational complexity of real powering and improved solving linear differential equations, in: Proc. 14th InternationalComputer Science Symposium in Russia, LNCS, Vol. 11532, 2019.

23.

A.G.

Kulikovskii,

N.V.

Pogorelov and

A.Y.

Semënov, Mathematical Aspects of Numerical Solution of Hyperbolic Systems, Chapman & Hall/CRC Press, Boca Raton, 2001.

24.

A.K.

Lenstra,

H.W.

Lenstra and

Lovasz, Factoring polynomials with rational coefficients, Math. Ann. 261 (1982), 515–534. doi:10.1007/BF01457454.

25.

Loos, Computing in algebraic extensions, in: Computer Algebra: Symbolic and Algebraic Computations, Springer-Verlag, 1982, pp. 115–138. doi:10.1007/978-3-7091-3406-1_9.

26.

Mizohata, The Theory of Partial Differential Equations, Cambridge Univ. Press, Cambridge, 1973.

27.

Pan and

Reif, The bit complexity of discrete solutions of partial differential equations: Compact multigrid, Computers Math. Applic. 20(2) (1990), 9–16. doi:10.1016/0898-1221(90)90235-C.

28.

V.Y.

Pan and

Z.Q.

Chen, The complexity of the matrix eigenproblem, in: Proceedings of the Thirty-First Annual ACM Symposium on Theory of Computing, ACM, 1999, pp. 507–516. doi:10.1145/301250.301389.

29.

Rellich, Störungstheorie der Spektralzerlegung I., Analytische Stor̈ung der isolierten Punkteigenwerte eines beschränkten Operators, Math. Ann. 113 (1937), 600–619, (in German). doi:10.1007/BF01571652.

30.

Schrijver, Theory of Linear and Integer Programming, Wiley Interscience, New York, 1986.

31.

Selivanova and

Selivanov, Computing solution operators of boundary-value problems for some linear hyperbolic systems of PDEs, Logical Methods in Computer Science 4(13) (2017), 1–31, Earlier version on

arXiv:1305.2494

(2013).

32.

Selivanova and

Ziegler, Turnkey solutions to PDEs in exact real computation, in: Book of Abstracts – 12th Summer Workshop on Interval Methods, Palaiseau, France, July 23–26, 2019, pp. 21–22, https://swim2019.ensta-paris.fr/SWIM-19_Book_of_Abstracts.pdf .

33.

S.V.

Selivanova and

V.L.

Selivanov, Computing solution operators of symmetric hyperbolic systems of PDEs, Journal of Universal Computer Science 15(6) (2009), 1337–1364.

34.

S.V.

Selivanova and

V.L.

Selivanov, On constructive number fields and computability of solutions of PDEs, Doklady Mathematics 477(3) (2017), 282–285.

35.

S.V.

Selivanova and

V.L.

Selivanov, Bit complexity of computing solutions for symmetric hyperbolic systems of PDEs (extended abstract), in: Proceedings of the Conference Computability in Europe,

Manea,

Miller and

Novotka, eds, LNCS, Vol. 10936, Springer, Berlin, 2018.

36.

Smale, On the topology of algorithms, J. of Complexity 3 (1987), 81–89. doi:10.1016/0885-064X(87)90021-5.

37.

J.C.

Strikverda, Finite Difference Schemes and Partial Differential Equations, SIAM, 2004.

38.

Szegö, Orthogonal Polynomials, AMS 531, West 116th Street New York, 1959.

39.

L.N.

Trefethen, Finite Difference and Spectral Methods for Ordinary and Partial Differential Equations, Cornell University, Department of Computer Science and Center for Applied Mathematics, 1996.

40.

B.L.

van der Waerden, Algebra, Springer, Berlin, 1967.

41.

Vassiliev, Cohomology of braid groups and the complexity of algorithms, Functional Anal. Appl 22 (1989), 182–190. doi:10.1007/BF01077624.

42.

Weihrauch, Computable Analysis, Springer, Berlin, 2000.

43.

Weihrauch, Computational complexity on computable metric spaces, Mathematical Logic Quarterly 49(1) (2003), 3–21.

44.

Yap,

Sagraloff and

Sharma, Analytic root clustering: A complete algorithm using soft zero tests, in: Conference on Computability in Europe, Springer, 2013, pp. 434–444.

45.

Y.S.

Zavyalov,

B.I.

Kvasov and

V.L.

Miroshnichenko, Methods of the Spline Functions, Fizmatgiz, Moscow, 1980, (in Russian).

46.

J.-P.

Zhou, On the degree of extensions generated by finitely many algebraic numbers, Journal of Number Theory 34 (1990), 133–141. doi:10.1016/0022-314X(90)90144-G.

47.

Ziegler, Real computation with least discrete advice: A complexity theory of nonuniform computability, Annals of Pure and Applied Logic 163(8) (2012), 1108–1139. doi:10.1016/j.apal.2011.12.030.

48.

Ziegler and

Brattka, A computable spectral theorem, in: Proc. CCA 2000, Lecture Notes in Computer Science, Vol. 2064, 2001, pp. 378–388.