Efficient and Generic Algorithm for Rigorous Integration Forward in Time of dPDEs: Part I

We propose an efficient and generic algorithm for rigorous integration forward in time of partial differential equations written in the Fourier basis. By rigorous integration we mean a procedure which operates on sets and return sets which are guaranteed to contain the exact solution. The presented algorithm generates, in an efficient way, normalized derivatives which are used by the Lohner algorithm to produce a rigorous bound. The algorithm has been successfully tested on several partial differential equations (PDEs) including the Burgers equation, the Kuramoto-Sivashinsky equation, and the Swift-Hohenberg equation. The problem of rigorous integration in time of partial differential equations is a problem of large computational complexity and efficient algorithms are required to deal with PDEs on higher dimensional domains, like the Navier-Stokes equation. Technicalities regarding the various optimization techniques implemented in the software used in this paper will be reported elsewhere.


Introduction
The main goal of this paper is to present an efficient algorithm for rigorous integration of dissipative partial differential equations (dPDEs).By rigorous integration we mean a process in which a compact and connected set containing the unique solution is propagated along the time interval.The rigorous integration process guarantees to include the exact solution in the output set.Whereas in any standard (non-rigorous) integration process a point, represented by floating point numbers, is being propagated without indication how close is the approximated solution to the real one.
Our main objective in the rigorous integration task is to perform computer assisted proofs of the existence of certain solutions: attractors, periodic orbits, connections between fixed points etc.
Sometimes, when one is willing to perform a computer assisted proof of the existence of a particular solution, the integration process can be avoided.Advantage of the methods which do not perform numerical integration is that the error involved in computing rigorously a solution using a boundary value setup is more distributed along the whole orbit.Examples of rigorous computational methods to study PDEs, different than the approach considered here-i.e.not involving the integration process include [3,9,10,15], and [25].
The main motivation behind this work was to reduce the actual execution time of a sequential algorithm for rigorous integration of dPDEs.We emphasize that the problem of rigorous integration of dPDEs is, by its nature, hard to parallelize, therefore any reduction of computational complexity is highly appreciated.
We have performed several rigorous numerics tests to show that -the presented algorithm is apparently more efficient than the currently used algorithm for rigorous numerical integration of dPDEs used e.g. in [7] and [35], -the presented algorithm produces negligibly worse bounds than the currently used algorithms.
The proposed approach in this paper combines several techniques.In this paper we briefly present issues arising when dealing with rigorous calculations and present methods which we used to circumvent the issues.We briefly present techniques used to obtain the algorithm generality, i.e. to make it easily adaptable to different equations and boundary conditions that could still be studied using Fourier series.More details and the proofs are presented in [8].
We restrict ourselves to the following problem u : [0, T ) × T → R, f, u 0 : T → R, u t = L(u) + N (u, u x , u x x , . . . ) + f, u(0, x) = u 0 (x), (1) where L is a linear elliptic operator (for instance L(u) = νu x x or L(u) = −νu x x x x − u x x ), N (u, u x , u x x , . . . ) is a proper polynomial of u and its partial derivatives, f is a constant in time forcing, T is a maximal time of existence of the solution, T is the one-dimensional torus.
We restrict N (u, u x , u x x , . ..) to be a proper polynomial, as not any equation ( 1) is a dPDE.It is additionally required from (1) to be called a dPDE that the linear part must "dominate" in some sense the nonlinear part, for the formal definition refer to [35].
The Fourier basis is introduced and the problem of solving (1) is reduced to the problem of solving the following infinite dimensional dynamical system where λ k are eigenvalues of the elliptic operator, c N , n, and r j depend on N (u, u x , u x x , . . .), {a k } k∈Z ⊂ C describes the evolution of the Fourier modes.Let us mention few important equations which fall into the dPDE class: the Burgers equation, the Kuramoto-Sivashinsky equation, the Ginzburg-Landau equation, the Swift-Hohenberg equation and, last but not least, the Navier-Stokes equation.
The system (2) is an infinite dimensional system of ordinary differential equations (ODEs), and rigorously integrating it on a computer requires special techniques.To the best of our knowledge, there exist two published theoretical approaches to the rigorous numerics solution for the non-stationary PDE problem-the method of self-consistent bounds, presented in the series of papers [32,[34][35][36], and the method presented in [1].
Both of the mentioned methods have been successfully applied to the study of the Kuramoto-Sivashinsky equation dynamics (KS equation).In fact, the existence of some time-periodic solutions has been proved.The focus of the previous research on one particular equation comes from the interest in establishing the first proof of the existence of spatio-temporal chaotic dynamics in a partial differential equation.
The presented algorithm is a first attempt to provide a generic and robust tool (a computer software) for the rigorous analysis of any dPDE, not only the KS equation.
Moreover, the proposed algorithm improves the efficiency of the rigorous integration task, which is of great importance for us, as the attempt to prove the existence of chaotic dynamics in the KS equation by using our method will require performing a lot of proofs in parallel.We hope that the presented work will permit computer assisted proofs concerning the dynamics of higher dimensional dPDEs, including the Navier-Stokes equations.We consider especially interesting 3D Navier-Stokes equations, because for this case even the fundamental question of the existence and uniqueness of solutions has no definitive answer yet.It is worth noting that there already exist results in literature about the rigorous numerics for higher dimensional PDEs, see e.g.[13,20], and [22], but in this work the stationary problem is considered only, the equations are not being numerically integrated at any point.
The conclusion that we draw is that the algorithm is a proper approach to the task of the rigorous integration of higher dimensional dPDEs, including the Navier-Stokes equations, which we will address in our forthcoming work.

Fast Fourier Transform algorithm
Any attempt to integrate in time an equation of the type (2) will require a significant number of convolution calculations.It is well known that calculation of this type of convolutions is done efficiently by using the Fast Fourier Transform (FFT) algorithm.Fields of applicability of the FFT are, for instance, numerical solving of partial differential equations for weather prediction, a signal analysis and image processing in computer graphics.To the best of our knowledge, the FFT algorithm was described for the first time in the paper [6].We refer to the papers [29] and [30] for an excellent review.In this case it is enough to know how to calculate a finite vector field approximately and the application is straight-forward.One would think that it is so in the case of rigorous methods too.In fact, in rigorous methods the situation is much more difficult.We address here the questions related to the application of the FFT algorithm within a rigorous method setting, including how the interval FFT algorithm behaves, whether the wrapping effect (see [33] for the statement of the wrapping effect problem) will be present, how to calculate the higher order time derivatives required by the Taylor method, and what optimizations to perform.We answer all these questions here, we introduce a new technique of eliminating the aliasing error, and in [8] we present the proofs and details.FFT was already used to study dPDEs by the rigorous numerics community in [11,12], and [16,17].

Outline of Approach
Notation From now on we assume that N is a positive natural number not to be confused with the polynomial N (u, u x , u x x , . ..) in (1).
When using a spectral method based on the Fourier series one often has to compute a considerable amount of convolutions of a finite sequence where n is the number of terms in the convolution u * It is known in the numerical analysis community [5] that computation of (3) could be considerably accelerated by means of the FFT algorithm.The idea of the approach is based on the following facts -the function u(x) = N k=−N u k exp (ikx) is in the unique way determined by u(2π j/M), where M ≥ 2N + 1 and j = 0, 1, . . ., M − 1, -passing from the values at the nodes u(2π j/M), j = 0, . . ., M − 1 to the Fourier coefficients u k , k = −N , . . ., N and vice-versa is done efficiently by using the FFT algorithm.

Uniform Matrix Form of the FFT
Definition 1 Let v : R → C be a 2π-periodic function such that v belongs to the subspace spanned by the functions {exp When it is clear from the context we will drop the notation of N , and use simply the l 2 coefficients of v. Definition 2 Let v : R → C be a 2π-periodic function such that v belongs to the subspace spanned by the functions {exp When it is clear from context we will drop the notation of M, and use simply the L 2 coefficients of v.
By l 2 (N ) → L 2 (M) transform (later on called l 2 → L 2 transform) we mean simply the evaluation of û j M−1 0 from {u k } N −N , being in fact the Discrete Fourier Transform From now on we assume M ≥ 2N + 1.By L 2 (M) → l 2 (N ) transform (later on called L 2 → l 2 transform) we mean the computation of the Fourier coefficients The idea behind the FFT algorithm is that (4) and ( 5) are unified in a matrix form, and then the resulting matrices are written as a product of sparse matrices.Let us fix N and M. We map the l 2 coefficients {u k } N k=−N and the L 2 coefficients { û j } M−1 j=0 to M dimensional vectors z and ẑ in the following way Then the matrix form is given by the following Lemma Lemma 1 (the FFT matrix form) where W M ( j, k) := exp(i jk2π/M), W M ( j, k) := 1 M exp(−i jk2π/M).To evaluate (7) efficiently one does not perform the multiplication by W M , which is a dense matrix but one factorizes W M into sparse matrices and perform the multiplication by sparse matrices.

Issues arising when the FFT algorithm is used within interval arithmetic framework
In this section we describe what the aliasing error is and we provide techniques that allow to eliminate it.
We emphasize that the technique presented in this section is apparently more efficient than the standard one, currently used in the rigorous numerics community, see [11] and references supplied there.Our approach is based on the technique introduced in [26].The authors of [26] introduced new shifted discrete grids in order to eliminate the aliasing error in three dimensional convolutions arising in the Navier-Stokes equations.As we present here, shifted discrete grids can be applied to rigorous methods for dPDEs in 1D and apparently yield a substantial improvement.

Aliasing Problem
When the convolutions like (3) are calculated using the approach outlined in Sect.2.1 the so-called aliasing problem appears.

Convolution of n terms We set
Consider now the convolution of n terms, u * • • • * u.Then the result of the multiplication of the L 2 coefficients is the following Observe that ŝ j M−1 j=0 are L 2 -coefficients of the following function is the aliasing error.The range of a appearing in the sum of the aliasing error is bounded, namely if a M > (n + 1)N , then for all k ∈ {−N , . . ., N − 1, N } holds k + a M > n N and k − a M < −n N .

Definition 3
We call sk 's as above the aliased convolutions.

Standard Technique of Aliasing Error Removal (Padding)
The aliasing error is eliminated simply by choosing the number of grid points M large enough.This is the approach used in [11] and other related works.The following condition has to be satisfied in the general case we take the worst case condition and obtain Therefore when M > (n + 1)N the aliasing error is not present.123

Technique Based on Phase Shifts
This technique is designed only for the functions, for which the Fourier coefficients {u k } are purely real or purely imaginary (for instance real odd or real even functions).To remove the aliasing error one may also use a phase shift, that is calculate L 2 coefficients at the points belonging to the uniform grid shifted to the right by a constant value . Let The particular choice of = π/(2M) is interesting for our purposes.The following lemma holds, a proof is provided in [8].= s k (for u odd and n odd), for all −N ≤ k ≤ N .
Now, let us put the presented techniques into the context of the FFT algorithm.Whatever aliasing elimination technique is used in the FFT approach to calculate convolutions, the M ≥ 2N + 1 relation has to be in any case satisfied to allow all modes (whose number is 2N + 1) to be mapped by the mapping (6).
Although the difference is asymptotically negligible, for the input sizes which are attainable in practical applications the number of operations in this technique compared to the technique from 3.2 is reduced by a significant factor.This is due to a reduction of the aliasing error influence in this technique-the values in the reduced number of nodes are calculated.

Overestimates
As we are interested in rigorous computations, therefore all computations have to be performed using the interval arithmetic [24,31].In our implementation we used the infimumsupremum representation of intervals.In this subsection we address the question if the rigorous evaluation of the convolutions (3) using the FFT algorithm produces more overestimates than the direct evaluation.We will refer to the interval FFT algorithm as the rigorous FFT algorithm.
We have not found a comprehensive study of the overestimates that arise in the rigorous FFT algorithm in the existing literature.We have performed numerical tests in order to provide clear answer to the question if the rigorous evaluation of the convolutions (3) using the FFT algorithm produces more overestimates compared to the direct evaluation.
In the test presented in this section the function u has been a real-valued and odd function belonging to the space spanned by a finite number of vectors, randomly picked for each single test, all u k have been intervals.We have fixed C > 0, s > 0 and required that {u in order to mimic regular functions.The diameter (width) was the same for all u k .In each test we used the techniques presented in Sects.3.2 and 3.3 in order to eliminate the aliasing error.We present example results in Tables 1, 2, and 3.
Apparently, the rigorous FFT algorithm produces more overestimates, especially when sets of large diameter (of order > 1) are used.The considerably larger overestimates produced by the FFT compared to the direct method are probably due to several matrix-vector multiplications performed by the FFT algorithm, it is known that matrix-vector multiplications are the main source of the so-called wrapping effect, for precise statement of the wrapping effect refer e.g.[23], and [33]. x The goal of a rigorous ODEs solver is to find a compact and connected set we denote the solution of (8) at the time t k with initial condition x 0 ∈ R n , and therefore ϕ(t k , x 0 ) denotes the set of all the values which are attained at the time t k by any solution of (8) with the initial condition in x 0 .
Notation We denote by There exist a few algorithms that offer reliable computations of the solution trajectories for ODEs that are based on interval arithmetic.The approach used in this paper is based on the Lohner algorithm, introduced in [21], see also [33].To obtain a rigorous bound for the solution and avoid at the same time the so-called wrapping effect we calculate a bound for the partial derivative with respect to the initial conditions.Basically, in the Lohner algorithm (see e.g.[33]) instead of calculating directly the expression the mean value form is used (10) h > 0 is a time step, Φ q is a q-th order numerical method (for the purpose of this paper, the Taylor method) for an equation (8), R Φ q is a remainder term calculated using a convex hull of an enclosure for the solution values in the time interval [0, h], i.e. a set Bounds for the set and bounds for the set containing partial derivatives with respect to the initial condition are both obtained by the Taylor method in our approach.To obtain bounds for (12) we will use a suitably modified Taylor method.Here we do not present in detail the Taylor method and the Lohner algorithm and we refer the reader to [21,27] and [33].
Obtaining Taylor coefficients of any order Now, we address the question of how to obtain bounds for the Taylor coefficients of any order efficiently.The algorithms from [33] used some explicit formulas for the Taylor coefficients of any order.These explicit formulas were derived for the second degree polynomials only.This approach was used in the following work considering rigorous dPDE integration, i.e. [32,[34][35][36].We found this approach awkward, mainly because of its troublesome implementation and limited applicability.Apparently, this approach works only for dPDEs with the non-linear term being a second degree polynomial.
Instead, we found a more suitable to use algorithm which combines the jet transport and the automatic differentiation techniques.These techniques have already proven to be extremely useful in many applications including a long-time integration of the solar system (see [18] and references given therein) or computer assisted proofs for ODEs.

Automatic Differentiation
Now, we address the question of how to obtain the derivatives with respect to time of a sufficiently smooth function efficiently.Here we present a convenient and yet efficient approach called the automatic differentiation, see e.g.[14,24].The automatic differentiation has many applications beside of our interest.It is a general technique of obtaining recursively a fixed number of normalized derivatives of a general operation/procedure.

Definition 4
Let n > 0, a : [0, t max ) → R be a sufficiently smooth function, we call its n-th normalized derivative.
Let b, c : [0, t max ) → R be sufficiently smooth functions, assume that a(t) = G(b(t), c(t)).Using the automatic differentiation the normalized derivatives of a(t) can be obtained assuming that G is a composition of some basic operations for which the recursive evaluation formulas are known.Recursive formulas for basic operations (such as multiplication, exponential, logarithm, trigonometric functions etc.) are easily derived, see for example [18].
Only particular systems originating from a dPDE with a polynomial non-linearity are in the scope of this paper.Rules for the addition and the multiplication of two functions are required only, which will became apparent later on.The scalar multiplication, being trivial, is not presented here.Assume that b [ j] , c [ j] are known for j = 0, . . ., r then Now let us consider a differential equation Given x [0] = x (0) we calculate x [ j] for j = 1, . . ., q − 1, by the following recursive formula where (F • x) [ j] is the j-th normalized derivative of the right-hand side function of (15) (F composed with x as a function of t).

Jet Transport
Now, we address the question how to obtain bounds for the partial derivative with respect to the initial condition (12).The most convenient approach is to modify the automatic differentiation procedure to make it calculate ( 11) and ( 12) simultaneously.This is done by the jet transport technique.
An algebra of the first order jets is formed with the addition and multiplication operations defined as follows Definition 6 Let n > 0, a(ξ ), b(ξ ) be first order jets, ξ = (ξ 0 , . . ., ξ n−1 ) are symbols.We define the addition + of first order jets as follows We define the multiplication • of first order jets as follows Definition 7 Let f : R n → R, t ≥ 0 be a time, x 0 ∈ R n be a fixed value, ξ = (ξ 0 , . . ., ξ n−1 ) are symbols.We define the following first order jet for the function f The definition above has a straightforward extension to the case when f is vector-valued.When f is vector-valued, the vector of jets for f is denoted by J f (ξ ) .Using the definition above we fix t ≥ 0, x 0 ∈ R n , define the first order jets for ϕ(t, x 0 ) and calculate ϕ [ j] (t, x 0 ), j = 1, . . ., q − 1 using the automatic differentiation procedure described in Sect.4.1, where ϕ(t, x 0 ) is the smooth solution of (15).Precisely In what follows in addition to the regular Taylor method (Φ q ), we are going to use an extended variant (Φ jet q ) whose coefficients are J ϕ [0] (ξ ) , . . ., J ϕ [q−1] (ξ ) .Φ jet q works as follows: the base part of J ϕ [0] (ξ ) is set to the value of the current initial condition and the variational part of J ϕ [0] (ξ ) is set to be the identity.The higher order coefficients are calculated performing the same operations as in the regular method, but on first order jets instead of scalars, and operations are performed with the rules of the algebra defined in Definition 6.This allows to solve (15) and the variational equations corresponding to it simultaneously.

An Algorithm Combining Automatic Differentiation with FFT
In this section we propose an algorithm for calculating both the solution of systems of ODEs originating from a dPDE and a matrix of partial derivatives with respect to the initial condition.
We emphasize the fact that the output of the presented algorithm is not enough to integrate efficiently a given system forward in time.This is not a stand-alone algorithm, nonetheless the output of the presented algorithm is necessary for any effective rigorous solver based on the Taylor method.Later we have combined it with the Lohner algorithm [21] in order to integrate example systems forward in time rigorously and we investigate how the techniques interplay with each other and what is the quality of the results.
The algorithm has been designed for dPDEs written using the truncated Fourier basis.We consider the following abstract system denote all values of the solution within the time interval [0, h], and starting with the initial condition in From now on let ϕ be regular with respect to time solution of (16).Observe that the coefficients {a k } N k=−N are complex in this case, and ϕ : [0, t max ) × C 2N +1 → C 2N +1 .First, let us comment on the possible approaches used to deal with the first order jets of complex functions.First possible approach is based on the following observation ) to be the following projection Re (a −N , . . ., a N ) = (Re (a 0 ) , Im (a 0 ) , . . ., Re (a N ) , Im (a N )) .

Observation 1 Any function
By using Observation 1 we are able to translate our system (16) and use real first order jets in order to calculate the partial derivative with respect to the initial condition of the numerical method.But then a lot of information about the structure of the system is lost, and we do not recommend to use this approach.
The approach that we use is the following.Instead we use "complex" first order jets, which have the following form, let We claim that using this kind of jets greatly facilitates computations and permits certain optimizations.This approach is strictly of computational interest.Below we provide a sketch of a justification of the presented approach.
Let a(ξ ) and b(ξ ) be arbitrary "complex" first order jets for two complex variables Let us take the corresponding real first order jets c = Re (a), d = Im (a), e = Re (b) and f = Im (b).The corresponding first order jets are A simple calculation shows that if and Let us comment on what the main ingredients of the described algorithm are.First, the information about the partial with respect to the initial condition derivative is obtained by introducing the first order jets in place of ordinary scalar values.Second, the first order jets are passed through the FFT algorithm in order to minimize the number of operations needed to calculate the convolutions.Third, in order to allow a fast calculation of the automatic differentiation convolutions the L 2 coefficients of the normalized derivatives are calculated and stored at each step.

Algorithm 1 Outline of an algorithm for efficient and rigorous calculation of one step of Φ
jet q , dedicated for any system of the type (16).Input -the Galerkin projection dimension N > 0, -number of points used by FFT, such that the aliasing problem is avoided (see Sect. 3.1) M ≥ 2N + 1, -an interval set of initial conditions [x 0 ] ⊂ C 2N +1 , -time step h > 0, -order of the Taylor method q > 0.

. , N [ j] M (ξ ) , (c) add the linear contribution F
[ j] −N (ξ ) , . . ., −N (ξ ) , . . ., F  (10).The matrix composed of the variational part of the vector of first order jets (10).The other terms in (10), i.e.Φ q (h, mid ([x])), R Φ q ([W ]) are obtained by two separate Taylor methods, as different initial conditions has to be used.The vectors Φ q (h, mid ([x])), and R Φ q ([W ]) are evaluated using scalar values, not first order jets, and when N is large the cost of calculating them is negligible compared to the cost of calculating Φ jet q (h, [x 0 ]).

Now we would like to comment on
Step (a) of the "for" loop in Algorithm 1, because that is where savings are introduced by the FFT algorithm.
The improvement of introducing the FFT algorithm into the standard procedure is as follows-at the given order j instead of calculating the double convolutions directly, the following operations, comprising steps 2a, 2b and 2e of Algorithm 1 or equivalently comprising steps (a), (b) and (e) of Algorithm 1 are used, compare with the outline presented in Sect.2.1: -one l 2 → L 2 transform, and one L 2 → l 2 transform, -calculation of Observe that in the FFT variant the loop over k is not present at all.We discuss in [8] several optimizations that can be introduced into this algorithm.

Rigorous Numerics Tests
Below we present a report from various tests we performed with the developed methods and algorithms.We emphasize that by tests we mean rigorous numerics tests.The purpose of the developed methods and algorithms are the computer assisted proofs for dPDEs.In the tests either a finite truncation (a Galerkin projection) or the full infinite system (giving a differential inclusion) were considered.In order to integrate in time the full infinite dimensional system of ODEs (2) we split it into two parts where F(x + y) is the right-hand side of (2), We assume the embedding X N {x −N , . . ., x 0 , . . ., x N } ≡ {0, . . ., 0, x −N , . . ., x 0 , . . ., x N , 0, . . ., 0} ∈ C ∞ , and use the same symbol to denote both of the elements of the finite dimensional space and the infinite dimensional space.We bound the solutions of ( 2) by considering the following differential inclusion where δ ⊂ X N describes the influence of y onto P N F(x + y).The Lohner-type algorithm for rigorous integration in time of differential inclusions was originally provided in [19].In this algorithm the solutions of a Cauchy problem associated with ( 21) are bounded, and then the influence of the perturbation δ is a-posteriori added.In order to calculate δ the FFT algorithm may also be used.We specified what system was being integrated for each test separately by stating a "Galerkin projection" or a "differential inclusion" respectively.We used the implementation of the part of the Lohner algorithm regarding the Lohner representation of the sets from the CAPD library [4].
In this section we present results from all the tests performed in the form of tables.Whenever a diameter of a finite or an infinite dimensional set is provided in a table, the meaning of this number is in fact the maximum over coordinates of all diameters.The labels "Direct approach" "FFT approach" and "FadBad++" appearing in the tables indicate different implementations of the Taylor method that were used."Direct approach" indicates that the normalized derivatives were obtained directly (18), whereas "FFT approach" indicates that the normalized derivatives were obtained by the FFT approach presented previously.To remove the aliasing error technique the "II technique based on phase shifts" was used, unless otherwise stated."FadBad++" indicates that the normalized derivatives were obtained by using the FadBad++ software library for automatic differentiation [2], without using the FFT at any point.The vector field, which was input to the FadBad++ library was optimized, by using Observation 1 and all possible symmetries.The part of the Lohner algorithm regarding the Lohner representation of the sets was the same for both.For a thorough description of the Lohner algorithm and the issues related to the wrapping effect and representation of sets refer to [33].The meaning of the symbol M depends on the context; in [35] and [7] it is used to denote the dimension of the so-called finite tail, when the full infinite dimensional system is considered, whereas in this paper it is used to denote the number of points of the discrete grid used by the FFT.To avoid ambiguities in this section, we denote the number of points of the discrete grid used by the FFT by M F FT .

Burgers Equation Fixed Point Test (Galerkin Projection)
When periodic boundary conditions are used and the forcing function is not present, the Burgers equation exhibits a globally attracting fixed point at zero.This is proved by the energy exponential decay, see [7].
With this test we address the question of what are, for the overall Lohner algorithm, the consequences of considerable overestimates produced by the FFT approach (see the results from Sect.3.4).This question is difficult to answer without performing actual calculations due to the sophisticated nature of the Lohner type rigorous integrator.The setting is The uniform zero-centred box (with the fixed point at the centre) was rigorously integrated forward in time on the fixed time interval with a fixed time step.The diameter of the initial box was the parameter of the tests (Fig. 1).
In Table 4 we present a report from the performed tests.The diameter of the result, which was, in fact, the diameter of the set at the time 10 (10, 000 constant time steps of length 0.001 were executed), was measured.Due to the fact that the fixed point is attracting, the initial box is expected to decrease in diameter.We were in particular interested in the difference between the quality of the FFT approach result, and the direct approach result.The results presented in Sect.3.4 predict much better quality of the direct approach result.
The parameters that were used for the tests were ν = 0.1, order of the Taylor method was 6, constant time-step h = 0.001, N = 19, M F FT = 40.
Total running time FadBad++ FFT approach 1849 sec 424 sec (23) See (23) for the comparison of the program running time.

Kuramoto-Sivashinsky Equation Attracting Periodic Orbit Test
The setting is We picked one of the periodic orbits proved to exist in [35] (the attracting periodic orbit for ν = 0.032).We took a box around an initial condition from [35] for this periodic orbit, and then the box was rigorously propagated forward in time using two implementations of the Lohner algorithm.In the tests the equations were integrated until the time of the full revolution of orbit was reached, which is approximately 0.4091.The results of these tests are contained in Tables 5 and 6.

I test-efficiency (Galerkin projection)
All of the techniques presented previously in this part of the article can be replaced by applying some of the available software libraries for In this test a blow-up of the set was experienced for the FFT approach.By a blow-up we mean that the rigorous bounds escaped the range of the representable numbers.Such bounds are then represented in the software by the unbounded interval ([−∞, ∞]).As soon as the unbounded interval appears in the calculation process, further calculations are meaningless.We found a simple method circumventing this problem at a little additional computational cost.To avoid the blow-up problem, the first order normalized derivative must be calculated by using the direct approach, nonetheless the higher order derivatives can be obtained by using the FFT approach.The cost of calculating the first order normalized derivative is negligible for the overall cost of calculating the normalized derivatives up to a reasonable fixed order The machine used was Linux 32-bit Intel Core i5-2430M CPU @ 2.40 GHz x 4 automatic differentiation, for instance Fad-Bad++ [2].The purpose of this test was to check if the presented approach will eventually be useful in calculating a rigorous bounds of the Poincaré map, which are needed e.g. to verify the Brouwer fixed point theorem assumptions.With these tests we would like to provide an argument that we did not perform redundant work developing this approach.The Lohner algorithm implementation based upon the FFT approach allowed to obtain bounds of the same order, and at the same time this bounds were obtained faster.
With this test we compared the efficiency of the Lohner algorithm implementation based upon the presented FFT approach with the implementation based on the direct approach implemented using the external software library for automatic differentiation FadBad++ [2] (the approach used to perform the computer assisted proof presented in [7]).
A Galerkin projection of the infinite dimensional system was considered in these tests in order to isolate this part of the overall algorithm to which the developed techniques contribute.
The initial condition was taken from [35], the diameter of the initial condition was 8e − 05 in each direction.Parameters used for this test were as follows ν = 0.032, order of the The full infinite dimensional system was considered (differential inclusion) a In this test a blow-up of the set was experienced for the FFT approach.By a blow-up we mean that the rigorous bounds escaped the range of the representable numbers.Such bounds are then represented in the software by the unbounded interval ([−∞, ∞]).As soon as the unbounded interval appears in the calculation process, further calculations are meaningless.We found a simple method circumventing this problem at a little additional computational cost.To avoid the blow-up problem, the first order normalized derivative must be calculated by using the direct approach, nonetheless the higher order derivatives can be obtained by using the FFT approach.The cost of calculating the first order normalized derivative is negligible for the overall cost of calculating the normalized derivatives up to a reasonable fixed order Taylor method was 5, constant time-step h = 0.00015, N = 23 (all the same as in [35]), M F FT = 48.
Observe that the time step had to be adjusted for each test, because of the stiffness problem exhibited by the Kuramoto-Sivashinsky equation.The results of this test are contained in Table 5.

II test-overestimates (differential inclusion)
The diameter of the propagated set was measured after the full revolution along the orbit.The full rigorous integrator (the Lohner algorithm for differential inclusions) of the infinite dimensional system was used in this test.
Parameters used for this test were as follows ν = 0.032, order of the Taylor method was 5, constant time-step h = 0.00015, N = 23, M = 92 (the dimension of the finite tail, for the exact meaning of this number refer to [35] or [7]), M F FT = 48.The results of this test are contained in Table 6.
We picked one of the unstable periodic orbits proved to exist in [35] (the unstable periodic orbit for ν = 0.02991).This orbit is on the apparently chaotic attractor.We took a box around an initial condition from [35] for this periodic orbit, and then the box was rigorously propagated forward in time using two implementations of the Lohner algorithm.In the tests the equations were integrated until the time of the full revolution of the orbit was reached, which is approximately 0.4490232.The results of these tests are contained in Table 7.
The diameter of the propagated set was measured after the full revolution along the orbit.The full rigorous integrator (the Lohner algorithm for differential inclusions) of the infinite dimensional system was used in this test.
Parameters used for this test were as follows ν = 0.02991, order of the Taylor method was 5, constant time-step h = 0.00015, N = 25, M = 75 (the dimension of the finite tail, for the exact meaning of this number refer to [35] or [7]), M F FT = 64.The orbit is on the apparently chaotic attractor.The full infinite dimensional system was considered (differential inclusion) a In this test a blow-up of the set was experienced for the FFT approach.By a blow-up we mean that the rigorous bounds escaped the range of the representable numbers.Such bounds are then represented in the software by the unbounded interval ([−∞, ∞]).As soon as the unbounded interval appears in the calculation process, further calculations are meaningless.We found a simple method circumventing this problem at a little additional computational cost.To avoid the blow-up problem, the first order normalized derivative must be calculated by using the direct approach, nonetheless the higher order derivatives can be obtained by using the FFT approach.The cost of calculating the first order normalized derivative is negligible for the overall cost of calculating the normalized derivatives up to a reasonable fixed order

Swift-Hohenberg Equation a Connection Between Fixed Points Test (Differential Inclusion)
The setting is Looking at the eigenvalues of the linear part of ( 25) it is immediately verified that zero is an unstable fixed point for (25a), for a sufficiently large ν.Apparently, there are also other fixed points, which are attracting.

Overestimates (differential inclusion)
As the initial condition, we took a small interval set of functions close to zero, i.e.This set of functions had been rigorously integrated until it was caught in the basin of attraction of an attracting fixed point for (25a).The process of integration of the set was stopped as soon as the set was attracted by and centred at an attracting fixed point location (Fig. 2).The order of the Taylor method used for the tests was 5.The results of this test are contained in Table 8.
Numerical data from all the tests presented in this section is available on-line [28].a In this test a blow-up of the set was experienced for the FFT approach.By a blow-up we mean that the rigorous bounds escaped the range of the representable numbers.Such bounds are then represented in the software by the unbounded interval ([−∞, ∞]).As soon as the unbounded interval appears in the calculation process, further calculations are meaningless.We found a simple method circumventing this problem at a little additional computational cost.To avoid the blow-up problem, the first order normalized derivative must be calculated by using the direct approach, nonetheless the higher order derivatives can be obtained by using the FFT approach.The cost of calculating the first order normalized derivative is negligible for the overall cost of calculating the normalized derivatives up to a reasonable fixed order b This is the same test as in 2(b).In this test a blow-up was experienced.It was possible to eliminate it by using the technique described in foot note a like in other cases, but with proper choice of M-the number of discrete points used by FFT (in this case for M = 32 the blow-up was eliminated, whereas for M = 30 the blow-up was still experienced) c This is the same test as in 2(a).In this test a blow-up was experienced.It was eliminated when further the padding technique was used instead of the phase-shift aliasing error elimination technique regardless the choice of M-the number of discrete points used by FFT

Complexity
In this section we analyse the complexity of the direct approach and the FFT approach, which were compared in numerical tests presented in Sect. 5. First, we compare the asymptotic complexity of one step (generating the coefficients and evaluating the polynomial) of the extended Taylor method of the order q (Φ jet q ).Coefficients of Φ jet q are the first order jets of a vector valued function.Second, we compare the actual number of elementary operations required to perform one step in both of the approaches.We calculated the actual number of operations, motivated by the fact that asymptotic complexity provides only a rough information about the magnitude of inputs that can be handled by using the current computing power, and we have been interested explicitly how much faster we can calculate.Fig. 3 The number of interval multiplications (divided by 10 6 ) required to perform one step of Φ jet q in the Burgers equation setting (22) with periodic boundary conditions We use the well known fact that the FFT algorithm's asymptotic complexity is O(n log n), where n is the size of the input.The asymptotic number of operations required to perform one step of Φ jet q using the FFT approach is FFT s on jets c 1 q N M log M + multiplications of the L 2 coefficients (jets) in FFT s c 2 q(q + 1)M N, (26) whereas the complexity of one step of Φ jet q using the direct approach is "double" convolution on jets where c 1 , c 2 , c 3 > 0 are constants, and M is a number depending linearly on N , such that M ≥ 2N + 1, see Sect.3.1.The term N M log M in (26) in place of M log M seems strange at the first sight, but it is due to the fact that the input and the output for all the FFT's are vectors of first order jets.The number of elementary operations required to perform an arithmetic operation on first order jet depends linearly on N , as one loop is requiredcompare Definition 6.For the meaning of "the FFT approach" and "the direct approach" refer to Sect. 5.
To confirm ( 26) and ( 27) empirically we counted the exact number of interval multiplications and interval additions, which are performed by the computer program during one step of Φ jet q in the Burgers equation setting (22).The results are presented in Fig. 3. Observe that the calculations are performed within the interval arithmetic framework.An interval arithmetic operation is considerably more costly than the corresponding ordinary floating-point arithmetic operation, e.g. to perform an interval multiplication usually four floating-point multiplications and two rounding mode switchings are required.Overall tests indicate that the execution time of an interval arithmetic operation is approximately of the order ten times the time of the corresponding floating-point arithmetic operation.
We also counted the number of interval additions and multiplications for the other two settings, i.e. (24) and (25).But the corresponding figures are similar in spirit to Fig. 3, so we Table 9 Exact number of elementary operations on complex numbers required to calculate the j-th coefficient of Φ jet q , which is one step of the "for" loop in Algorithm 1 The calculations were performed in the Burgers equation setting (22), p is a number such that M = 2 p do not present them here and refer the interested reader to the corresponding numerical data [28].

Approach Additions Multiplications
In order to get an idea about the magnitude of constants c 1 , c 2 , c 3 in ( 26) and ( 27) we provided in [8] the exact number of elementary operations on complex numbers required to calculate the j-th coefficient of Φ jet q .The calculations were performed in the Burgers equation setting (22), and assuming that M = 2 p .We present the results in Table 9.

Conclusion
We would like to present some concluding remarks from the results presented above.
The diameters of the sets obtained in the performed tests (Tables 4, 6, and 8) indicate that, surprisingly, the Lohner algorithm seems to be immune to the considerable overestimates produced by the FFT approach.Consequently, we are convinced that the Lohner algorithm based on the developed here FFT approach can be successfully applied to the problem of rigorous integration of higher dimensional PDEs, including the Navier-Stokes equations.And the presented algorithm perfectly suits performing proofs of an existence of solution on a fixed time interval, for instance when one needs to establish existence of a Poincaré map only.
The different test cases have demonstrated that the software based on the developed techniques has improved the efficiency of the rigorous integration task, without severe quality drop.
All the test cases were successfully finished when the new approach was used (no interval blow-ups were experienced), and there was not any considerable quality drop (measured by diameter of the propagated interval sets).We consider this a significant achievement in the context of eventual computer assisted proofs for higher dimensional PDEs.The task of rigorous integration of a higher dimensional PDEs like the Navier-Stokes equation, which is our goal, will require a large amount of computing time.Any reduction of the time spent on computations is of great importance here, because this could allow some phenomena to be actually proven.
Different equations and different boundary conditions were used in each of the presented tests (the Burgers equation with a pure periodic boundary conditions, the KS equation with a periodic-odd boundary conditions and the SH equation with a periodic-even boundary conditions) in order to demonstrate that the developed software is generic and robust.The software can be applied for many problems involving dPDEs.
The goal of achieving both generic and efficient software for rigorous integration of dPDEs has been achieved to some extent.It is a difficult task in general.We have been able to achieve this for periodic boundary conditions in 1D.Still, such software had not existed, and all the results published so far (e.g.[1,35]) have focused on a particular case of the Kuramoto-Sivashinsky equation, and the software used there was hard-coded for the KS equation with periodic-odd boundary conditions.
The tests discussed in previous section are in principle repeatable, files with data from the presented proofs and all the source code used are available on-line [28].

Remark 1
the value of Φ jet q (h, [x 0 ]) using the coefficients from the previous step and Horner's rule.end Let us relate the output results of Algorithm 1 to the mean value form formulation

Fig. 1
Fig.1The setting for the Burgers equation

2
Let M be the number of the points on the grid.Let u : R →

Table 1
Tests for N = 21, s k

Table 2
Tests for N = 21, s k

Table 3
Tests for N = 21, s k

Table 4
Fixed point for Burgers equation test results

Table 5
Total running time of programs propagating a box along an orbit for the KS equation (an attracting periodic orbit for ν = 0.032)

Table 6
Diameter (max over coordinates) of boxes obtained by programs propagating a box along orbit for KS equation (an attracting periodic orbit for ν = 0.032)

Table 7
Diameter (max over coordinates) of boxes obtained by programs propagating a box along orbit for KS equation (an unstable periodic orbit for ν = 0.02991)

Table 8
Diameter (max over coordinates) of boxes obtained by programs propagating a box along connection between fixed points for SH equation