Efficient high order algorithms for fractional integrals and fractional differential equations

We propose an efficient algorithm for the approximation of fractional integrals by using Runge–Kutta based convolution quadrature. The algorithm is based on a novel integral representation of the convolution weights and a special quadrature for it. The resulting method is easy to implement, allows for high order, relies on rigorous error estimates and its performance in terms of memory and computational cost is among the best to date. Several numerical results illustrate the method and we describe how to apply the new algorithm to solve fractional diffusion equations. For a class of fractional diffusion equations we give the error analysis of the full space-time discretization obtained by coupling the FEM method in space with Runge–Kutta based convolution quadrature in time.


Introduction
Fractional Differential Equations (FDEs) have nowadays become very popular for modeling different physical processes, such as anomalous diffusion [26] or viscoelasticity [1,25].In the present paper we develop a fast and memory efficient method to compute the fractional integral for a given α ∈ (0, 1).A standard discretization of ( 1) is obtained by convolution quadrature (CQ) based on a Runge-Kutta scheme [18,5] where the convolution weights ω n can be expressed as, see Lemma 9, with e n (•) a function that depends on the Runge-Kutta scheme.For discretizations based on linear multistep methods, see [14].
To compute up to time T = N h using formula (2) requires O(N ) memory and O(N 2 ) arithmetic operations.Algorithms based on FFT can reduce the computational complexity to O(N log N ) [16] or O(N log 2 N ) [9], but not the memory requirements; for an overview of FFT algorithms see [7].Here we develop algorithms that reduce the memory requirement to O(| log ε| log N ) and the computational cost to O(| log ε|N log N ), with ε the accuracy in the computation of the convolution weights.Hence, our algorithm has the same complexity as the fast and oblivious quadratures of [19] and [23], but as we will see, a simpler construction.
The algorithms will depend on an efficient quadrature of (3) for n ≥ n 0 , with a very moderate threshold value for n 0 , say n 0 = 5.As e n (z) = r(z) n q(z) and r(z) = e z + O(z p+1 ), where p is the order of the underlying RK method, this is intimately related to the construction of an efficient quadrature for the integral representation of the convolution kernel with t ∈ [n 0 h, T ].Note that as Γ(1 − α)Γ(α) = π/ sin(πα), h −1 ω n is an approximation of 1 Γ(α) t α−1 , i.e., the kernel of (1).Even though we eventually only require the quadrature for (3), we begin with developing a quadrature formula for (4) for a number of reasons: the calculation for ( 4) is cleaner and easier to follow, such a quadrature allows for efficient algorithms that are not based on CQ, and finally once this is available the analysis for (3) is much shorter.The quadrature we develop for (4) is closely related to the one developed in [13], the main difference being our treatment of the singularity at x = 0 by Gauss-Jacobi quadrature and the restriction of t to the finite interval rather than semi-infinite as used in [13].Both these decisions allow us to substantially reduce constants in the above asymptotic estimates of memory and computational costs.Recent references [27,12,2] also consider fast computation of (1), but do not address the approximation of the convolution quadrature approximation exploiting (3).Our main contribution here is the development of an efficient quadrature to approximate (3) and its use in a fast and memory efficient scheme for computing the discrete convolution (2).
The stability and convergence properties of RK convolution quadrature are well understood, see [18,6].This allows us to apply convolution quadrature not only to the evaluation of fractional integrals, but also to the solution of fractional subdiffusion or diffusion-wave equations of the form ) (0) = 0, j = 0, . . ., m − 1, with β ∈ (0, 2).Here, ∂ β t = I m−β ∂ m t , with m = ⌈β⌉, denotes the Caputo fractional derivative.Solutions of such equations typically have low regularity at t = 0, but a discussion of adaptive or modified quadratures for this case is beyond the scope of the current paper.For a careful analysis of BDF2 based convolution quadrature of fractional differential equations see [8].
To our knowledge, underlying high order solvers for ODEs have been considered for the approximation of (1) only at experimental level in [2,3] and in [23], where a fast and oblivious implementation of RK based CQ is considered for more general applications than (1).The fast and oblivious quadratures of [19] and [23] have the same asymptotic complexity as our algorithm, but have a more complicated memory management structure and require the optimization of the shape of the integration contour.Our new algorithm has the advantage of being much easier to implement, as it does not require sophisticated memory management and the optimization of quadrature parameters is much simpler, and furthermore only real arithmetic is required.The new method is also much better suited for the extension to variable steps -this will be investigated in a follow up work.On the other hand, the present algorithm is specially tailored to the application to (1) and related FDEs, whereas the algorithms in [19,23] allow for a wider range of applications.
The paper is organized as follows.In Section 2 we develop and fully analyze a special quadrature for (1), which uses the same nodes and weights for every t ∈ [n 0 h, T ].In Section 3, we recall Convolution Quadrature based on Runge-Kutta methods and derive the special representation of the associated weights already stated in (3).In Section 4 we derive a special quadrature for (3), which uses the same nodes and weights for every n ∈ [n 0 , N ], with T = hN .In Section 5 we explain how to turn our quadrature for the CQ weights into a fast and memory saving algorithm.In Section 6 we test our algorithm with a scalar problem and in Section 7 we consider the application to a fractional diffusion equation.We provide a complete error analysis of the discretization in space and time of a class of fractional diffusion equations.
2 Efficient quadrature for t α−1 In the following we fix an integer n 0 > 0, time step h > 0, and the final computational time T > 0. Throughout, the parameter α is restricted to the interval (0, 1).We develop an efficient quadrature for (4) accurate for t ∈ [n 0 h, T ].

Truncation
First of all we truncate the integral where τ (L) denotes the truncation error.
Lemma 1.For t ≥ n 0 h and L = A/h we have that Proof.
Assuming A ≥ 1 we can choose However, in practice it is advantageous to use the bound (5) to numerically find the optimal A.

Gauss-Jacobi quadrature for the initial interval
We choose an initial integration interval along which we will perform Gauss-Jacobi integration.
Recall the Bernstein ellipse E ̺ , which is given as the image of the circle of radius ̺ > 1 under the map z → (z + z −1 )/2.The largest imaginary part on E ρ is (̺ − ̺ −1 )/2 and the largest real part is (̺ + ̺ −1 )/2.Theorem 3. Let f be analytic inside the Bernstein ellipse E ̺ with ̺ > 1 and bounded there by M .Then the error of Gauss quadrature with weight w(x) is bounded by where is the corresponding Gauss formula, with weights w j > 0.
Proof.A proof of this result for w(x) ≡ 1 can be found in [24,Chapter 19].The same proof works for the weighted Gauss quadrature as well.We give the details next.
First of all note that we can expand f in Chebyshev series As I Q is exact for polynomials of degree 2Q − 1, we have that where we have used the fact the weights are positive and integrate constants exactly.
We apply Theorem 3 to the case and denote Theorem 4. For t ∈ [0, T ] and any Q ≥ 1 we have the bound Proof.Since f 0 in ( 7) is an entire function, by Theorem 3 we can estimate Let ̺ = e δ with δ > 0. Then the error bound can be written as We now choose δ so that it maximises the function we have a maximum at Using the identities sinh −1 y = log y + 1 + y 2 , cosh x = 1 + sinh 2 x, we derive an error estimate with the above choice of δ: where in the last step above we have used that −1 + √ 1 + x 2 ≤ x for x > 0. This gives the stated result.

Gauss quadrature on increasing intervals
We next split the remaining integral as where −α e −t(y+1)∆L j /2 dy, where ∆L j = L j − L j−1 , j = 1, . . ., J, with L J = L.The intervals are chosen so that for some B ≥ 1, ∆L j = BL j−1 , i.e., L j = (B + 1)L j−1 and J = ⌈log B+1 L/L 0 ⌉.To each integral we apply standard, i.e., w(x) ≡ 1 in Theorem 3, Gauss quadrature with Q nodes and denote the corresponding error by Theorem 5.For any Q ≥ 1 and t ≥ 0 Proof.Note that the integrand f j in ( 8) is now not entire and there will be a restriction ̺ < ̺ max on the choice of the Bernstein ellipse E ̺ in order to avoid the singularity of the fractional power.In particular we require which is satisfied for 1 < ̺ < ̺ max and Setting we see that ε ∈ (0, 1) for ̺ ∈ (1, ̺ max ) and that Hence The result is now obtained by using Remark 6. Choosing for instance ε = 0.1 and B = 1 in the above estimate, we obtain and As we will require a uniform bound for t ∈ [t n 0 , T ], we can substitute t = t n 0 in this estimate.
3 Runge-Kutta Convolution Quadrature Let us consider an s-stage Runge-Kutta method described by the coefficient matrix O ι = (a ij ) s i,j=1 ∈ R s×s , the vectors of weights b = (b 1 , . . ., b s ) T ∈ R s and the vector of abcissae c = (c 1 , . . ., c s ) T ∈ [0, 1] s .We assume that the method is A-stable, has classical order p ≥ 1, stage order q and satisfies a s,j = b j , j = 1, . . ., s, [10].The corresponding stability function is given by where Our assumptions imply the following properties: Important examples of RK methods satisfying our assumptions are Radau IIA and Lobatto IIIC methods.
Let us consider the convolution where K(z) denotes the Laplace transform of the convolution kernel k(t).K is assumed to be analytic for Re z > 0 and bounded there as |K(z)| ≤ |z| −µ for some µ > 0. The operational notation K(∂ t )f introduced in [15], is useful in emphasising certain properties of convolutions.Of particular importance is the composition rule, namely, if .. This will be used when solving fractional differential equations in Section 7.2.If µ < 0, the convolution is defined by where K m (z) = z −m K(z) and m smallest integer such that m > −µ.
For K(z) = z −α the convolution coincides with the fractional integral of order α, i.e., according to the operational notation, we can write For β > 0, ∂ β t is equivalent to the Riemann-Liouville fractional derivative, see definition (28).
Runge-Kutta convolution quadrature has been derived in [18] and applied to (10) provides approximations at time-vectors t n = (t n,j ) s j=1 , with t n,j = t n + c j h and t n = nh, defined by where ) and the weight matrices W j are the coefficients of the power series with The notation in (11) again emphasises that the composition rule holds also after discretization: if (11) defines the approximation at the time grid t n+1 , since c s = 1.Denoting ω j (K) the last row of W j (K), the approximation reads For the rest of the paper we will denote by W j = W j (K) and ω j = ω j (K) the weights for the fractional integral case, i.e., for K(z) = z −α .
Remark 7 (Notation).We have defined the discrete convolution K(∂ h t )f for functions f .For a sequence f 0 , . . ., f N ∈ R s , we use the same notation and similarly for K(∂ h t )f (t n ) with the meaning FFT techniques based on ( 12) can be applied to compute at once all the required W j , j = 0, . . ., N , with N = ⌈T /h⌉, [16].The computational cost associated to this method is O(N log(N )).It implies precomputing and keeping in memory all weight matrices for the approximation of every [4] for details and many experiments.
The following error estimate for the approximation of ( 1) by ( 14) is given by [18,Theorem 2.2].Notice that we allow K(z) to be a map between two Banach spaces with appropriate norms denoted by • in the following.This will be needed in Section 7.

Real integral representation of the CQ weights
The convolution quadrature weights ω j can also be expressed as [23] ω for e n (λ) a function which depends on the ODE method underlying the CQ formula and an integration contour Γ which can be chosen as a Hankel contour beginning and ending in the left half of the complex plane.
Lemma 9.The weights are given by and where and e n (z) is the last row of E n (z).
Explicit formulas for E n and e n are given by ) where r is the stability function of the method and q(z) = b T (I − zO ι) −1 .
Proof.Since z −α is analytic in the whole complex plane but for the branch cut on the negative real axis, the Hankel contour Γ can be degenerated into negative real axis as in the derivation of the real inversion formula for the Laplace transform [11,Section 10.7] to obtain The expression for W n is obtained in the same way and the explicit formulas for E n and e n can be found in [23].
The next properties will be used later in Section 4 Lemma 10.There exist constants γ > 1, b > 0 and C q > 0 such that |r(z)| ≤ e γ Re z , for 0 ≤ Re z ≤ b, and q(z) ≤ C q , for Re z ≤ b, where C q depends on the choice of the norm • .
Proof.Fix a b > 0 such that all the poles of r(z) (and hence q(z)) belong to Re z > b.Define now where we have used the properties of r(z) to see that sup Re z=0 Recall that q(z) = b T (I − zO ι) −1 .As all the singularities of q are in the half-plane Re z > b and q(z) → 0 as |z| → ∞, we have that q(z) is bounded in the region Re z ≤ b.
(b) For the 2-stage Radau IIA method we have As the poles of r and q are at z = 2 ± √ 2i, we can choose any b ∈ (0, 2) and obtain the optimal γ numerically using (21).For example for b = 1, we can choose γ ≈ 1.0735.Similarly we can compute C q by computing C q = sup Re z=0 or Re z=b q(z) .
For b = 1 and the Euclidian norm we have C q ≈ 1.6429.Using the same procedure, for b = 3/2, we have γ ≈ 1.2617 and C q ≈ 3.3183.
(c) For the 3-stage Radau IIA method the poles of r(z) and q(z) belong to Re z ≥ 9 2/3 6 − 9 1/3 2 + 3 ≈ 2.681.Choosing b = 1 gives γ ≈ 1.0117 and C q ≈ 1.1803, whereas for b = 1.5 we obtain γ ≈ 1.0521 and C q ≈ 1.7954.Lemma 12.There exist constants c > 0 and x 0 > 0 such that for Re z < 0. In the rest of the Section our goal is to derive a good quadrature for the approximation of ω n and W n .We will perform the same steps as in Section 2 for the ω n .The same quadrature rules will give essentially the same error estimates for the W n ; see Remark 21.

Efficient quadrature for the CQ weights
Analogously to the the continuous case (1), we fix n 0 , h, and T and develop an efficient quadrature for the CQ weights representation (17), for nh ∈ [(n 0 + 1)h, T ] and α ∈ (0, 1).

Truncation of the CQ weights integral representation
Again we truncate the integral and give a bound on the truncation error τ (L).
Lemma 14.With the choice L = Ah −1 , the truncation error is bounded as or more explicitly Proof.From Lemma 12 we have that ensures τ (L) ≤ tol.The estimate becomes uniform in n > n 0 by setting n = n 0 + 1 in the above error bound.
Proof.We have from above from which the result follows.
Remark 16.In practice we find that instead of using Corollary 15, better results are obtained if a simple numerical search is done to find the optimal A such that the right-hand side in (22) with n = n 0 + 1 is less than tol.To do this, we start from A = 0 and iteratively approximate the integral in (22) for increased values of A (A ← A + 0.125 in our code) until the resulting quantity is below our error tolerance.The approximation of the integrals is done by the MATLAB built-in routine integral.Notice that this has to be done only once for each RK-CQ formula and value of α ∈ (0, 1).

Gauss-Jacobi quadrature for the CQ weights
In a similar way as in Section 2.2, we consider the approximation of the integral and investigate in the first place the approximation of for some suitable L 0 > 0 by using Gauss-Jacobi quadrature.Changing variables as in Section 2.2 we obtain (y + 1) −α e n (−h(y + 1)L 0 /2)dy and apply Theorem 3 to estimate the error with the weight w(x) = (x + 1) −α and integrand Theorem 17.Let with b and γ from Lemma 10, and Otherwise we have the bound Proof.We again consider the Bernstein ellipse E ̺ around [−1, 1], but now in order to be able to use Lemma 10 and avoid the singularities of e n (z) in the right-half plane we have a restriction on ̺.Namely, the maximal value of ̺ is given by h and thus giving the expression for ̺ max from the statement of the theorem.The error estimate for Gauss-Jacobi quadrature then reads, by using Lemma 10, Proceeding as in Section 2.2 with γT in place of T we obtain the bound provided that the optimal value for ̺ is within the accepted interval otherwise we make the choice ̺ = ̺ max .
Remark 18.In all our numerical experiments, we have found that ̺ opt < ̺ max .

Gauss quadrature on increasing intervals for the CQ weights
We next split the remaining integral into the sum where The intervals are again chosen so that for some B ≥ 1, L j = (B + 1)L j−1 .
To each integral we apply standard Gauss quadrature, i.e., w(x) ≡ 1 in Theorem 3, with Q nodes and denote the corresponding error by τ n,j (Q).
with constants C q , c, x 0 from Lemmas 10 and 12 and Proof.The proof is the same as the proof of Theorem 5, we only need to combine the facts that |r(z)| ≤ 1 for Re z ≤ 0, the bound q(z) ≤ C q from Lemma 10, and the bound from Lemma 12.
Remark 20.To obtain a uniform bound for t n ∈ [t n 0 +1 , T ], we replace n by n 0 + 1 in the above bound.
Remark 21.We have developed the quadrature for the weights ω n .However, up to a small difference in constants, the same error estimates hold for the matrix weights W n .Certainly, due to Lemma 12, the truncation estimate is the same.The main estimate used in the proof of Theorem 17 is the bound on the stability function r(z) and on q(z).The additional terms in E n (z) would only contribute to the constant.Similar comment holds for Theorem 17.

Fast summation and computational cost
Now the efficient quadrature is available we explain how to use it to develop a fast algorithm for computing the corresponding discrete convolution.In order to do this, we split the convolution as where as before (f n ) ℓ = f (t n,ℓ ); see (11).The first term is computed exactly, whereas for the second we can use the quadrature.Let N Q be the total number of quadrature nodes and let (w k , x k ) denote the quadrature weights and nodes with the weights including the values of x −α k in the region L 0 to L where Gauss quadrature is used.Then our approximation of I 2 n has the form we see that Hence the convolution can be approximated as with the Q n,k satisfying the above recursion.Notice that for each k = 1, . . ., N Q , Q n,k is the RK approximation at time t n−n 0 −1 of the ODE: q = −x k q + f, q(0) = 0.
Thus, from one step to the next one we only need updating Q n,k , for k = 1, . . ., N Q , N Q being the total number of quadrature nodes.Set ε the target accuracy of the quadrature.Then, from the results in Section 4 it follows that the total computational cost is O(N N Q ) with For n ≥ 5, Corollary 15 implies L ∼ h −1 and from Theorem 17 a reasonable choice for L 0 is L 0 = 4/(eT ), which leads to Therefore, the computational complexity is O(| log ε|N log N ), whereas the storage requirement scales as

Numerical experiments
Given a tolerance tol > 0, time step h > 0, minimal index n 0 , final time T > 0, and the fractional power α ∈ (0, 1) we use the above estimates to choose the parameters in the quadrature.
In particular we choose L 0 = 4/T and L = Ah −1 with A such that the upper bound for the trunction error in (22) is less than tol/3.We set B = 3 and Note that in general this choice results in fewer integration intervals than when fixing B and setting J to the smallest integer such that L ≤ L 0 (1+B) J .Next, we set L j = L 0 (1 + B) j , for j = 0, . . ., J and let Q 0 denote the number of quadrature points in the Gauss-Jacobi quadrature on [0, L 0 ] and Q j , j = 0, . . ., J − 1, the number of Gauss quadrature points in the interval [L j , L j+1 ].We choose the smallest Q 0 so that the bound on τ GJ,n (Q 0 ) in Theorem 17 is less than tol/3; note that in all of the experiments below we had ̺ opt < ̺ max .By doing a simple numerical minimization on the bound in Theorem 19, we find the optimal Q j such that the error τ n,j (Q j ) < tol J −1 /3.With this choice of parameters each weight ω j , j > n 0 , is computed to accuracy less than tol.
In Figure 1 we show the error ω n − ω n , where ω n is the nth weight computed using the new quadrature scheme and ω n is an accurate approximation of the weight computed by standard means.We see that the error is bounded by the tolerance and that for the initial weights the error is close to this bound.The error for larger n is considerably smaller than the required tolerance.This is expected, as in Corollary 15 we need to use the worst case n = n 0 + 1 to determine the trunction parameter A.
We also investigate the number of quadrature points in dependence on h, T in Table 1 and on α and tol in Table 2.We observe only a moderate increase with decreasing h, tol and increasing T .The dependence on α is mild.

Fractional integral
Let us now consider the evaluation of a fractional integral where First we investigate the behaviour of the standard implementation of CQ, based on FFT.In the particular case of the two-stage Radau IIA, p = 3 and q = 2, from Theorem 8 we would expect full order convergence.We set T = 128, and h = 2 3−j , j = 0, . . ., 7, α = 1/4.We do not have access to the exact solution u(t), so its role is taken by an accurate numerical approximation.
In Figure 2 we show the convergence of the error max n |u(t n ) − u n | using the standard implementation of CQ.We compare it with the theoretical reference curve 10 −2.5 (h 3 + |log(h)| h 3+α ), which fits the results better in this pre-asymptotic regime than the dominant term h 3 on its own.Next, we apply our new quadrature implementation of CQ.We set tol = 10 −6 , and the rest of the parameters as in the above Section.We denote by ũn the new approximation of u n and plot the error |u n − ũn | in Figure 3.We see that the error is bounded by 10 −6 for all n, showing that the final perturbation error introduced by our approximation of the CQ weights remains bounded with respect to the target accuracy in our quadrature, cf.[23].We also compare computational times in Figure 3.For the implementation of the standard CQ we have used the O(N log N ) FFT based algorithm from [16].We see that for larger time-steps the FFT method is faster due to a certain overhead in constructing the quadrature points for the new method.For smaller time steps however the new method is even marginally faster.The main advantage of the new method is the O(log N ) amount of memory required compared to O(N ) amount of memory by the standard method.For example, in this computation with the smallest time step, there are N = 2048 time steps and the total number of quadrature points is 37.As each quadrature point carries approximately the same amount of memory as one directly computed time-step, we see that the memory requirement is around 50 times smaller with the new method for this example.Such a difference in memory requirements becomes of crucial importance when faced with non-scalar examples coming from discretizations of PDE.The next section considers this case.
Remark 22.For simplicity we avoid the integer case β = 1 as it is just the standard heat equation and in some places this case would have to be treated slighlty differently.
The application of CQ based on BDF2 to integrate (28) in time has been analyzed in [8,Section 8].A related problem with a fractional power of the Laplacian has been studied in [20], but not with a CQ time discretization.Here we apply Runge-Kutta based CQ.The analysis of the application of RK based CQ to (27) is not available in the literature, hence we give the analysis here for sufficiently smooth and compatible right-hand side f .We first analyze the error of the spatial discretization.

Space-time discretization of the FPDE: error estimates
Let X ∆x ⊂ H 1 0 (Ω) be a finite element space of piecewise linear functions and let ∆x be the meshwidth.Applying the Galerkin method in space we obtain a system of fractional differential equations: Find u(t) ∈ X ∆x such that Theorem 23.Let f ∈ C m ([0, T ]; L 2 (Ω)) with f (k) (0) = 0, k = 0, . . ., m − 1 and let u(t) be the solution of (29) and u(t) the solution of (27).Then if m > β we have If further m > 2β we have Proof.Consider the Laplace transform of ( 27) for some fixed δ > 0, and the bilinear form Hence, a(u, v) = Ω f vdx is the weak form of (30).The bilinear form is continuous L 2 (Ω) and using the Poincaré inequality we obtain coercivity in H 1 0 (Ω).Lax-Milgram gives us that there exists a unique û ∈ H 1 0 (Ω) solution of (30) and that Dividing by |z| β û L 2 (Ω) and using the fact that For the finite element solution, Céa's lemma gives that where û denotes the Laplace transform of u.Using the Aubin-Nitche trick, we obtain the estimate in the weaker norm where ϕ g is the solution of the dual problem Recalling that Ω is assumed to be convex, we can use standard elliptic regularity results together with where the final inequality follows from (31).Similarly and using standard approximation results we have that The proof is completed by applying Parseval's theorem.
The fully discrete system is now obtained by simply discretizing the fractional derivative at stage level using RK-CQ: for n = 1, . . ., N − 1, v ∈ X h .
Theorem 24.Let an A-stable, s-stage Runge-Kutta method of order p and stage order q be given which satisifes the assumptions of Section for any δ > 0, following from (31); see also [21] and [8].The same estimate holds for the solution operator K ∆x : f → u of the Galerkin discretization in the space X ∆x .Note also that f.The second to last equality above follows from standard properties of convolution quadrature, see for example [17,Section 4] and [22,Chapter 9], whereas the last one is simply the use of operational notation explained in Remark 7. The result now follows from Theorem 8, Theorem 23, and the triangle inequality.
or using the definition of V n At each time step this system needs to be solved, where the expensive part is the computation of the discrete convolution in the right-hand side and the storage of all the vectors V j .This problem is resolved by our fast method of evaluation of discrete convolutions with the following variation with respect to Section 5 in order to deal with the stages: with satisfying the recursion As a final point let us note that due to ( 12) and hence As the spectrum of O ι −β is in the right-half complex plane the problem to be solved at each time-step has a unique solution.
For the numerical experiments we let Ω be the square with corners (−1, −1) and (1, 1) and choose f so that the exact solution is u(x, t) = sin 3 3  2 πt cos 1 2 πx 1 cos 1 2 πx 2 .We let the final time be T = 7, fix the finite element space on a triangular mesh with meshwidth ∆x = 5 × 10 −3 and compute the error in the L 2 (Ω) norm at t = T .The error and memory requirements as the number of time-steps is increased are given in Table 3 for our new method and for the standard implementation of the CQ.We have used as tolerance tol = 10 −4 and the 2-stage RadauIIA based CQ, for which the theory predicts convergence of order O(h 3 ).We see that the error is the same for the two implementations of the CQ, achieving in both cases the predicted order 3, but that the memory requirements for the new method stay almost constant whereas for the standard implementation they grow linearly.

Table 1 :
Dependence of the total number of quadrature points on time step h and final time T .The other parameters are fixed at n 0 = 5, B = 3, tol = 10 −6 , α = 0.5.On the left the data is for backward Euler and on the right for the 2-stage Radau IIA CQ.

Figure 1 :Figure 2 :
Figure 1: We show the error in the computation of the 2-stage Radau IIA weights ω n for n > n 0 with two different tolerances.The number of quadrature points is also shown.The results are for α = 0.5, T = 5, h = 10 −2 , n 0 = 5, and B = 3.

Figure 3 :
Figure 3:  We plot in the left graph the difference |ũ n − u n | against t n , where u n is computed using the standard implementation of CQ and ũn with the new method at tol = 10 −6 .On the right we plot the required time for the two methods.

Table 3 :
We show the error and the memory requirements for the new method and the standard implementation of CQ.