Pseudospectral Roaming Contour Integral Methods for Convection-Diffusion Equations

Abstract

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ópez-Fernández and Nino 2020, see Guglielmi et al. (Math Comput 89:1161–1191, 2020). 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 [20], where all algorithms described here are implemented.

1 Introduction

We consider convection diffusion PDEs of the form

$$\begin{aligned} \frac{\partial U}{\partial t}(x,t)= & {} {{{\mathcal {A}}}} (x) U(x,t) + f(x,t),+ \mathrm {B.C.} \nonumber \\ U(x,0)= & {} U_0(x) \end{aligned}$$

(1)

with ${{{\mathcal {A}}}}$ a linear second order elliptic operator. After discretizing the problem in space, we study efficient numerical integrators for the Cauchy problem

$$\begin{aligned} {\dot{u}} = A u + b(t), \qquad u(0)=u_0, \qquad t > 0, \end{aligned}$$

(2)

with A representing a suitable discretization of the elliptic operator ${{\mathcal {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<\ldots <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 $\left\| \left( zI-A\right) ^{-1}\right\|$ 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 $\varepsilon$-pseudospectrum of A (see [24])

$$\begin{aligned} \sigma _\varepsilon (A) = \left\{ z \in {\mathbb {C}}\ : \left\| \left( z\mathrm {I}-A\right) ^{-1}\right\| \le \frac{1}{\varepsilon }\right\} , \end{aligned}$$

(3)

for suitable values $\varepsilon > 0$. Since A is in general non-normal, due to the convection terms in the operator ${\mathcal {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, ${\hat{u}}={{{\mathcal {L}}}}(u)$,

$$\begin{aligned} {\hat{u}}(z)=\left( zI-A\right) ^{-1}\left( u_0+{\hat{b}}(z)\right) , \end{aligned}$$

(4)

where ${\hat{b}}={{{\mathcal {L}}}}(b)$ and I stands for the identity matrix.

After solving (4), we reobtain u by considering the inverse Laplace transform

$$\begin{aligned} u(t)=\frac{1}{2\pi {\mathrm {i}}}\int _{\Gamma }\mathrm {e}^{zt}{\hat{u}}(z)\,dz, \end{aligned}$$

(5)

being the contour $\Gamma$ an open piecewise smooth curve running from $-{\mathrm {i}}\infty$ to $+{\mathrm {i}}\infty$ surrounding all singularities of ${\hat{u}}$. To approximate the Bromwich integral (5), we parameterize the integration contour $\Gamma$ by $z=z(x)$, $x\in {\mathbb {R}}$, for a suitable mapping z(x), so that

$$\begin{aligned} \int _{\Gamma }\mathrm {e}^{zt}{\hat{u}}(z)\,dz = \int \limits _{{\mathbb {R}}} G(x) dx, \end{aligned}$$

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\pi ,c\pi ]$, this is,

$$\begin{aligned} I = \int \limits _{{\mathbb {R}}} G(x) dx \approx \int \limits _{-c \pi }^{c \pi } G(x) dx, \end{aligned}$$

for a certain truncation parameter $c \in (0,c_{\text {max}})$, which we determine by the estimate

$$\begin{aligned} | G(c \pi ) | = tol\end{aligned}$$

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, \Lambda t_0]$, $\Lambda > 1$, without need of computing it at intermediate time instants. An application of the trapezoidal rule

$$\begin{aligned} I_N=\frac{2 c \pi }{N}\sum _{j=1}^{N-1}G(\xi _j)\, \quad \text{ with } \ \xi _j=-c \pi + j \frac{2 c \pi }{N},\quad j=1,\ldots ,N-1. \end{aligned}$$

(6)

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(\xi _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 $\Gamma$. We refer the reader to the recent article [6] for a detailed review of the literature concerning the crucial choice of the profile $\Gamma$. 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 ${\text {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 Sect. 2 we describe the three contours we consider (elliptic, parabolic and hyperbolic). In Sect. 3 we present a new method to obtain approximations of pseudospectral level sets, which does not require making use of Eigtool. In Sect. 4 we study in full detail the pseudospectra of the 1D Black and Scholes operator. In Sect. 5 we provide the determination of the parameters characterizing the whole procedure based on sharp error estimates. In Sect. 6 we compare the profiles and present some numerical illustrations. Finally, in Sect. 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.

2 The Integration Contours

We propose contours $\Gamma$ in (5) which are either elliptic, parabolic or hyperbolic arcs, possibly linked to half-lines.

2.1 Elliptic Profile: A Review From [6]

In [6], $\Gamma$ is parameterized by

$$\begin{aligned} z(x)=\left\{ \begin{array}{ll} \ell _1(x), &{}\quad x \le -\frac{\pi }{2}, \\ z(x), &{}\quad -\frac{\pi }{2} \le x \le \frac{\pi }{2}, \\ \ell _2(x), &{}\quad x\ge \frac{\pi }{2}, \\ \end{array}\right. \end{aligned}$$

(7)

where, for constant parameters $A_1, A_2, A_3$ to be determined,

$$\begin{aligned} z(x) = A_1\cos x+{\mathrm {i}}A_2\sin x+A_3 \end{aligned}$$

parametrizes an elliptic arc and

$$\begin{aligned} \ell _{1}(x)= & {} A_3 + x+\frac{\pi }{2}-{\mathrm {i}}\left( A_2 - d\left( x+\frac{\pi }{2}\right) \right) ,\\ \ell _2(x)= & {} A_3-x+\frac{\pi }{2}+{\mathrm {i}}\left( A_2+d\left( x-\frac{\pi }{2}\right) \right) \end{aligned}$$

parametrize two half-lines, see Fig. 1.

Fig. 1

The integration profile $\Gamma$ from [6]

The choice of the parameters $A_1, A_2$ ad $A_3$ is discussed on [6] and is fundamentally based on the knowledge of the $\varepsilon$-pseudospectrum of A in a rectangular region surrounding the rightmost section of the spectrum.

The only section of the integration contour $\Gamma$ (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

$$\begin{aligned} z(x+iy)=A_1(y)\cos x+iA_2(y)\sin x+A_3(y), \quad x\in \left[ -\frac{\pi }{2},\frac{\pi }{2}\right] , y\in [-a,a], \end{aligned}$$

(8)

for a certain parameter $a>0$ to be determined. We require (8) to be holomorphic in the rectangle

$$\begin{aligned} R = [-\pi /2,\pi /2] \times [-\mathrm {i}a, \mathrm {i}a] \end{aligned}$$

and thus impose the Cauchy-Riemann equations. In this way we obtain that $A_3$ is necessarily a constant,

$$\begin{aligned} A_1(y)= & {} a_1\mathrm {e}^y+a_2\mathrm {e}^{-y}, \end{aligned}$$

(9)

$$\begin{aligned} A_2(y)= & {} a_2\mathrm {e}^{-y}-a_1\mathrm {e}^y, \end{aligned}$$

(10)

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].

Fig. 2

The ellipse $\Gamma _{left}$ (in black) and the integration profile (in blue) (Color figure online)

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 $\Gamma _{left}$ (black) in Fig. 2 and is selected in a way that its rightmost section (continuous line) is external to the $\varepsilon$-pseudospectrum for a suitable value $\varepsilon$. In order to select $\Gamma _{left}$, 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 $\mathrm{e}^{z^L t} \approx \mathrm{eps}$ (the working precision) and $z^R$ is the rightmost intersection point of the $\varepsilon$-pseudospectrum of A and the real axis. The interpolation point $z^B=d + \mathrm {i}r$ is chosen in such a way that the ellipse encloses the $\varepsilon$-pseudospectrum of A for a suitable $\varepsilon$, as well as the possible singularities of ${\hat{b}}$. Next the half-elliptic integration profile

$$\begin{aligned} \Gamma : \quad z(x)=(a_1+a_2)\cos x+{\mathrm {i}}(a_2-a_1)\sin x+A_3, \end{aligned}$$

(11)

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 $\Gamma _{left}$ to be centered at $z^L$ and to pass through the points $z^R$ and $z^B=d+{\mathrm {i}}r$, we get

$$\begin{aligned} a_1\mathrm {e}^a+a_2\mathrm {e}^{-a}=z^R-z^L,\quad a_2\mathrm {e}^{-a}-a_1\mathrm {e}^a=\frac{r}{\sin (\theta )},\quad A_3=z^L, \end{aligned}$$

(12)

where

$$\begin{aligned} \theta =\arccos \left( \frac{d-z^L}{z^R-z^L}\right) . \end{aligned}$$

Solving (12) for $a_1,a_2, A_3$ yields

$$\begin{aligned} a_1= & {} \frac{\mathrm {e}^{-a}}{2}\left( z^R-z^L-\frac{r}{\sin (\theta )}\right) \end{aligned}$$

(13)

$$\begin{aligned} a_2= & {} \frac{\mathrm {e}^{ a}}{2}\left( z^R-z^L+\frac{r}{\sin (\theta )}\right) \\ A_3= & {} z^L \nonumber \end{aligned}$$

(14)

which only depend on the real parameter a.

2.2 Parabolic Profile

For y fixed, the mapping

$$\begin{aligned} z(x + \mathrm {i}y) = - x^2 - 2 \mathrm {i}x A_2(y) + A_1(y), \quad x\in {\mathbb {R}}, \end{aligned}$$

(15)

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

$$\begin{aligned} A_1(y)= & {} y^2 + 2 a_1 y + a_2, \\ A_2(y)= & {} y + a_1, \end{aligned}$$

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

$$\begin{aligned} z:{\mathbb {R}}\times [-{\mathrm {i}}a,{\mathrm {i}}a]\rightarrow {\mathbb {C}} \end{aligned}$$

for a certain positive a. Proceeding analogously as for the elliptic profile, we call $\Gamma _{left} = z({\mathbb {R}}+ {\mathrm {i}}a)$, this is, the parabola limiting to the left the image by z of the horizontal strip $|{\text {Im}}y| \le a$. We impose the vertex of $\Gamma _{left}$ to be the rightmost intersection point of the $\varepsilon$-pseudospectrum of A, that we call $z^R$, and a further interpolation control point $z^B = d + \mathrm {i}r$. In this way we obtain the uniparametric family of values

$$\begin{aligned} a_1= & {} -\frac{r}{2 \sqrt{ z^R - d}} - a, \end{aligned}$$

(16)

$$\begin{aligned} a_2= & {} z^R - a^2 - 2 a a_1 , \end{aligned}$$

(17)

depending on the free parameter a, which is the band width of the analyticity domain of z.

2.3 Hyperbolic Profile

As it was done in [18], we set $w=x+iy$ and notice that the function of w

$$\begin{aligned} -\sin (\alpha +{\mathrm {i}}w) = \phi +{\mathrm {i}}\psi , \quad \text{ with } \ \phi , \psi \in {\mathbb {R}}, \end{aligned}$$

maps the horizontal line ${\text {Im}}w=y$ into the left branch of the hyperbola

$$\begin{aligned} \left( \frac{\phi }{\sin (\alpha -y)} \right) ^2- \left( \frac{\psi }{\cos (\alpha -y)} \right) ^2 =1, \end{aligned}$$

whose asymptotes are given by $\phi + {\mathrm {i}}\psi$ with

$$\begin{aligned} \psi = \pm \cot (\alpha -y) \phi . \end{aligned}$$

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

$$\begin{aligned} z(x+{\mathrm {i}}y) = A_3(y) - A_2(y)\sin (A_1(y) -y + {\mathrm {i}}x ). \end{aligned}$$

(18)

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

$$\begin{aligned} z(x+{\mathrm {i}}y) = a_3 - a_2 \sin (a_1-y)\cosh x - {\mathrm {i}}a_2 \cos (a_1-y) \sinh x. \end{aligned}$$

(19)

To obtain optimal error estimates, we need to control the image of the horizontal strip $|{\text {Im}}y|\le a$ under z, which is the region in the complex planed limited by the two branches of hyperbola $z(x-{\mathrm {i}}a)$ (the one closest and furthest to the left) and $z(x+{\mathrm {i}}a)$, $x\in {\mathbb {R}}$. Similarly to what we do for the other type of contours, we prescribe $z(x-{\mathrm {i}}a)$, $x\in {\mathbb {R}}$, to be an appropriate critical hyperbola $\Gamma _{left}$, with vertex at $z^R$, center at $z^C$ and passing through a third point $d+{\mathrm {i}}r$. In this way we obtain

$$\begin{aligned} a_3&= z^C \end{aligned}$$

(20)

$$\begin{aligned} a_3-a_2 \sin (a_1+a)&= z^R \quad \Rightarrow \quad a_2=\frac{z^C-z^R}{\sin (a_1+a)} \nonumber \\ a_3 - a_2 \sin (a_1+a)\cosh x&=d \quad \,\,\,\Rightarrow \quad z^C-(z^C-z^R) \cosh x = d \nonumber \\ -a_2 \cos (a_1+a) \sinh x&=r \!\!\qquad \Rightarrow \quad -(z^C-z^R)\cot (a_1+a) \sinh x = r, \end{aligned}$$

(21)

so that

$$\begin{aligned} \tan (a_1+a) = \frac{z^R-z^C}{r}\sqrt{\cosh ^2 x-1}=\frac{z^R-z^C}{r}\sqrt{\left( \frac{d-z^C}{z^R-z^C} \right) ^2 -1} \end{aligned}$$

and

$$\begin{aligned} a_1= \arctan \left( \frac{1}{r} \sqrt{(d-z^C)^2-(z^R-z^C)^2} \right) -a. \end{aligned}$$

(22)

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 $-\infty$ to $\infty$. This is resolved by simply taking the conjugate of (19) as parametrization.

3 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+\mathrm {i}r$ (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 = \phi _k +{\mathrm {i}}\psi _k$, $k=1,\ldots ,M,$ on this curve, ordered in such a way that ${\text {Re}}z_k > {\text {Re}}z_{k+1}$. The parametric form of the inner parabola is

$$\begin{aligned} z(x) = -x^2 + z^R + \frac{\mathrm {i}r x}{\sqrt{z^R-d}}, \qquad x \in {\mathbb {R}}. \end{aligned}$$

(23)

Setting as in the previous section $z(x)=\phi +{\mathrm {i}}\psi$, with $\psi >0$, this means that if we fix the abscissa $\phi ={\text {Re}}{(z)}$ we obtain

$$\begin{aligned} \phi = z^R - x^2 \qquad \Longrightarrow \qquad x = \sqrt{z^R -\phi } \end{aligned}$$

which uniquely defines the argument of the parametrization $x > 0$, and consequently

$$\begin{aligned} \psi = \frac{r x}{\sqrt{z^R-d}}, \end{aligned}$$

which depends on r and d. We easily obtain

$$\begin{aligned} \frac{\partial \psi }{\partial d}= & {} \frac{x r}{2 (z^R-d)^{3/2}}\\ \frac{\partial \psi }{\partial r}= & {} \frac{x }{\sqrt{z^R-d}}. \end{aligned}$$

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 + \mathrm {i}r$, 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. [8, 17]).

Lemma 1

Let D(t) be a differentiable matrix-valued function in a neighborhood of $t_0$. Let

$$\begin{aligned} D(t) = U(t) \Sigma (t) V(t)^{*} = \sum \limits _{i} u_i(t) \sigma _i(t) v_i(t)^* \end{aligned}$$

(24)

be a smooth (with respect to t) singular value decomposition of the matrix D(t) and $\sigma (t)$ be a certain singular value of D(t) converging to a simple singular value ${{\hat{\sigma }}} \ne 0$ of $D_0 = D(t_0)$.

If ${{\hat{u}}}, {{\hat{v}}}$ are the associated left and right singular vectors, respectively, the function $\sigma (t)$ is differentiable near $t=t_0$ with

$$\begin{aligned} {\dot{\sigma }}(t_0) = {\text {Re}}\bigl ( {{\hat{u}}}^{*} {\dot{D}}_0 {{\hat{v}}} \bigr ) \qquad \text{ with } \ {\dot{D}}_0 = {\dot{D}}(t_0) . \end{aligned}$$

(25)

Proof

We have that $\sigma (t)^2$ is an eigenvalue of $D(t) D(t)^* = U(t) \Sigma (t)^2 U(t)^*$.

At $t=t_0$ the left and right eigenvectors associated to ${{\hat{\sigma }}}^2$ coincide and are equal to ${{\hat{u}}}$, having unit norm. Note that ${{\hat{u}}}$ is a certain column of $U(t_0)$ determined by the position of ${{\hat{\sigma }}}^2$ in the diagonal matrix $\Sigma (t_0)^2$.

Then—omitting the ubiquitous dependence on t - by [8, Theorem 6.3.12] we get for $t=t_0$

$$\begin{aligned} \frac{d}{dt} \sigma ^2(t) \Big |_{t=t_0} = 2 {{\hat{\sigma }}} {\dot{\sigma }}(t_0) = \frac{{{\hat{u}}}^* \bigl ( {\dot{D}}_0 D^*_0 + D_0 \dot{D}^*_0 \bigr ) {{\hat{u}}}}{{{\hat{u}}}^* {{\hat{u}}}} = 2 {\text {Re}}\bigl ( {{\hat{u}}}^* {\dot{D}}_0 D^*_0 {{\hat{u}}} \bigr ). \end{aligned}$$

Now using the fact (see (24)) that $D^*_0 {{\hat{u}}} = {{\hat{v}}} {{\hat{\sigma }}}$ we get (25). $\square$

3.1 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 rely 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 $\varepsilon$-pseudospectrum as the set

$$\begin{aligned} \sigma _{\varepsilon ,t}(A)=\left\{ z\in {\mathbb {C}} : \mathrm {e}^{{\text {Re}}(z)t}\Big \Vert (zI-A)^{-1}\Big \Vert \ge \frac{1}{\varepsilon }\right\} , \end{aligned}$$

(26)

which is equivalent to

$$\begin{aligned} \sigma _{\varepsilon ,t}(A)=\left\{ z\in {\mathbb {C}} : \mathrm {e}^{-{\text {Re}}(z)t}\sigma _{\min }\Big (zI-A\Big )\le \varepsilon \right\} . \end{aligned}$$

(27)

with $\sigma _{\min }\Big (zI-A\Big )$ denoting the smallest singular value of $zI-A$. In particular, we are interested on the boundary of the weighted $\varepsilon$-pseudospectral level set, this is

$$\begin{aligned} \partial \sigma _{\varepsilon ,t}(A)=\left\{ z\in {\mathbb {C}} : \mathrm {e}^{{\text {Re}}(z)t}\Big \Vert (zI-A)^{-1}\Big \Vert =\frac{1}{\varepsilon }\right\} . \end{aligned}$$

(28)

Note that for $t=0$ we recover the standard definition of pseudospectrum and $\varepsilon$-pseudospectrum.

For a fixed target $\varepsilon$, we impose the set in (28) to lay internal to the parabola. A natural approach to achieve this goal is defining the functional

$$\begin{aligned} {\mathcal {F}}(d,r) = \frac{1}{2} \sum \limits _{k=1}^{M} \Bigl (\varepsilon - {\tilde{\sigma }}_k(d,r) \Bigr )_{+}^2 \end{aligned}$$

where

$$\begin{aligned} {\tilde{\sigma }}_k(d,r)=\mathrm {e}^{-{\text {Re}}(z_k)t}\sigma _{\min }\Big (A - z_k \mathrm {I}\Big ), \end{aligned}$$

(29)

with $z_k$ depending on d and r, and $(\tau )_+ = \max \{ \tau ,0 \}$. In this way, values of ${\tilde{\sigma }}_k(d,r)$ which are larger than $\varepsilon$ do not contribute to the functional. The goal is to compute a solution (d, r) to

$$\begin{aligned} \min \limits _{d,r}&{\mathcal {F}}(d,r)\\ {{\mathrm{s. t.}}}&\min _{k=1, \ldots , M} {\tilde{\sigma }}_k(d,r) = \varepsilon \end{aligned}$$

where - for numerical convenience - the constraint may be treated as a penalization term. This can be done by computing the gradient of ${\mathcal {F}}$, where we use Lemma 1. The gradient is continuous and has the form

$$\begin{aligned} G(d,r) = \sum _{k=1}^{M} (\varepsilon -{\tilde{\sigma }}_k(d,r))_{+}\, {\text {Re}}(\mathrm {i}u_k^* v_k) \left( \begin{array}{cc} \displaystyle {\frac{x_k r}{2 (z^R-d)^{3/2}}} \\ \displaystyle {\frac{x_k}{\sqrt{z^R-d}}} \end{array} \right) , \end{aligned}$$

(30)

with $x_k$ such that $z(x_k)=z_k$. Note that ${\text {Re}}(z_k)$ does not depend on d and r, see (23). We add to the functional the penalization term

$$\begin{aligned} P(d,r) = \frac{1}{2} \left( \min _{k=1,\ldots ,M} {\tilde{\sigma }}_k(d,r) - \varepsilon \right) ^2, \end{aligned}$$

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

$$\begin{aligned} {\mathcal {F}}(d,r) + c P(d,r), \end{aligned}$$

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.

3.2 Selecting Points Internal to the Weighted $\varepsilon$-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 = \phi _1 + \mathrm {i}\psi _1$ and compute

$$\begin{aligned} \left( A - z_1 \mathrm {I}\right) = U_1 \Sigma _1 V_1, \end{aligned}$$

setting

$$\begin{aligned} {\tilde{\sigma }}_1 = \mathrm {e}^{-{\text {Re}}{(z_1)}t}\min {{\mathrm{diag}}} \left( \Sigma _1 \right) , \end{aligned}$$

according to the definition in (29). Then we check the difference $\delta \varepsilon = \varepsilon - {\tilde{\sigma }}_1$. If $\delta \varepsilon \le 0$ we proceed by considering $z_2$ and repeating the same steps, otherwise it means that we are inside the weighted $\varepsilon$-pseudospectrum and therefore we need to update the internal parabola. If we do not find any $z_k$ such that ${\tilde{\sigma }}_k < \varepsilon$, we may consider the point $z \in \{ z_k \}_{k=1}^{M}$ for which the corresponding ${\tilde{\sigma }}$ is minimal and proceed in order to find a closer parabola to the weighted $\varepsilon$-pseudospectral level set.

The algorithm we adopt tunes the interpolation point $w=d+\mathrm {i}r$, so that the updated curve is external weighted $\varepsilon$-pseudospectrum of A at z. We indicate by $p = p(d,r)$ a selected point $z_k=z(x_k)=\phi _k+{\mathrm {i}}\psi _k$ laying in the wrong weighted pseudospectral level set and write ${\tilde{\sigma }}_k$ as

$$\begin{aligned} {\tilde{\sigma }}(d,r)=\mathrm {e}^{-{\text {Re}}(p(d,r))t}\sigma _{\min }\Big (A-p(d,r)\mathrm {I}\Big ). \end{aligned}$$

(31)

In principle we want to solve the equation

$$\begin{aligned} {\tilde{\sigma }}(d,r) - \varepsilon = 0 \qquad \text{ w.r.t. } \ r. \end{aligned}$$

(32)

It seems natural to fix the parameter d as the mean of the abscissas of the support points $z_k$, i.e.

$$\begin{aligned} d = \frac{1}{M} \sum \limits _{k=1}^{M} \phi _k \end{aligned}$$

and solve the scalar equation

$$\begin{aligned} {\tilde{\sigma }}(d, r) - \varepsilon = 0 \qquad \text{ w.r.t. } \ r. \end{aligned}$$

(33)

Applying Lemma 1 to ${\tilde{\sigma }}(d,r)$- with u and v left and right associated singular vectors - we get

$$\begin{aligned} \frac{d}{d r} {\tilde{\sigma }} \left( d,r\right) = -\mathrm {e}^{-{\text {Re}}(p(d,r))t}{\text {Re}}(\mathrm {i}u^* v)\,g \end{aligned}$$

with

$$\begin{aligned} g = \frac{x_k }{\sqrt{z^R-d}}. \end{aligned}$$

In order to accurately compute r such that ${\tilde{\sigma }}(d,r)= \varepsilon$ we make a few (say m) Newton iterations

$$\begin{aligned} r^{\ell + 1} = r^{\ell } + \frac{\mathrm {e}^{-{\text {Re}}(p(d,r^\ell ))t} \sigma _{\min }\left( A-p(d,r^\ell )\mathrm {I}\right) -\varepsilon }{\mathrm {e}^{-{\text {Re}}(p(d,r^\ell ))t}{\text {Re}}\Bigl (\mathrm {i}(u^{\ell })^* v^{\ell }\Bigr )\,g}, \quad \ell =1,\ldots ,m-1 \end{aligned}$$

(34)

with $u^{\ell }$ and $v^{\ell }$ singular vectors associated to $\sigma _{\min }(A- p(d, r^{\ell }) \mathrm {I})$ and $r^{\ell }$ the actual ordinate of the interpolation point w.

Then we compute a new parabola, which interpolates $d + \mathrm {i}r^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 ${\tilde{\sigma }}(d,r)$ smaller than the target value $\varepsilon$. This would provide a simple variant to the method described above. We would first determine two points $z_1$ and $z_2$ such that ${\tilde{\sigma }}_j < \varepsilon$, for $j=1,2$ and then solve equations (in analogy to (33))

$$\begin{aligned}&{\tilde{\sigma }}_1(d,r) - \varepsilon = 0 \\&{\tilde{\sigma }}_2(d,r) - \varepsilon = 0, \end{aligned}$$

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.

4 A Case Study: The 1D Black and Scholes Equation

The well known (deterministic) Black-Scholes equation [3] has the following form:

$$\begin{aligned} \frac{\partial u}{\partial \tau }=\frac{1}{2}\sigma ^2s^2\frac{\partial ^2u}{\partial s^2}+rs\frac{\partial u}{\partial s}-ru,\quad \;\;s>L,\;\;0<\tau \le t, \end{aligned}$$

(35)

for L, t given, where the unknown function $u(s,\tau )$ stands for the fair price of the option when the corresponding asset price at time $t-\tau$ is s and t is the maturity time of the option. Moreover, $r\ge 0$, $\sigma >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

$$\begin{aligned} L<s<S \end{aligned}$$

for a sufficiently large S. We take (35) together with the following conditions, typical for the European call option, cf. [12]:

$$\begin{aligned} \begin{aligned} u(s,0)&=\max (0,s-K),\\ u(L,\tau )&=0,\;0\le \tau \le t,\\ u(S,\tau )&=S-\mathrm {e}^{-r\tau }K,\;\;0\le \tau \le t, \end{aligned} \end{aligned}$$

(36)

being K the reference strike price. In this Section we extend the theory developed in [21] for the 1D convection-diffusion operator

$$\begin{aligned} {\mathcal {L}}u = u_{xx}+u_x. \end{aligned}$$

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 Sect. 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

$$\begin{aligned} v(s,\tau )=u(s,\tau ) - y(s,\tau ), \end{aligned}$$

with

$$\begin{aligned} y(s,\tau )= \frac{s-L}{S-L}\left( S - \mathrm {e}^{-r\tau }K \right) . \end{aligned}$$

The differential equation for v reads

$$\begin{aligned}&\frac{\partial v}{\partial \tau }{=}\frac{1}{2}\sigma ^2s^2\frac{\partial ^2v}{\partial s^2}{+}rs\frac{\partial v}{\partial s}{-}rv {-} {\frac{s}{S{-}L}}r\mathrm {e}^{{-}r\tau } K{+}{\frac{Lr}{S{-}L}S}{,} s>L,\;\;0<\tau \le t, \end{aligned}$$

(37)

with initial and boundary data

$$\begin{aligned} \begin{aligned} v(s,0)&= \max (0,s-K)-\frac{s-L}{S-L}\left( S - K \right) =: v_0(s),\\ v(L,\tau )&=0,\;0\le \tau \le t,\\ v(S,\tau )&=0,\;\;0\le \tau \le t, \end{aligned} \end{aligned}$$

(38)

We can now apply the Laplace transform to both sides of (37) with respect to $\tau$. This leads to the following equation for V(s, z), the Laplace transform of $v(s,\tau )$:

$$\begin{aligned} V(s,z) = (zI-{\mathcal {L}})^{-1} \left( v_0(s) - \frac{srK}{S-L}\frac{1}{z+r}+\frac{Lr}{S-L}\frac{1}{z} \right) , \end{aligned}$$

with ${\mathcal {L}}$ the differential operator for the Black-Scholes problem with homogeneous boundary conditions.

4.1 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={\mathcal {L}}u$, with

$$\begin{aligned} {\mathcal {L}}u=\frac{1}{2}\sigma ^2 u_{xx}+ \left( r-\frac{1}{2}\sigma ^2 \right) u_x-ru,\qquad 0\le x\le \log (S), \end{aligned}$$

(39)

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 ${\mathcal {L}}$ by applying it to a mapping of the form $\varphi (x)=\mathrm {e}^{\alpha x}$. In this way, we obtain

$$\begin{aligned} {\mathcal {L}}\varphi =(\nu \alpha ^2+(r-\nu )\alpha -r)\mathrm {e}^{\alpha x} =\lambda \varphi ,\qquad \text {with}\;\;\nu =\frac{1}{2}\sigma ^2 \end{aligned}$$

(40)

and

$$\begin{aligned} \lambda = \nu \alpha ^2+(r-\nu )\alpha -r. \end{aligned}$$

(41)

Then, for each $\lambda$ real we have two associated values of $\alpha$, namely

$$\begin{aligned} \alpha _{\pm }=\frac{-(r-\nu )\pm \sqrt{(r+\nu )^2+4\lambda \nu }}{2\nu }. \end{aligned}$$

(42)

For any $\lambda$ and corresponding $\alpha _+$ and $\alpha _-$, the function

$$\begin{aligned} \phi (x)=\frac{\mathrm {e}^{\alpha _+x}-\mathrm {e}^{\alpha _-x}}{\alpha _+-\alpha _-} \end{aligned}$$

(43)

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 $\mathrm {e}^{\alpha _+ \log (S)}=\mathrm {e}^{\alpha _- \log (S)}$, that is, $(\alpha _+-\alpha _-)\log (S)=2\pi n {\mathrm {i}}$ for some nonzero $n\in {\mathbb {Z}}$. By (42), this amounts to the condition $(\log (S)/\nu )\sqrt{(r+\nu )^2+4\lambda \nu }=2\pi n {\mathrm {i}}$, and upon squaring we obtain the following eigenvalues

$$\begin{aligned} \lambda _n=-\left( \frac{r+\nu }{2}\right) ^2\frac{1}{\nu } -\frac{\pi ^2n^2\nu }{\log (S)^2},\;\;\;\;n=1,2,3,... \end{aligned}$$

(44)

Thus $\Lambda ({\mathcal {L}})$ is a discrete set of negative real numbers in the interval $(-\infty ,-\frac{1}{4})$. Note that, for our problem, there are choices of $\lambda$ for which both $\alpha _+$ and $\alpha _-$ lie in the left half-plane and thus both $\mathrm {e}^{\alpha _+x}$ and $\mathrm {e}^{\alpha _-x}$ are decreasing functions. For the eigenfunctions associated to (44), this occurs with ${\text {Re}}(\alpha _+)={\text {Re}}(\alpha _-)=-\frac{(r-\nu )}{2\nu }$ under the assumption $r>\nu$. More generally, it occurs if and only if $\alpha$ belongs to the strip $B=\{\alpha \in {\mathbb {C}}:\;-\frac{(r-\nu )}{\nu }\le {\text {Re}}(\alpha )\le 0\}$, since if $\alpha$ is one solution of (40), the other is $-\frac{(r-\nu )}{\nu }-\alpha$. The corresponding region in the $\lambda$-plane is the image of B under the function $\lambda =\nu \alpha ^2+(r-\nu )\alpha -r$, which we denote by $\Pi$:

$$\begin{aligned} \Pi =\left\{ \lambda \in {\mathbb {C}}:\;\lambda =\nu \alpha ^2+(r-\nu )\alpha -r,\; -\frac{(r-\nu )}{\nu }\le \text {Re}(\alpha )\le 0\right\} . \end{aligned}$$

(45)

The “critical parabola” that bounds $\Pi$ is the image of the boundary of B under the same function, which we can simply represent by

$$\begin{aligned} P=\{\lambda \in {\mathbb {C}}:\;\lambda =\nu \alpha ^2 +(r-\nu )\alpha -r,\;{\text {Re}}(\alpha )=0\} \end{aligned}$$

(46)

since $\text {Re}(\alpha )=-\frac{(r-\nu )}{\nu }$ maps onto the same parabola as $\text {Re}(\alpha )=0$. Suppose now that $\lambda$ is any complex number in the interior of $\Pi$ so that ${\text {Re}}(\alpha _+)<0$ and ${\text {Re}}(\alpha _-)<0$. Then $\phi (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 $\mathrm {e}^{\mu \log (S)}=S^{\mu }$, where $\mu =\max \{{\text {Re}}(\alpha _+),{\text {Re}}(\alpha _-)\}$. Thus $\phi (x)$ is nearly an eigenfunction of ${\mathcal {L}}$, though $\lambda$ 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.

Theorem 3

Let $\lambda$ be an arbitrary point in the interior of $\Pi$, with $\alpha _{\pm }$ and $\phi (x)$ defined by (42) and (43), and write $\alpha _+-\alpha _-=(1/\nu )\sqrt{(r+\nu )^2+4\lambda \nu }=\beta +{\mathrm {i}}\tau$. Then

$$\begin{aligned} \Vert (\lambda I-{\mathcal {L}})^{-1}\Vert \sim \frac{\Vert \phi \Vert ^2}{\phi (\log (S))} \end{aligned}$$

(47)

If in addition $\lambda \not \in (-\infty ,\;-\left( \frac{r+\nu }{2}\right) ^2\frac{1}{\nu }]$, then $\phi (\log (S))\sim S^{\mu }/|\alpha _+-\alpha _-|$ and therefore

$$\begin{aligned} \Vert (\lambda I-{\mathcal {L}})^{-1}\Vert \sim {S^{-\mu }(\beta ^2+\tau ^2)^{1/2}\Vert \phi \Vert ^2}, \end{aligned}$$

(48)

where $\mu =\max \{{\text {Re}}(\alpha _+),{\text {Re}}(\alpha _-)\}<0$.

This result tells us that the resolvent norm changes exponentially along any vertical line inside $\Pi$. Indeed we know, by construction, that the critical parabola is the curve such that ${\text {Re}}(\alpha )=0$ and therefore $S^{-\mu }=1$, while on the real axis, for ${\text {Re}}(\lambda )$ sufficiently small the real part of $\alpha _-$ and $\alpha _+$ are the same and equal to $-\frac{r-\nu }{2\nu }$, that corresponds to the case where $S^{-\mu }$ is maximized.

4.2 Symmetrizability and a Further Estimate

As done in [21] we can explicitly symmetrize the BS operator ${\mathcal {L}}$. First let’s define $\rho =\frac{r-\nu }{2\nu }$ and $u(x)=\mathrm {e}^{-\rho x}v(x)$, which implies

$$\begin{aligned} u'&=\mathrm {e}^{-\rho x}\left( -\rho v+v'\right) , \\ u''&=\mathrm {e}^{-\rho x}\left( \rho ^2v-2\rho v'+v''\right) , \end{aligned}$$

and therefore

$$\begin{aligned} {\mathcal {L}}u&=\nu u''+(r-\nu )u'-ru\\&=\mathrm {e}^{-\rho x}\left[ (\rho ^2v-2\rho v'+v'')\nu +(-\rho v+v')(r-\nu )-rv\right] \\&=\mathrm {e}^{-\rho x}\left[ \rho ^2\nu v-(r-\nu )v'+\nu v''-\rho (r-\nu )v+(r-\nu )v'-rv\right] \\&=\mathrm {e}^{-\rho x}\left[ v''\nu +v(\rho ^2\nu -\rho (r-\nu )-r)\right] \\&=\mathrm {e}^{-\rho x}\left[ \nu v''-\left( \frac{(r-\nu )^2}{4\nu }-r\right) v\right] . \end{aligned}$$

Thus if we define ${\mathcal {K}}v=\nu v''-\left( \frac{(r-\nu )^2}{4\nu }-r\right) v$, ${\mathcal {M}}v=\mathrm {e}^{-\rho x}v(x)$, then we have

$$\begin{aligned} {\mathcal {L}}=\mathcal {MKM}^{-1}. \end{aligned}$$

(49)

Here ${\mathcal {K}}$ is a self-adjoint operator and ${\mathcal {M}}$ is a diagonal operator with $\Vert {\mathcal {M}}\Vert =1$, $\Vert {\mathcal {M}}^{-1}\Vert =\mathrm {e}^{\rho \log (S)}$, and consequently

$$\begin{aligned} \kappa ({\mathcal {M}})=\Vert {\mathcal {M}}\Vert \Vert {\mathcal {M}}^{-1}\Vert =S^{\frac{r-\nu }{2\nu }}. \end{aligned}$$

(50)

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>\nu$ and $\lambda \in {\mathbb {C}}$,

$$\begin{aligned} \Vert (\lambda I-{\mathcal {L}})^{-1}\Vert \le \frac{S^{\frac{r-\nu }{2\nu }}}{\mathrm {dist}(\lambda ,P)}\le \frac{S^{\frac{r-\nu }{2\nu }}}{|{\text {Im}}( \lambda )|}. \end{aligned}$$

(51)

The discussion done in the previous subsection this theorem suggested that the pseudospectra level curves of ${\mathcal {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 $|\lambda |\rightarrow 0$ but, for any fixed $\varepsilon$ and S, $\Lambda _{\varepsilon }({\mathcal {L}})$ the $\varepsilon$-pseudospectra is contained in a strip of finite (though typically large) width:

$$\begin{aligned} \Lambda _{\varepsilon }({\mathcal {L}})\subset \{\lambda \in {\mathbb {C}}: \;|{\text {Im}}(\lambda )|\le \varepsilon \; S^{\frac{r-\nu }{2\nu }}\}. \end{aligned}$$

(52)

4.3 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 ${\mathcal {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 Fig. 3 where we set $S=200$ and used 2000 points for the space discretization; the resolvent norm is plotted in a subset of the $\Pi$ 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.

Fig. 3

Magnitude of the estimate of the resolvent norm (48) (left) and magnitude of the computed resolvent norm (right)

5 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 Sect. 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 $\mathrm {e}^{z^L t} = \varepsilon$, with $\varepsilon$ the working precision.

• 2. Compute the critical curve $\Gamma _{left}$ according to Sect. 3. This procedure provides an explicit parametrization of $\Gamma _{left}$, of the form $\psi (x)$, $x\in {\mathbb {R}}$.

• 3. Compute $c_{\max }$ as the unique value that satisfies ${{\text {Re}}(\Gamma (c_{\max } \pi ) )= z^L}$. This gives, according to the discussion in Sect. 2:

  $$\begin{aligned} c_{\max }(a) = \left\{ \begin{array}{ll} \displaystyle \frac{1}{2}, &{}\quad \text{ for } \text{ the } \text{ ellipse },\\ \displaystyle {\frac{1}{\pi }\sqrt{z^R-z^L+a^2+\frac{ar}{\sqrt{z^R-d}}}}, &{}\quad \text{ for } \text{ the } \text{ parabola }, \\ \displaystyle {\frac{1}{\pi } \log \left( b+\sqrt{b^2-1 } \right) },\;\text {with}\;\\ b=\frac{(z^C-z^L)\sin (a_1+a)}{(z^R-z^L)\sin (a_1)},&{} \quad \text{ for } \text{ the } \text{ hyperbola }. \end{array} \right. \end{aligned}$$

  (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

  $$\begin{aligned} D(a) = \left\{ \begin{array}{ll} \displaystyle z(-{\mathrm {i}}a) = z^L+ \cosh (2a)(z^R-z^L)\\ +\sinh (2a) \frac{r}{\sin (\theta )}, &{}\quad \text{ for } \text{ the } \text{ ellipse },\\ \displaystyle z(-{\mathrm {i}}a) = a^2-2 a_1(a)a +a_2(a), &{}\quad \text{ for } \text{ the } \text{ parabola }, \\ \displaystyle z({\mathrm {i}}a) = z^C-a_2(a)\sin (a_1(a) -a),&{}\quad \text{ for } \text{ the } \text{ hyperbola }, \end{array} \right. \end{aligned}$$

  (54)

  where r in the first expression above is the imaginary part of the control point $d+{\mathrm {i}}r$ from Sect. 2 and $\theta =\arccos \left( \frac{d-z^L}{z^R-z^L}\right)$, as in (13) and (14). In order to have a more robust estimate we take into account the constant $M_{left}$, named $M_{+}$ in [6, equation (23)], which is defined as

  $$\begin{aligned} M_{left}=\frac{1}{2\pi }\max _{z\in {\tilde{\Gamma }}_{left}} \Big \Vert \mathrm {e}^{zt}(zI-A)^{-1}\left( u_0+{\hat{b}}(z)\right) z'\Big \Vert , \end{aligned}$$

  (55)

  where ${\tilde{\Gamma }}_{left}$ is a suitable restriction of $\Gamma _{left}$.

  We also consider a different estimate for $M_{right}$, named $M_{-}$ in [6, equation (24)], which is defined as

  $$\begin{aligned} M_{right}=\frac{1}{2\pi }\max _{z\in {\tilde{\Gamma }}_{right}} \Big \Vert \mathrm {e}^{zt}(zI-A)^{-1}\left( u_0+{\hat{b}}(z)\right) z'\Big \Vert , \end{aligned}$$

  (56)

  where ${\tilde{\Gamma }}_{right}$ is a suitable restriction of $\Gamma _{right}$. Note that $M_{left}$ can be bounded as

  $$\begin{aligned} M_{left}\le \frac{1}{2\pi }\max _{z\in {\tilde{\Gamma }}_{left}} \mathrm {e}^{{\text {Re}}(z)t}\Big \Vert (zI-A)^{-1}\Big \Vert \Big \Vert \left( u_0+{\hat{b}}(z) \right) \Big \Vert |z'|={\tilde{M}}_{left}, \end{aligned}$$

  (57)

  and we can take advantage of the computation done in step 2 for the weighted $\varepsilon$-pseudospectrum to approximate ${\tilde{M}}_{left}$. Concerning $M_{right}$ we first assume that the maximum is reached on $z_v(a)$, the vertex of ${\tilde{\Gamma }}_{right}$; in this way we have

  $$\begin{aligned} M_{right}\le \frac{1}{2\pi }\mathrm {e}^{D(a)t}\Big \Vert (z_v(a)I-A)^{-1}\Big \Vert \left\| u_0+{\hat{b}}\left( z_v(a)\right) \right\| |z'_v(a)|={\tilde{M}}_{right}. \end{aligned}$$

  (58)

  Therefore repeating the calculations of [6, Section 3.2] and including $M_{left}$ and $M_{right}$ we have, for a fixed target accuracy $tol$,

  $$\begin{aligned} N=&\frac{c}{a}\Big (\log \left( 2\pi c{\tilde{M}}_{right}+\pi {\tilde{M}}_{left} \right) -\log \left( tol\right) \Big )\nonumber \\ \le&\frac{c_{\max }(a)}{a}\Big (\log \left( 2\pi c_{\max }(a){\tilde{M}}_{right} +\pi {\tilde{M}}_{left} \right) -\log \left( tol\right) \Big ). \end{aligned}$$

  (59)

  

  The value of a is then computed by minimizing the right hand side of (59) which requires an interval of the form $a\in [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 Sect. 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

    $$\begin{aligned} f(a)=\frac{c_{\max }(a)}{a}\Big (\log \left( 2\pi c_{\max }(a){\tilde{M}}_{right}+\pi {\tilde{M}}_{left} \right) -\log \left( tol\right) \Big )+wp(a), \end{aligned}$$

    (60)

    where w is a positive scalar and

    $$\begin{aligned} p(a)={\left\{ \begin{array}{ll} 0 &{}\quad err^{num}_N(a) < tol\\ 1 &{}\quad err^{num}_N(a) \ge tol\end{array}\right. } \end{aligned}$$

    (61)

  • (ii) We note that evaluating (58) for every $a\in [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.

• 5. Compute a truncation parameter $c\le c_{\max }$. The actual numerical integration is performed for $x\in [-c\pi , c\pi ]$. 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.

• 6. Set

  $$\begin{aligned} N= \left\lceil \frac{c}{a}\Big (\log \left( 2\pi c{\tilde{M}}_{right}+\pi {\tilde{M}}_{left}\right) -\log \left( tol\right) \Big )\right\rceil . \end{aligned}$$

  (62)

6 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, Algorithm 3.2].

We also notice that in all our tests we construct the inner curves as explained in Sect. 3, taking the weighted $\varepsilon$-pseudospectral level set with $\varepsilon =10^{-7}$.

6.1 Black-Scholes Equation

Following the same strategy adopted in [6, 13], 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 ${\mathcal {O}}(\Delta x^2)$, with $\Delta x$ the diameter of the spatial grid, being around 0.001 for our choice of $N_h$. This can be observed in Fig. 4, where we consider a reference solution computed with $N_h=2\cdot 10^4$ spatial nodes and approximated in time by using expmv.

We set $r=0.06$, $\sigma =0.05$, and $K=80$. We plot the error for a selection of tolerances for the cases $t=1$ (Fig. 5) and $t=10$ (Fig. 6). In Fig. 7 we show the selected profiles of integration for the tolerance $tol=5\cdot 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\cdot 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)).

Fig. 4

Error vs number of nodes for Black-Scholes, $t=1$. The asymptotic value reached corresponds to the error in space when $N_h=2000$

Fig. 5

Error vs number of nodes for Black-Scholes, $t=1$. Comparison for different values of the tolerance among the elliptic, parabolic and hyperbolic contours

Fig. 6

Error vs number of nodes for Black-Scholes, $t=10$. Comparison for different values of the tolerance among the elliptic, parabolic and hyperbolic contours

Fig. 7

Example of integration profiles for the Black-Scholes problem for tolerance $tol=5\cdot 10^{-5}$ at time $t=1$ (left) and $t=10$ (right)

6.2 Heston Equation

The Heston equation [7] is given by

$$\begin{aligned} \frac{\partial u}{\partial \tau }=\frac{1}{2}s^2v\frac{\partial ^2u}{\partial s^2}+\rho \sigma sv\frac{\partial ^2 u}{\partial s \partial v} +\frac{1}{2}\sigma ^2v\frac{\partial ^2u}{\partial v^2}+(r_d-r_f)s \frac{\partial u}{\partial s}+\kappa (\eta -v)\frac{\partial u}{\partial v}-r_du. \end{aligned}$$

(63)

The unknown function $u(s, v,\tau )$ represents the price of a European option when at time $t-\tau$ the corresponding asset price is equal to s and its variance is v. We consider the equation on the unbounded domain

$$\begin{aligned} 0\le \tau \le t,\;\;s>0,\;v>0, \end{aligned}$$

where the time t is fixed. The parameters $\kappa >0$, $\sigma >0$, and $\rho \in [-1,1]$ are given. Moreover Eq. (63) is usually considered under the condition $2\kappa \eta >\sigma ^2$ that is known as the Feller condition (see [15]). We take equation (63) together with the initial condition

$$\begin{aligned} u(s,v,0)=\max (0,s-K), \end{aligned}$$

where $K>0$ is fixed a priori (and represents the strike price of the option), and boundary condition

$$\begin{aligned} u(L,v,\tau )=0,\;\;0\le \tau \le t. \end{aligned}$$

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

$$\begin{aligned} 0\le s\le S,\;\; 0\le v\le V. \end{aligned}$$

On the new boundary, we need to add two more conditions (specific for the European call option),

$$\begin{aligned} \frac{\partial u}{\partial s}(S,v,\tau )&=\mathrm {e}^{-r_f\tau },\;\;0\le \tau \le t,\\ u(s,V,\tau )&=s\mathrm {e}^{-r_f\tau },\;\;0\le \tau \le t, \end{aligned}$$

which are treated analogously to the boundary condition in (36). We use the spatial discretization proposed in [10] for $k=1.5$, $\eta =0.04$, $\sigma =0.3$, $\rho =-0.9$, $r_d=0.025$, $r_f=0$, $K=100$, $L=0$, $S=8K$, $V=5$. In Fig. 8, we plot the error for a selection of tolerances for $t=1$ and in Fig. 9 we do the same for $t=10$. In Fig. 10 we show the selected profiles of integration for $tol=5\cdot 10^{-5}$ at $t=1$ and $t=10$.

Table 1 Estimate (62) for the Black-Scholes test problem at time $t=1$ and $t=10$. Results for the elliptic, parabolic and hyperbolic profiles

Fig. 8

Error vs number of nodes for Heston, $t=1$. Comparison for different values of the tolerance among the elliptic, parabolic and hyperbolic contours

Fig. 9

Error vs number of nodes for Heston, $t=10$. Comparison for different values of the tolerance among the elliptic, parabolic and hyperbolic contours

Fig. 10

Example of integration profiles for the Heston problem for tolerance $tol=5\cdot 10^{-5}$ at time $t=1$ (left) and $t=10$ (right)

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.

6.3 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

$$\begin{aligned} t_1=\Lambda t_0,\;\Lambda >1. \end{aligned}$$

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 $\Gamma _{left}$ as explained in Sect. 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 Sect. 5. An application of the construction described in Sect. 3 gives the interpolation point $d+ir$, which uniquely defines $\Gamma _{left}$. Then we determine a by following the procedure described in Sect. 5 with $t=t_1$, the upper end of the time interval. We use the profile determined in this way for all $t\in [t_0,t_1]$.

Table 2 Estimate (62) for the Heston test problem at time $t=1$ and $t=10$. Results for the elliptic, parabolic and hyperbolic profiles

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 combination

$$\begin{aligned} {\tilde{I}}_N=\frac{c}{Ni}\sum _{j=1}^{N-1}\mathrm {e}^{z(x_i)t}{\hat{u}}_jz'(x_j), \end{aligned}$$

(64)

where ${\hat{u}}_j={\hat{u}}(z(x_j))+\rho _j$ and $\rho _j$ is the error in the numerical solution of the linear system

$$\begin{aligned} ((z(x_j))I-A){\hat{u}}=u_0+{\hat{b}}(z(x_j)), \end{aligned}$$

(65)

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 by

$$\begin{aligned} {\tilde{err}}_N=\big |u(t)-{\tilde{I}}_N\big |, \end{aligned}$$

(66)

that we can estimate in the following way:

$$\begin{aligned} {\widetilde{err}}_N=\Big |u(t)-\frac{c}{Ni}\sum _{j=1}^{N-1}\mathrm {e}^{z(x_i)t} {\hat{u}}_jz'(x_j)\Big |\le err_T+err_N+err^{num}_N, \end{aligned}$$

(67)

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

$$\begin{aligned} {err}_N^{num}&=\Bigg |\frac{c}{Ni}\sum _{j=1}^{N-1}\mathrm {e}^{z(x_i)t}({\hat{u}} (z(x_j))-{\hat{u}}_j)z'(x_j)\Bigg |\nonumber \\&\le \frac{c_{\max }}{N}\sum _{j=1}^{N-1}\mathrm {e}^{{\text {Re}}{(z(x_i))}t}|\rho _j||z'(x_j)|. \end{aligned}$$

(68)

Therefore, once the integration contour for a time window $[t_0,\;t_1]$ is fixed, we can observe that $err_N^{num}$ does not remain constant but instead changes exponentially with respect to time. Since we want that ${\widetilde{err}}_N \le tol$ for all $t\in [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\in [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, \Lambda 0.1]$ and $[1, \Lambda 1]$ for $\Lambda =10$. In the plots of Figs. 11, 12, 13, 14 we show the numerical results for Black-Scholes and Heston equations. The target tolerance chosen is $tol=5\cdot 10^{-8}$ for Black-Scholes and $tol=5\cdot 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.

Fig. 11

Black-Scholes equation in time interval [0.1, 1], $tol=5\cdot 10^{-8}$, $z^L=-400$, $z^R=0.01$

Fig. 12

Black-Scholes equation in time interval [1, 10], $tol=5\cdot 10^{-8}$, $z^L=-40$, $z^R=0.01$

Fig. 13

Heston equation in time interval [0.1, 1], $tol=5\cdot 10^{-4}$, $z^L=-400$, $z^R=0.06$

Fig. 14

Heston equation in time interval [5.5, 10], $tol=5\cdot 10^{-4}$, $z^L=-40$, $z^R=0.06$

7 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 term ${\hat{b}}$ 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.

7.1 Construction of the Inner Curve and Integration Map

The first step is the computation of the inner curve $\Gamma _{left}$, as described in Sect. 3; then we determine the bandwidth of the map a and the truncation parameter c. Concerning $\Gamma _{left}$ 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^{\ell +1}\gg r^{\ell }$. If this is the case, we do the following update: $r^{\ell +1}=r^{\ell }+pr^{\ell }$ 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 $\delta \varepsilon$ in the coarser grid; if $\delta \varepsilon \le 0$ we consider the successive point in the same coarse grid. Otherwise we pass to the finer one and we consider points there until $\delta \varepsilon >0$. This is all implemented in the functions

• 1. Elliptic_Map,

• 2. Parabolic_Map,

• 3. Hyperbolic_Map.

The three functions have the same structure, it only changes the type of contour, therefore we limit our description to Parabolic_Map.

7.1.1 An Example of Implementation: The Function Parabolic_Map

The following are the input and output arguments:

• Input:

  • 1. A the discrete operator;

  • 2. ${\hat{b}}(z)$ Laplace transform of the known term;

  • 3. $u_0$ initial solution;

  • 4. $A_r$ operator of smaller size than A used to compute the internal curve $\Gamma _{left}$;

  • 5. $\varepsilon _1$ target value of the weighted pseudospectral level curve;

  • 6. T time where to evaluate the solution;

  • 7. $n_X$ maximum number of points where I compute the pseudospectra;

  • 8. $z^L$ minimum real value;

  • 9. $z^R$ the internal parabola vertex position;

  • 10. $tol$ accuracy required.

• Output:

  • 1. $a_p$ uniquely defines the map, see Sect. 2;

  • 2. c truncation value of the integral, see Sect. 5;

  • 3. $a_1$ and $a_2$ defined in (16) and (17) respectively;

  • 4. N number of nodes required to integrate with accuracy $tol$.

The structure of the function is reported in Algorithm 3.

7.2 The Integral Approximation

Once the map and the truncation value c are known, the integration related to the numerical inversion of the Laplace transform can be performed. This is implemented in function InLa_Quadrature whose structure is reported in Algorithm 4.

The new variables introduced are: u the output vector that contains the nodal values of the solution at time t and flag, an input variable that specify which profile of integration has to be used: $flag=1$ to the elliptic profile, $flag=2$ the parabolic one and finally $flag=3$ the hyperbolic.

7.3 Amplitude of the Time Window

In Sect. 6.3 we have established that the amplitude for a time window $[t_0, t_1]$ is acceptable if (68) is smaller than $tol$ when evaluated at $t_0$ and $t_1$. To compute (68) we first need to approximate the error in the numerical solution $\rho _j$, which is not simple since it depends on many factors: the condition number of matrix $A_j=z(x_j)I-A$, the perturbations in the matrix and in the r.h.s., the stability of the algorithm to solve the linear system and machine precision of the solver. Assuming the matrix and the r.h.s. are exact, we can approximate $\rho _j$ by means of the residual associated to the solution of the linear system, i.e.

$$\begin{aligned} A_j\rho _j= A_j({\hat{u}}(z(x_j))-{\hat{u}}_j)=b-A_j({\hat{u}}_j)=r_j, \end{aligned}$$

(69)

so that

$$\begin{aligned} \Vert \rho _j\Vert \le \Vert A_j^{-1}\Vert \Vert r_j\Vert . \end{aligned}$$

(70)

Expression (70) needs the evaluation of the numerical solution at the nodes $x_j$ and the computation of the smallest non zero singular value of $A_j$. Since $\kappa (A_j)$ (and $\Vert A_j^{-1}\Vert$) may change significantly when evaluated at different nodes, we can expect that (68) depends only weakly on N. This fact allows us to use only few nodes to evaluate (68), which results into a smaller computational effort with respect Algorithm 4. However, the extra computational effort is justified by the fact that the solution is made available on a continuous time window with an error below the prescribed accuracy for all $t\in [t_0, t_1]$.

8 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.