Pseudospectral roaming contour integral methods for convection-diffusion equations

We generalize ideas in the recent literature and develop new ones in order to propose a general class of contour integral methods for linear convection-diffusion PDEs and in particular for those arising in finance. These methods aim to provide a numerical approximation of the solution by computing its inverse Laplace transform. The choice of the integration contour is determined by the computation of a few suitably weighted pseudo-spectral level sets of the leading operator of the equation. Parabolic and hyperbolic profiles proposed in the literature are investigated and compared to the elliptic contour originally proposed by Guglielmi, L\'opez-Fern\'andez and Nino. In summary, the article (i) provides a comparison among three different integration profiles; (ii) proposes a new fast pseudospectral roaming method; (iii) optimizes the selection of time windows on which one may arbitrarily approximate the solution by no extra computational cost with respect to the case of a fixed time instant; (iv) focuses extensively on computational aspects and it is the reference of the MATLAB code https://github.com/MattiaManucci/Contour_Integral_Methods.git, where all algorithms described here are implemented.


Introduction
We consider convection diffusion PDEs of the form U (x, 0) = U 0 (x) with A a linear second order elliptic operator. After discretizing the problem in space, we study efficient numerical integrators for the Cauchy probleṁ with A representing a suitable discretization of the elliptic operator A and b is a source term which possibly includes boundary contributions. We are particularly interested in equations arising in mathematical finance, such as Black-Scholes, Heston or Heston-Hull-White equations [3,7,9], but our approach is by no means restricted to them.
In order to approximate the solution u(t) to (2) one may use Runge-Kutta methods, multistep integrators as well as splitting schemes. The drawback of these time-stepping schemes is that in order to approximate the solution at a certain time T = t n , one needs to compute an approximation of the solution at grid points 0 < t 1 < t 2 < . . . < t n , which would be particularly demanding if T is large. As an alternative, it is possible to derive methods based on the Laplace transform and its numerical inversion, which do not advance on a grid. In the literature this approach has been widely studied for pure diffusion equations (see e.g. [5,18,19,22]) and for convection diffusion equations recently in [6]. An important case is when the time T at which one is interested to determine the solution is not known exactly but is uncertain although it belongs to a certain time window of moderate size. In such case it would be convenient to develop methods which do not require substantial additional computations with respect to the model case when T is fixed a priori. This is another goal of this article, i.e. to discuss and analyze methods able to approximate the solution on suitable time windows.
In the sequel of the article -when not indicated differently -the considered norm is the spectral one, that is the matrix norm induced by the vector Euclidean norm.
The magnitude of the resolvent norm (zI − A) −1 has a crucial role in the rate of convergence of any contour integral method based on Laplace transformation. Due to this, the choice and parametrization of the integration contour is of main importance. In a recent paper [6], an elliptic profile has been proposed, in connection to the knowledge of the ε-pseudospectrum of A (see [24]) for suitable values ε > 0. Since A is in general non-normal, due to the convection terms in the operator A, the pseudospectrum may increase fast around the spectrum of A, making the problem challenging. We assume the existence of the Laplace transform of b and that it admits a bounded analytic extension to a suitable region of the complex plane outside the spectrum of A. We then apply the Laplace transform to (2), which yields, for the Laplace transform of u,û = L(u), whereb = L(b) and I stands for the identity matrix.
After solving (4), we reobtain u by considering the inverse Laplace transform being the contour Γ an open piecewise smooth curve running from −i∞ to +i∞ surrounding all singularities ofû. To approximate the Bromwich integral (5), we parameterize the integration contour Γ by z = z(x), x ∈ R, for a suitable mapping z(x), so that with G appropriately defined. Since we are interested in approximating u(t) within precision tol, we will only consider the portion of the Bromwich integral parameterized in [−cπ, cπ], this is, for a certain truncation parameter c ∈ (0, c max ), which we determine by the estimate |G(cπ)| = tol for tol the desired accuracy. Finally, the application of a quadrature formula to approximate (5) provides a numerical approximation of u, for a given time t, or even time windows of the form [t 0 , Λt 0 ], Λ > 1, without need of computing it at intermediate time instants. An application of the trapezoidal rule provides the desired approximation I N of I. Note that the computation of each term in the summation (6) requires the solution of a shifted linear system A − z(ξ j )I. An advantage of the method we propose is that these computations can be done in parallel. Furthermore, if the integrand is conjugate symmetric, the number of addends and thus, the number of linear systems, can be halved. This is often the case in applications, such as the ones we consider here.
Assuming that the Laplace transform can be analytically extended to the left half of the complex plane and that this extension is properly bounded with respect to z, several authors have proposed different contour profiles and parametrizations for Γ. We refer the reader to the recent article [6] for a detailed review of the literature concerning the crucial choice of the profile Γ. In this paper we extend the results of [6] by considering not only elliptic but also parabolic and hyperbolic profiles, which we compare. The parametrization of all contours is optimized by using the knowledge of the pseudospectrum of A on a region of the complex plane surrounding the spectrum of A. A novel Newton iteration is developed to obtain the required knowledge of pseudospectral level sets. In this way we are able to determine more specifically, accurately and efficiently the required pseudospectral level curves and avoid the use of the software eigtool [27] as it was done in [6]. We notice that since the exponential factor in (5) reduces the norm of the integrand function when Re(z) is sufficiently large and negative, we have to control the pseudospectrum of A only in a vertical strip of the complex plane. A main advantage of the method we discuss is that it provides an approximation of the solution by a prescribed accuracy tol, simply increasing the number of quadrature points on the integration contour, without changing the integration profile and taking advantage of previous computations. The paper is organized as follows. In Section 2 we describe the three contours we consider (elliptic, parabolic and hyperbolic). In Section 3 we present a new method to obtain approximations of pseudospectral level sets, which does not require making use of Eigtool. In Section 4 we study in full detail the pseudospectra of the 1D Black and Scholes operator. In Section 5 we provide the determination of the parameters characterizing the whole procedure based on sharp error estimates. In Section 6 we compare the profiles and present some numerical illustrations. Finally, in Section 7 we focus our attention on implementation issues and present a Matlab code aimed to approximate the solution of the problem by means of any of the considered methods.

The integration contours
We propose contours Γ in (5) which are either elliptic, parabolic or hyperbolic arcs, possibly linked to half-lines.

Elliptic profile: a review from [6]
In [6], Γ is parameterized by where, for constant parameters A 1 , A 2 , A 3 to be determined, parametrizes an elliptic arc and parametrize two half-lines, see Figure 2.1. The choice of the parameters A 1 , A 2 ad A 3 is discussed Figure 1: The integration profile Γ from [6].
on [6] and is fundamentally based on the knowledge of the ε-pseudospectrum of A in a rectangular region surrounding the rightmost section of the spectrum.
The only section of the integration contour Γ (7) that it is actually used in practise is the arc of ellipse parameterized by z. In order to improve the performance of the quadrature, z is extended to a rectangle in the complex plane by for a certain parameter a > 0 to be determined. We require (8) to be holomorphic in the rectangle and thus impose the Cauchy-Riemann equations. In this way we obtain that A 3 is necessarily a constant, with a 1 and a 2 real constants. The resulting mapping turns out to be entire. The holomorphy of z in the rectangle leads to the exponential convergence of the trapezoidal rule when it is applied to the integral resulting after parametrizing the elliptic contour by z(x), being the rate of convergence increasing with a, see [6]. The rectangle R is mapped into an elliptic ring-shaped region. In particular the upper horizontal side of the rectangle is mapped into the inner ellipse Γ lef t (black) in Figure 2 and is selected in a way that its rightmost section (continuous line) is external to the ε-pseudospectrum for a suitable value ε. In order to select Γ lef t , we fix the center of the ellipse z L , its right intersection with the real axis z R and one interpolation point z B . In particular z L is such that e z L t ≈ eps (the working precision) and z R is the rightmost intersection point of the ε-pseudospectrum of A and the real axis. The interpolation point z B = d + ir is chosen in such a way that the ellipse encloses the ε-pseudospectrum of A for a suitable ε, as well as the possible singularities ofb. Next the half-elliptic integration profile is determined, with coefficients a 1 , a 2 , A 3 depending on the unique free parameter a (see [6] for the details). Indeed, imposing the ellipse Γ lef t to be centered at z L and to pass through the points z R and z B = d + ir, we get where Solving (12) for a 1 , a 2 , A 3 yields A 3 = z L which only depend on the real parameter a.

Parabolic profile
For y fixed, the mapping defines a parabola symmetric with respect to the real axis in the complex plane. In order to obtain a holomorphic parametrization of the parabola, we impose the Cauchy-Riemann equations to determine A 1 and A 2 . This yields with real constants a 1 and a 2 . Since the parabolic profile is symmetric with respect to the real axis we only need of two points in order to determine it uniquely. We proceed by constructing the conformal mapping for a certain positive a. Proceeding analogously as for the elliptic profile, we call Γ lef t = z(R + ia), this is, the parabola limiting to the left the image by z of the horizontal strip | Im y| ≤ a. We impose the vertex of Γ lef t to be the rightmost intersection point of the ε-pseudospectrum of A, that we call z R , and a further interpolation control point z B = d + ir. In this way we obtain the uniparametric family of values depending on the free parameter a, which is the band width of the analyticity domain of z.

Hyperbolic profile
As it was done in [18], we set w = x + iy and notice that the function of w − sin(α + iw) = φ + iψ, with φ, ψ ∈ R, maps the horizontal line Im w = y into the left branch of the hyperbola φ sin(α − y) 2 − ψ cos(α − y) 2 = 1, whose asymptotes are given by φ + iψ with In order to better control the position of the center, the length of the horizontal semi axis and the angle of the asymptotes, we start here from the more general expression After imposing the Cauchy-Riemann equations we actually obtain that A 1 , A 2 and A 3 must be constant with respect to y, leaving a map of the form To obtain optimal error estimates, we need to control the image of the horizontal strip | Im y| ≤ a under z, which is the region in the complex planed limited by the two branches of hyperbola z(x−ia) (the one closest and furthest to the left) and z(x + ia), x ∈ R. Similarly to what we do for the other type of contours, we prescribe z(x − ia), x ∈ R, to be an appropriate critical hyperbola Γ lef t , with vertex at z R , center at z C and passing through a third point d + ir. In this way we obtain Equations (22), (21) and (20) give a uniparametric family of solutions for the parameters a 1 , a 2 and a 3 , respectively, depending on the width a of the strip of analyticity of the mapping z. We notice that the orientation we need to invert the Laplace transform is actually the opposite to the one in (19), as x runs from −∞ to ∞. This is resolved by simply taking the conjugate of (19) as parametrization.

Roaming pseudospectral sets
We consider here the case of a parabolic contour. The idea can be extended in a straightforward way to the elliptic and the hyperbolic contour. We start from an initial internal parabola uniquely identified by a prescribed interpolation point w = d + ir (together with the vertex z R , which we consider fixed). Since the parabola is symmetric with respect to the real axis, we only work with its upper half, with positive imaginary part, and consider a set of M points z k = φ k + iψ k , k = 1, . . . , M, on this curve, increasingly ordered according to their real parts φ k . The parametric form of the inner parabola is Setting as in the previous section z(x) = φ + iψ, with ψ > 0, this means that if we fix the abscissa which uniquely defines the argument of the parametrization x > 0, and consequently which depends on r and d. We easily obtain We discuss two possible approaches and after numerical simulations we have adopted the second in our code. The idea is that of modifying the interpolation point w = d + ir, which determines the parabola (together with the vertex z R , which we consider fixed) by using variational results for simple singular values. In the sequel we shall make use of the following classical result on the derivative of a simple eigenvalue (see e.g. [17,8]).
be a smooth (with respect to t) singular value decomposition of the matrix D(t) and σ(t) be a certain singular value of D(t) converging to a simple singular valueσ = 0 of D 0 = D(t 0 ).
Ifû,v are the associated left and right singular vectors, respectively, the function σ(t) is differen- At t = t 0 the left and right eigenvectors associated toσ 2 coincide and are equal toû, having unit norm. Note thatû is a certain column of U (t 0 ) determined by the position ofσ 2 in the diagonal matrix Σ(t 0 ) 2 .

An optimal although expensive approach
Differently from what proposed in [6], where both pseudospectral and a weighted version of pseudospectral level sets were considered, here we will only consider the weighted ones. This choice is motivated by the fact that we do not relay anymore on Eigtool for the computation of pseudospectral level sets. Instead, we will directly compute the value of a weighted version of the pseudospectrum, which is the one that really matters in the error estimate, see [6,Theorem 2]. More precisely, we define the weighted ε-pseudospectrum as the set which is equivalent to with σ min zI − A denoting the smallest singular value of zI − A. In particular, we are interested on the boundary of the weighted ε-pseudospectral level set, this is Note that for t = 0 we recover the standard definition of pseudospectrum and ε-pseudospectrum.
For a fixed target ε, we impose the set in (28) to lay internal to the parabola. A natural approach to achieve this goal is defining the functional with z k depending on d and r, and (τ ) + = max{τ, 0}. In this way, values ofσ k (d, r) which are larger than ε do not contribute to the functional. The goal is to compute a solution (d, r) to where -for numerical convenience -the constraint may be treated as a penalization term. This can be done by computing the gradient of F, where we use Lemma 1. The gradient is continuous and has the form with x k such that z(x k ) = z k . Note that Re(z k ) does not depend on d and r, see (23). We add to the functional the penalization term , whose gradient is obtained in a straightforward way. Then to compute a solution we can apply any gradient based method for unconstrained optimization to the functional for a sufficiently large c (in the context of a penalization methodology). The method turns out to be effective but appears to be computationally expensive due to the fact that at every step of the gradient descent method we have to compute several singular values and the associated singular vectors.

Selecting points internal to the weighted ε-pseudospectrum
A cheaper (and preferred) alternative to the previous method, is obtained by treating the points {z k } singularly, one after the other. In this way, we start by considering the first point z 1 = φ 1 +iψ 1 and compute , according to the definition in (29). Then we check the difference δε = ε −σ 1 . If δε ≤ 0 we proceed by considering z 2 and repeating the same steps, otherwise it means that we are inside the weighted ε-pseudospectrum and therefore we need to update the internal parabola. If we do not find any z k such thatσ k < ε, we may consider the point z ∈ {z k } M k=1 for which the correspondingσ is minimal and proceed in order to find a closer parabola to the weighted ε-pseudospectral level set.
The algorithm we adopt tunes the interpolation point w = d + ir, so that the updated curve is external weighted ε-pseudospectrum of A at z. We indicate by p = p(d, r) a selected point z k = z(x k ) = φ k + iψ k laying in the wrong weighted pseudospectral level set and writeσ k as In principle we want to solve the equatioñ It seems natural to fix the parameter d as the mean of the abscissas of the support points z k , i.e.
Applying Lemma 1 toσ(d, r)-with u and v left and right associated singular vectors -we get In order to accurately compute r such thatσ(d, r) = ε we make a few (say m) Newton iterations with u and v singular vectors associated to σ min (A − p(d, r )I) and r the actual ordinate of the interpolation point w.
Then we compute a new parabola, which interpolates d + ir m , reparametrize it and compute a new set of points. Iterating a few times this procedure we compute the desired parabolic profile.
Remark 2. Since we have 2 free real parameters to determine, d and r, we may consider at the same time two points z 1 and z 2 to which correspond values ofσ(d, r) smaller than the target value ε. This would provide a simple variant to the method described above. We would first determine two points z 1 and z 2 such thatσ j < ε, for j = 1, 2 and then solve equations (in analogy to (33)) with respect to r and d by Newton method. It would be natural to expect that this method would result into a fewer number of iterations.

A case study: the 1D Black and Scholes equation
The well known (deterministic) Black-Scholes equation [3] has the following form: for L, t given, where the unknown function u(s, τ ) stands for the fair price of the option when the corresponding asset price at time t − τ is s and t is the maturity time of the option. Moreover, r ≥ 0, σ > 0 are given constants (representing the interest rate and the volatility, respectively). In practice, for the sake of numerical approximation, we consider a bounded spatial domain, setting L < s < S for a sufficiently large S. We take (35) together with the following conditions, typical for the European call option, cf. [12]: being K the reference strike price. In this Section we extend the theory developed in [21] for the 1D convection-diffusion operator to equation (35). In this way we are able to theoretically determine a region in the complex plane where the norm of the resolvent of the Black-Scholes differential operator grows exponentially. This knowledge allows us to use (35) as a benchmark problem to test the new pseudospectral roaming strategy in Section 3. Our goal is to solve (35) with (36) by applying the Laplace transform method. To do this we first transform the problem to an equivalent one with homogeneous boundary conditions. This is easily achieved by considering v(s, τ ) = u(s, τ ) − y(s, τ ), The differential equation for v reads with initial and boundary data We can now apply the Laplace transform to both sides of (37) with respect to τ . This leads to the following equation for V (s, z), the Laplace transform of v(s, τ ): with L the differential operator for the Black-Scholes problem with homogeneous boundary conditions.

Pseudospectra of the Black-Scholes equation
For our analysis we set L = 1, which is reasonable if S, K >> 1 in (36) and allows to apply the change of coordinates x = log(s) while keeping the domain bounded. After this change of variable we obtain the evolution problem u t = Lu, with a second order diffusion-convection-reaction differential operator with constant coefficients on a bounded domain with homogeneous boundary conditions of Dirichlet type. We thus can compute explicitly the eigenvalues and eigenfunctions of L by applying it to a mapping of the form ϕ(x) = e αx . In this way, we obtain and Then, for each λ real we have two associated values of α, namely For any λ and corresponding α + and α − , the function satisfies (40) in the interior of [0, log(S)] and the boundary condition at x = 0. It also satisfies the boundary condition at x = log(S) provided e α + log(S) = e α − log(S) , that is, (α + − α − ) log(S) = 2πni for some nonzero n ∈ Z. By (42), this amounts to the condition (log(S)/ν) (r + ν) 2 + 4λν = 2πni, and upon squaring we obtain the following eigenvalues The corresponding region in the λ-plane is the image of B under the function λ = να 2 + (r − ν)α − r, which we denote by Π: The "critical parabola" that bounds Π is the image of the boundary of B under the same function, which we can simply represent by since Re(α) = − (r−ν) ν maps onto the same parabola as Re(α) = 0. Suppose now that λ is any complex number in the interior of Π so that Re(α + ) < 0 and Re(α − ) < 0. Then φ(x) decreases exponentially with x, so if log(S) is reasonably large, the boundary condition u(log(S)) = 0 is nearly satisfied, with an error of order e µ log(S) = S µ , where µ = max{Re(α + ), Re(α − )}. Thus φ(x) is nearly an eigenfunction of L, though λ may be far from any of the exact eigenvalues. Then we can just repeat the arguments and passages of [21] to get their results. The main difference is that, in our case, we also consider a reaction term, anyway results of [21] for the convection-diffusion operator can be extended with their same procedure to the convection-diffusion-reaction operator. We now state our version of Theorem 5 of [21] which deal with the Black-Scholes differential operator.
This result tells us that the resolvent norm changes exponentially along any vertical line inside Π. Indeed we know, by construction, that the critical parabola is the curve such that Re(α) = 0 and therefore S −µ = 1, while on the real axis, for Re(λ) sufficiently small the real part of α − and α + are the same and equal to − r−ν 2ν , that corresponds to the case where S −µ is maximized.

Symmetrizability and a further estimate
As done in [21] we can explicitly symmetrize the BS operator L. First let's define ρ = r−ν 2ν and u(x) = e −ρx v(x), which implies From here, we follow the analysis in [21] and derive the following bound for the resolvent norm of Black-Scholes operator.
Theorem 4. For any d > 0, r > ν and λ ∈ C, The discussion done in the previous subsection this theorem suggested that the pseudospectra level curves of L are bounded approximately by parabolas. On the light of Theorem 4 we can say that the exponential bound by parabola does not hold as |λ| → 0 but, for any fixed ε and S, Λ ε (L) the ε-pseudospectra is contained in a strip of finite (though typically large) width:

Numerical validation
We have to keep in mind that the results previously exposed hold for the continuous operator. Since our aim is to have a numerical validation we have to deal with the discrete version L h . The way properties of the pseudospectrum of the discrete operator converge to the properties of the continuous one is not treated here and, to our knowledge, is an open research problem. Nevertheless it is reasonable to expect for small h a behaviour close to the one exhibited by the differential operator. Indeed this is what we observe in Figure 3 where we set S = 200 and used 2000 points for the space discretization; the resolvent norm is plotted in a subset of the Π region. The comparison between the magnitude of the resolvent norm estimate (48) and the resolvent norm of the discrete operator indicates a very similar behaviour of the two operators. We clearly see that both decrease when approaching the critical parabola and are maximal close to the real axis.

Choice of the parameters
In [6] error estimates are developed to provide a practical strategy to optimize the elliptical integration contour and minimize the required number N of quadrature nodes, for a prescribed target accuracy. However, the main results in [6], namely Theorems 1 and 2, do not depend on the specific choice of an ellipse and apply in a straight forward way to the parabolic and hyperbolic contours we have described in Section 2. The steps to determine the integration contour are thus common for the three types of contours under study and are the following: 1. Compute z L from e z L t = ε, with ε the working precision.
2. Compute the critical curve Γ lef t according to Section 3. This procedure provides an explicit parametrization of Γ lef t , of the form ψ(x), x ∈ R.
3. Compute c max as the unique value that satisfies Re(Γ(c max π)) = z L . This gives, according to the discussion in Section 2: for the ellipse, for the parabola, , for the hyperbola.
(53) 4. Compute a by following the same steps as in [6,Section 3.2]. This requires the identification of the right-most point of the external ellipse/parabola/hyperbola, this is , for the ellipse, for the parabola, where r in the first expression above is the imaginary part of the control point d + ir from Section 2 and θ = arccos d−z L z R −z L , as in (13) and (14). In order to have a more robust estimate we take into account the constant M lef t , named M + in [6, equation (23)], which is defined as whereΓ lef t is a suitable restriction of Γ lef t .
We also consider a different estimate for M right , named M − in [6, equation (24)], which is defined as whereΓ right is a suitable restriction of Γ right . Note that M lef t can be bounded as and we can take advantage of the computation done in step 2 for the weighted ε-pseudospectrum to approximateM lef t . Concerning M right we first assume that the maximum is reached on z v (a), the vertex ofΓ right ; in this way we have Algorithm 1 Numerical algorithm for computing a.
Input: a max , Γ lef t ,M lef t ,M right = 1, a 1 = a max , err = 1 and j = 0 Output: a 1: while err ≥ prec do 2: Find a j+1 ∈ [0, a max ] such that a j+1 = arg min f (a), see (60); 3: if err ≥ prec; then The value of a is then computed by minimizing the right hand side of (59) which requires an interval of the form a ∈ [a min , a max ], with prescribed bounds a min and a max . Clearly a min has to be set equal to 0 while to determine a max we need a further discussion. The procedure described above does not take into account the numerical error due to the conditioning of zI − A, whose effect is also amplified by the multiplication with the exponential term. This contribution becomes relevant when one wants to reach high accuracies or has to deal with ill conditioned systems. To keep under control this error we estimate it on the vertex of the integration profile defined by a = a max . If this is higher than the required tolerance we reduce a max and we check again until the estimate is below the required tolerance. Reducing a max has the effect of bringing the vertex profile closer to the imaginary axis, which reduces the amplification effect of the exponential. We remand to Subsection 6.3 for the description on how we estimate the numerical error due to the solution of linear systems. Two aspects are remarkable.
(i) It may be too restrictive to select the first acceptable a max , since we may exclude some acceptable a. Therefore, we select the value computed in the iteration before convergence and add a penalization term to (59). This penalization term is based on the estimation of the numerical error due to the conditioning of zI − A on the vertex of the integration profile. The new function to minimize becomes Algorithm 2 Numerical algorithm for approximating c, K.
Input: K (1) given, Find c (j) such that Re(z(c (j) π)) = 1 t log tol K (j) ; 3: where w is a positive scalar and (ii) We note that evaluating (58) for every a ∈ [a min , a max ] is computationally expensive due to the presence of the resolvent norm. We thus implement the iteration procedure described in Algorithm 1, which we have observed to convergence in very few iterations (from 2 to 6 depending on which type of contour we use) in all the numerical tests done.

Compute a truncation parameter c ≤ c max . The actual numerical integration is performed for
x ∈ [−cπ, cπ]. The computation of c requires an iterative procedure resumed in Algorithm 2 and follows precisely [6, Section 3.5], for the three types of integration contours.

Comparison of the integration profiles: numerical illustrations
In this section we apply the considered contour integral methods to two illustrative problems arising from finance. The first problem is the Black-Scholes equation, while the second is the Heston equation. Both models are the same as in [6].
We show the absolute error rather than the relative error in order to check the match with the target accuracy tol, that is the accuracy we want to reach in the approximation of (5). Similarly to what has been done in other publications presented in the literature, we compute the time approximation error for a semidiscretization in space of the PDE. We do not address here specific estimates of the spatial discretization error, but rely on the referred literature for every example. We measure the error in time against a reference solution computed by using the MATLAB function expmv, see [1, We also notice that in all our tests we construct the inner curves as explained in Section 3, taking the weighted ε-pseudospectral level set with ε = 10 −7 .

Black-Scholes equation
Following the same strategy adopted in [13,6], we discretize (35) in space on a uniform space grid of N h = 2000 points in [L, S] for L = 0, S = 200, by using the classical centered finite difference scheme. The error associated to the spatial discretization is then O(∆x 2 ), with ∆x the diameter of the spatial grid, being around 0.001 for our choice of N h . This can be observed in Figure 4, where we consider a reference solution computed with N h = 2 · 10 4 spatial nodes and approximated in time by using expmv.
We set r = 0.06, σ = 0.05, and K = 80. We plot the error for a selection of tolerances for the cases t = 1 ( Figure 5) and t = 10 ( Figure 6). In Figure 7 we show the selected profiles of integration for the tolerance tol = 5 · 10 −5 at time t = 1 and t = 10. In the pictures we also highlight, for each profile, the estimated number of quadrature nodes N according to (62), to reach the target accuracy. All contour methods are effective when applied to the Black-Scholes test problem. The hyperbolic contour is slightly faster in reaching the target accuracy for t = 1, while for t = 10 the parabolic contour seems to provide the best results. In Table 1 we report the values of N estimated by (62), which are similar for the three types of contour. Also, as expected, the estimate returns larger N when a higher accuracy is required. We note that for t = 10 and tol = 5 · 10 −9 the number of estimated quadrature points increases significantly. This is due to the fact that when considering larger times, the exponential term may considerably amplify the numerical error introduced by the solution of linear systems. To avoid this, one has to select integration profiles which lie in a region with small positive real part; this results in lower value of a and thus in a larger number of quadrature nodes in order to reach the prescribed accuracy (see (59)).

Heston equation
The Heston equation [7] is given by The unknown function u(s, v, τ ) represents the price of a European option when at time t − τ the corresponding asset price is equal to s and its variance is v. We consider the equation on the unbounded domain where the time t is fixed. The parameters κ > 0, σ > 0, and ρ ∈ [−1, 1] are given. Moreover equation (63) is usually considered under the condition 2κη > σ 2 that is known as the Feller condition (see [15]). We take equation (63) together with the initial condition where K > 0 is fixed a priori (and represents the strike price of the option), and boundary condition For the numerical solution of (63), we need to choose a bounded domain of integration, we follow [10] for this issue. In particular, we fix two positive constants S, V and we let the two variables s, v vary in the set On the new boundary, we need to add two more conditions (specific for the European call option), which are treated analogously yo the boundary condition in (36). We use the spatial discretization proposed in [10] for k = 1.5, η = 0.04, σ = 0.3, ρ = −0.9, r d = 0.025, r f = 0, K = 100, L = 0, S = 8K, V = 5. In Figure 8, we plot the error for a selection of tolerances for t = 1 and in Figure 9 tol N at t = 1 N at t = 10 Ellipse Parabola Hyperbola Ellipse Parabola Hyperbola 5 · 10 −3 13   we do the same for t = 10. In Figure 10 we show the selected profiles of integration for tol = 5· 10 −5 at t = 1 and t = 10.
We obtain similar results for the Heston test problem to the ones reported for the Black-Scholes problem. The hyperbolic contour is again the fastest one in reaching the target accuracy for t = 1. For t = 10 elliptic, parabolic and hyperbolic contours show almost the same behaviour. In Table  2 we report the values of N estimated by (62) and again we do not observe remarkable differences for the three types of contour.

Extension to the case of time intervals
This situation is particularly important when the time T at which the solution is required is only known approximately. In this subsection we extend the results in [6, Section 7] to determine the solution on an entire time window [t 0 , t 1 ], with Most part of the computational effort in the evaluation of (6) is devoted to the solution of linear systems with matrices zI − A, for z the quadrature nodes. However this inversion does not involve the time t, that only appears in the exponential term. Thus, it is of high interest being able to use a unique integration contour on a whole time interval [t 0 , t 1 ]. This requires a suitable modification of our strategy to construct the integration contour. We also aim to determine an acceptable amplitude of the time window.
Concerning the profile of integration, which is uniquely defined by the parameter a, we start by constructing Γ lef t as explained in Section 3, for t = t 0 (lower end of the time interval). The choice of t 0 reflects into the setting of the parameter z L as explained in Section 5. An application of the construction described in Section 3 gives the interpolation point d+ir, which uniquely defines Γ lef t . Then we determine a by following the procedure described in Section 5 with t = t 1 , the upper end of the time interval. We use the profile determined in this way for all t ∈ [t 0 , t 1 ].
Despite the fact that, theoretically, the amplitude of the window can be arbitrarily large, we need to keep into account the role of the exponential term in amplifying the error introduced by the conditioning of zI − A. We approximate the exact solution u(t) by the linear combinatioñ whereû j =û(z(x j )) + ρ j and ρ j is the error in the numerical solution of the linear system for x j our quadrature nodes and with the assumption that the nodes x j , the parametrization z(x), and its derivative z (x) are computed exactly. The actual error in our computation is given bỹ that we can estimate in the following way: where err T is the component of the total error due to the truncation of the integral, err N is the component due to the approximation of the integral by the quadrature rule and finally err num N is the component due to the fact that we operate with finite precision arithmetic. This last component reads as Therefore, once the integration contour for a time window [t 0 , t 1 ] is fixed, we can observe that err num N does not remain constant but instead changes exponentially with respect to time. Since we want that err N ≤ tol for all t ∈ [t 0 , t 1 ], we need to check that (68) is below tol at t 0 and t 1 . If so, we proceed with the computation, otherwise we halve the amplitude of the window by increasing t 0 or decreasing t 1 , depending on which of the two time instants fails to satisfy the inequality (68) smaller than tol. If for both t 0 and t 1 (68) is greater than tol we deduce that the problem is possibly ill conditioned and we suggest to increase tol.
About the truncation parameter c we observe that this value decreases as time increases, therefore we run Algorithm 2 for t = t 0 and we use the computed value of c in t 0 for every t ∈ [t 0 , t 1 ]. Doing so we expect to have the final error proportional to tol for t = t 0 and smaller than tol when t > t 0 .
We show few numerical experiments for both Black-Scholes and Heston equations. In particular, we make the experiments on the intervals [0.1, Λ0.1] and [1, Λ1] for Λ = 10. In the plots of Figures  11, 12, 13, 14 we show the numerical results for Black-Scholes and Heston equations. The target tolerance chosen is tol = 5 · 10 −8 for Black-Scholes and tol = 5 · 10 −4 for Heston. We also fix z R = 0.06 for both problems. We note that for the Heston problem the time window [1,10] has been found to be too large by the procedure previously described, therefore the initial time was automatically increased. For the Black-Scholes problem a slow-down in convergence is observed as t increases, we see this phenomena also in the Heston problem even if in a less remarkable way. This because, as t increases, the truncation value used c is larger than what should be required to achieve the prescribed tolerance tol. This results into a smaller error for N large enough but also produces a slower convergence.

Implementation and codes
The aim of this Section is to explain the structure and the main features of the code available at [20]. We assume that a space discretization is already available and thus we have the operator A, the Laplace transform of the known termb and the initial solution u 0 . The code provided includes three test problems: a classical convection-diffusion problem, the Black and Scholes model and the Heston model.

Construction of the inner curve and integration map
The first step is the computation of the inner curve Γ lef t , as described in Section 3; then we determine the bandwidth of the map a and the truncation parameter c. Concerning Γ lef t there are some practical implementation remarks we want to mention.
First, given a discrete problem, we can use operators arising from coarser discretization to save computational effort in evaluating the pseudospectra. Not so much is known about the behavior of the pseudospectra as the size of the operator increases. It is reasonable to expect a behavior similar to the one of the condition number, i.e. that it increases as size increases. Therefore the magnitude of pseudospectra is most likely to be underestimated considering operators of smaller size. Anyway, since we have exponential convergence with respect to N , it is sufficient to add few nodes to get a final solution under the required accuracy.
Second, we consider the Newton iteration (34). At the beginning of the process it often happens that the perturbation due to the gradient is such that r +1 r . If this is the case, we do the following update: r +1 = r + pr with p a suitable positive real number chosen by the user (0.5 in our problems).
Finally let us comment the choice of the set of points {z k }. We define two grids of points, one finer than the other. We start by computing δε in the coarser grid; if δε ≤ 0 we consider the successive point in the same coarse grid. Otherwise we pass to the finer one and we consider points there until δε > 0. This is all implemented in the functions 1. Elliptic_Map, 2. Parabolic_Map, 3. Hyperbolic_Map.

Conclusions
In this article we present significant algorithmic developments of the method in [6] for the approximation of convection-diffusion problems by means of the inverse Laplace transform. The main achievements are the extension to more general quadratic contour curves, the setting of a novel method to roam pseudospectral level sets -which makes the whole algorithm faster and independent of the code Eigtool-and the extension of the method to approximate the solutions to the PDEs at time windows of suitable length to achieve a given target accuracy.