Following this procedure for $p_{k,s-j}, j = 2, \dots, s-1$, we arrive at the Runge-Kutta formulation (4.26) of the double adjoint of scheme (4.25), where we go back to the notation $p_{ki}$ instead of $\hat{p}_{ki}$. $\square$

**Remark 4.3.4.** A straightforward calculation yields that without damping (for $\eta = 0$), we have $\alpha_{s-j} = j + 1$, and method (4.25)-(4.26) reduces to the following recurrence

$$
\begin{align*}
y_{k0} &= y_k, \quad y_{k1} = y_{k0} + \frac{h}{s^2} f(u_{k0}, y_{k0}), \\
y_{ki} &= \frac{2h}{s^2} f(u_{k,i-1}, y_{k,i-1}) + 2y_{k,i-1} - y_{k,i-2}, \quad i = 2, \dots, s, \\
y_{k+1} &= y_{ks}, \\
p_N &= \nabla \Psi(y_N), \quad p_{ks} = p_{k+1}, \\
p_{k,s-1} &= p_{ks} + \frac{h}{s^2} \nabla_y H(u_{k,s-1}, y_{k,s-1}, p_{ks}), \\
p_{k,s-j} &= \frac{2j}{(j+1)s^2} h \nabla_y H(u_{k,s-j}, y_{k,s-j}, p_{k,s-j+1}) + \frac{2j}{j+1} p_{k,s-i+1} \\
&\quad + \frac{1-j}{j+1} p_{k,s-j+2}, \quad j = 2, \dots, s-1, \\
p_{k0} &= \frac{h}{s} \nabla_y H(u_{k0}, y_{k0}, p_{k1}) + sp_{k1} + (1-s)p_{k2}, \\
p_k &= p_{k0}, \\
\nabla_u H(u_{k,s-j}, y_{k,s-j}, p_{k,s-j+1}) &= 0, \quad j = 1, \dots, s.
\end{align*}
\tag{4.30}
$$

where $k = 0, \dots, N-1$. In Section 4.3.4, we shall study the stability of (4.30) (without damping) and of (4.25)-(4.26) (with damping).

### 4.3.3 RKC method of order 2

We consider the new implementation (4.21) of the RKC method applied to (4.1) for which the internal stages behave similarly to that of the order one method, given by

$$
\begin{align}
\min \Psi(y_N), \quad \text{such that} \nonumber \\
y_{k0} &= y_k, \quad y_{k1} = y_{k0} + \mu_1 h f(u_{k0}, y_{k0}), \nonumber \\
y_{ki} &= \mu_i h f(u_{k,i-1}, y_{k,i-1}) + \nu_i y_{k,i-1} + (1-\nu_i) y_{k,i-2}, \quad i = 2, \dots, s \tag{4.31} \\
y_{k+1} &= a_s y_{k0} + b_s T_s(\omega_0) y_{ks}, \nonumber
\end{align}
$$

where, $k = 0, \dots, N-1$, $\eta = 0.15$, and again all the parameters are defined as for the Chebyshev method using $\omega_2$ instead of $\omega_1$. The order two RKC method with formulation (4.31) for the state equation can be implemented using Algorithm 4.2.3, the costate equation being implemented using the recurrence relations given by the following Theorem 4.3.5. Its proof is analogous to that of Theorem 4.3.3 and thus omitted.
<figure><img src="image_42.png" /><figcaption>Figure 4.4: Internal stages (thin curves) and stability polynomial (bold curve) of the double adjoint of the Chebyshev method (4.26) of order one and the RKC method (4.32) of order two for $s = 10$ internal stages.</figcaption></figure>

**Theorem 4.3.5.** The double adjoint of the scheme (4.31) is given by the recurrence

$$ 
\begin{aligned}
p_N &= \nabla \Psi(y_N), \quad p_{ks} = p_{k+1}, \\
p_{k,s-1} &= p_{ks} + \frac{\mu_s}{\nu_s} h \nabla_y H(u_{k,s-1}, y_{k,s-1}, p_{ks}), \\
p_{k,s-j} &= \frac{\mu_{s-j+1}\alpha_{s-j+1}}{\alpha_{s-j}} h \nabla_y H(u_{k,s-j}, y_{k,s-j}, p_{k,s-j+1}) \\
&\quad + \frac{\nu_{s-j+1}\alpha_{s-j+1}}{\alpha_{s-j}} p_{k,s-j+1}, \\
&\quad + \frac{(1-\nu_{s-j+2})\alpha_{s-j+2}}{\alpha_{s-j}} p_{k,s-j+2}, \quad j = 2, \dots, s-1, \\
p_{k0} &= \mu_1 \alpha_1 h \nabla_y H(u_{k0}, y_{k0}, p_{k1}) + \alpha_1 p_{k1} + (1-\nu_2) \alpha_2 p_{k2} + a_s p_{k+1}, \\
p_k &= p_{k0}, \\
\nabla_u H(u_{k,s-j}, y_{k,s-j}, p_{k,s-j+1}) &= 0, \quad j = 1, \dots, s,
\end{aligned}
\quad (4.32) $$

where the coefficients $\alpha_j$ are defined using the induction

$$ 
\begin{aligned}
\alpha_s &= b_s T_s(\omega_0), \quad \alpha_{s-1} = \nu_s \alpha_s, \\
\alpha_{s-j} &= \nu_{s-j+1} \alpha_{s-j+1} + (1-\nu_{s-j+2}) \alpha_{s-j+2}, \quad j = 2 \dots s-1.
\end{aligned}
\quad (4.33) $$

In Figure 4.4, we plot the stability function and the internal stages of the double adjoint (4.26) of Chebyshev (4.11) and the double adjoint (4.32) of RKC (4.21). Comparing with Figures 4.1b and 4.2b, we observe that the internal stages are not the same for the double adjoint methods compared to the (4.11) and (4.21), while the stability function itself is identical as shown in Theorem 4.3.2 for a general Runge-Kutta method.
### 4.3.4 Stability and convergence analysis

In this section, we study the stability of the double adjoint of the Chebyshev method (4.25) and the RKC method (4.31). We recall that for the Chebyshev method of order one (resp. RKC of order two) the stability domain contains the interval $[-\beta(s, \eta), 0]$ where $\beta_{\text{Cheb}}(s, \eta) \approx (2 - 4\eta/3)s^2$ (resp. $\beta_{\text{RKC}}(s, \eta) \approx 0.653s^2$ for $\eta = 0.15$).

**Theorem 4.3.6.** Consider the Chebyshev (4.25) and the RKC (4.31) methods. For $\eta = 0$, the stability functions of the internal stages $R_{s,i}(z)$ of the Chebyshev (resp. RKC) double adjoint (4.26) (resp. (4.32)), are bounded by 1 for all $z \in [-2s^2, 0]$ (resp. $[-\frac{2}{3}s^2 + \frac{2}{3}, 0]$) and all $s \in \mathbb{N}$.

The proof of Theorem 4.3.6 relies on the following lemma.

**Lemma 4.3.7.** Let $s \ge 1$, and consider the double sequence $\tilde{\gamma}_j^i$ indexed by $i$ and $j$,

$$ 
\begin{aligned}
\tilde{\gamma}_j^i &= 0 \quad \forall j > i, & i &= 0, \dots, s-1, \\
\tilde{\gamma}_0^0 &= 1, & \tilde{\gamma}_0^1 &= 0, & \tilde{\gamma}_1^1 &= 2,
\end{aligned}
\qquad (4.34) $$

$$ 
\begin{aligned}
\tilde{\gamma}_0^i &= \tilde{\gamma}_1^{i-1} - \tilde{\gamma}_0^{i-2}, & \tilde{\gamma}_1^i &= 2\tilde{\gamma}_0^{i-1} + \tilde{\gamma}_2^{i-1} - \tilde{\gamma}_1^{i-2} & i &= 2, \dots, s-1, \\
\tilde{\gamma}_j^i &= \tilde{\gamma}_{j-1}^{i-1} + \tilde{\gamma}_{j+1}^{i-1} - \tilde{\gamma}_j^{i-2} & i &= 2, \dots, s-1, j = 2, \dots, i,
\end{aligned}
\qquad (4.35) $$

Then,

$$ 
\begin{aligned}
\tilde{\gamma}_j^i &= 0 \quad \forall j > i, \\
\tilde{\gamma}_0^i &= \begin{cases} 1 & \text{if } i \text{ is even} \\ 0 & \text{if } i \text{ is odd} \end{cases} \qquad \tilde{\gamma}_j^i = \begin{cases} 2 & \text{if } j - i \text{ is even} \\ 0 & \text{otherwise} \end{cases} \qquad (4.36)
\end{aligned}
$$

where $i = 0, \dots, s-1$, $j = 0, \dots, i$.

*Proof.* It can be checked that the coefficients defined in (4.36) verify the induction (4.35). Hence using the fact that they have the same initial terms (4.34), we conclude that they coincide by induction on $i$ and $j$. $\square$

*Proof of Theorem 4.3.6.* We first consider the Chebyshev method without damping applied to the linear test problem $y' = \lambda y$, $\lambda \in \mathbb{C}$, $t \in (0, T]$, $y(0) = 1$, with a uniform subdivision $x_0 = 0 < x_1 < \dots < x_N = T$ of stepsize $h$. Using Remark 4.3.4, we obtain for $k = 0$:

$$ 
\begin{aligned}
y_{k0} &= 1, & y_{k1} &= y_{k0} + \frac{h\lambda}{s^2}y_{k0}, \\
y_{ki} &= \frac{2ih\lambda}{(i+1)s^2}y_{k,i-1} + \frac{2i}{i+1}y_{k,i-1} + \frac{1-i}{i+1}y_{k,i-2}, & i &= 2, \dots, s-1, \\
y_{ks} &= \frac{h\lambda}{s}y_{k,s-1} + sy_{k,s-1} + (1-s)y_{k,s-2}, \\
y_1 &= y_{ks}.
\end{aligned}
\qquad (4.37) $$

First, notice that since $\frac{2i}{i+1} + \frac{1-i}{i+1} = 1$, we have that for all $i = 0, \dots, s$, $y_{ki} = (1 + O(h))$ (by induction). Setting $z = h\lambda$, it is sufficient to prove the identity

$$ y_{ki} = \sum_{j=0}^{i} \gamma_{j}^{i} T_{j} \left(1 + \frac{z}{s^{2}}\right) \qquad (4.38) $$
where $y_{ki}$ is a convex combination of the polynomials $T_j(1 + \frac{z}{s^2})$,

$$ \sum_{j=0}^{i} \gamma_{j}^{i} = 1 \text{ and } \gamma_{j}^{i} \geq 0 \ \forall \ i, j = 1, \dots, s-1, \qquad (4.39) $$

because $|T_j(1 + \frac{z}{s^2})| \le 1$ for all $j = 0, \dots, s$ and $z \in [-\beta(s, 0), 0] = [-2s^2, 0]$.

Since the Chebyshev polynomials form a basis of the vector space of polynomials, this already justifies the existence of the expansion (4.38) with some real coefficients $\gamma_j^i$. The identity $\sum_{j=0}^i \gamma_j^i = 1$ follows from the fact that $y_{ki} = 1 + O(h)$ and $T_j(1 + \frac{z}{s^2}) = 1 + O(h)$ for all $i, j = 0, \dots, s$. Now we can calculate these coefficients for the first two internal stages, $R_{s,0}(z) = y_{k0} = y_k = 1 = T_0(1 + \frac{z}{s^2})$, thus $\gamma_0^0 = 1$. Analogously,

$$ R_{s,1}(z) = y_{k1} = y_{k0} + \frac{h\lambda}{s^2}y_{k0} = 1 + \frac{z}{s^2} = T_1(1 + \frac{z}{s^2}) \text{ we obtain } \gamma_0^1 = 0, \gamma_1^1 = 1. $$

It remains to prove the positivity of the coefficients $\gamma_j^i$. Coupling (4.38) and (4.39), we obtain

$$ \begin{aligned} R_{s,i}(z) &= y_{ki} = \frac{2i}{i+1}\left(1+\frac{z}{s^2}\right)y_{k,i-1} + \frac{1-i}{i+1}y_{k,i-2} \\ &= \frac{2i}{i+1}\left(1+\frac{z}{s^2}\right)\sum_{j=0}^{i-1}\gamma_j^{i-1}T_j\left(1+\frac{z}{s^2}\right) + \frac{1-i}{i+1}\sum_{j=0}^{i-2}\gamma_j^{i-2}T_j\left(1+\frac{z}{s^2}\right) \\ &= \frac{2i}{i+1}\gamma_0^{i-1}T_1\left(1+\frac{z}{s^2}\right) + \sum_{j=2}^{i}\frac{i}{i+1}\gamma_{j-1}^{i-1}T_j\left(1+\frac{z}{s^2}\right) \\ &\quad + \sum_{j=2}^{i}\left(\frac{i}{i+1}\gamma_{j-1}^{i-1} - \frac{i-1}{i+1}\gamma_{j-2}^{i-2}\right)T_{j-2}\left(1+\frac{z}{s^2}\right) \end{aligned} $$

where we used (4.14). By comparison with (4.38), we obtain the following recurrence

$$ \begin{aligned} \gamma_j^i &= 0 \quad \forall j > i, \qquad i = 0, \dots, s-1, \\ \gamma_0^0 &= 1, \quad \gamma_1^0 = 0, \quad \gamma_1^1 = 1, \quad \gamma_0^i = \frac{i}{i+1}\gamma_1^{i-1} - \frac{i-1}{i+1}\gamma_0^{i-2} \quad i = 2, \dots, s-1, \\ \gamma_1^i &= \frac{2i}{i+1}\gamma_0^{i-1} + \frac{i}{i+1}\gamma_2^{i-1} - \frac{i-1}{i+1}\gamma_1^{i-2} \quad i = 2, \dots, s-1, \\ \gamma_j^i &= \frac{i}{i+1}\gamma_{j-1}^{i-1} + \frac{i}{i+1}\gamma_{j+1}^{i-1} - \frac{i-1}{i+1}\gamma_j^{i-2} \quad i = 2, \dots, s-1, \quad j = 2, \dots, i. \end{aligned} $$

Now defining $\tilde{\gamma}_j^i = (i+1)\gamma_j^i$, the above induction relations simplify to (4.34) and (4.35). The positivity of $\tilde{\gamma}_j^i$, and hence of $\gamma_j^i$, follows from Lemma 4.3.7. For $i = s$, the stability is a consequence of Theorem 4.3.2.
Analogously, the RKC method reads for $\eta = 0$,

$$ 
\begin{aligned}
y_{k0} &= 1, \quad y_{k1} = y_{k0} + \frac{3h\lambda}{s^2 - 1}y_{k0}, \\
y_{ki} &= \frac{6ih\lambda}{(i+1)(s^2-1)}y_{k,i-1} + \frac{2i}{i+1}y_{k,i-1} + \frac{1-i}{i+1}y_{k,i-2}, \quad i = 2, \dots, s-1, \\
y_{ks} &= \frac{h\lambda}{s}y_{k,s-1} + \frac{s^2-1}{3s}y_{k,s-1} - \frac{s^3-s^2-s+1}{3s^2}y_{k,s-2} + \frac{2s^2+1}{3s^2}y_0, \\
y_1 &= y_{ks},
\end{aligned} 
\quad (4.40) 
$$

and we follow the same methodology as above. Note that for RKC we have $|T_j(1+\frac{3}{s^2-1}z)| \le 1$ for all $j = 0, \dots, s$ and $z \in [-\beta(s, 0), 0] = [-\frac{2}{3}(s^2 - 1), 0]$. Using the same notations we search for coefficients satisfying the following

$$ y_{ki} = \sum_{j=0}^{i} \gamma_{j}^{i} T_{j} \left(1 + \frac{3}{s^2 - 1} z\right), \quad \text{where} \quad \sum_{j=0}^{i} \gamma_{j}^{i} = 1 \text{ and } \gamma_{j}^{i} \ge 0 \quad \forall i, j = 1, \dots, s-1. \quad (4.41) $$

Remark that $R_{s,0}(z) = y_{k0} = y_k = 1 = T_0(1 + \frac{3}{s^2-1}z)$, hence $\gamma_0^0 = 1$. Analogously, $R_{s,1}(z) = y_{k1} = y_{k0} + \frac{3h\lambda}{s^2-1}y_{k0} = 1 + \frac{3z}{s^2-1} = T_1(1+\frac{3}{s^2-1}z)$, and we deduce that $\gamma_1^1 = 0$, $\gamma_1^1 = 1$. Again we find a relation between these new coefficients to prove their positivity using $y_{ki} = \frac{2i}{i+1}(1 + \frac{3}{s^2-1}z)y_{k,i-1} + \frac{1-i}{i+1}y_{k,i-2}$. We get a recurrence of the same form as in the Chebyshev double adjoint method (4.37) but with different parameter, proceeding in the same we obtain exactly the same coefficients $\gamma_j^i$, which concludes the proof. $\square$

**Remark 4.3.8.** For the case of positive damping $\eta > 0$, the coefficients get very complicated and it is difficult to find a recurrence relation between them in order to prove their positivity. However, observing that all the coefficients in the recurrence relations of the internal stages of the methods are continuous functions of $\eta$, then for all $s$, there exists $\eta_0(s)$ such that the internal stages are stable (bounded) for all $\eta \in [0, \eta_0(s)]$. Numerical investigations suggest that Theorem 4.3.6 remains valid for all $\eta > 0$ i.e the methods remain stable with the stability functions of the internal stages bounded by 1, for all integers $s \ge 1$ for Chebyshev and $s \ge 2$ for RKC, and all $\eta > 0$. We have verified this numerically for $s \le 200$.

We conclude this section by the following convergence theorem for the new explicit stabilized methods for stiff optimal control problems.

**Theorem 4.3.9.** The method (4.25)-(4.26) (resp. (4.31)-(4.32)) has order 1 (resp. 2) of accuracy for the optimal control problem (4.1).

*Proof.* The proof follows immediately from Theorem 4.2.1 with $p_{oc} = p_{ode} = 2$ for the RKC method. $\square$

**Remark 4.3.10.** The proposed explicit stabilized integrators for optimal control problems could be combined with the idea of implicit-explicit (IMEX) integrators as proposed in [30], where RKC type methods would replace the implicit part in the IMEX integrator. This idea is already proposed in [58, 10] in the context of advection-diffusion-reaction problems. In [58], the diffusion part is discretized with an RKC method which typically has a large number of internal stages, and the advection-reaction part is integrated using a 4-stage explicit
Runge-Kutta method. In [10], the method integrates the diffusion term using ROCK2 method, the advection term using a 3-stage explicit method, and the nonlinear reaction term is solved implicitly. Such an extension is however out of the scope of the thesis.

## 4.4 Numerical experiments

In this Section, we illustrate numerically our theoretical findings of convergence and stability of the new fully explicit methods for stiff optimal control problems, first on a stiff three dimensional problem and second on a nonlinear advection-diffusion PDE (Burgers equation).

### 4.4.1 A linear quadratic stiff test problem

We start this section by a simple test problem taken from [27]:

$$ \begin{gathered} \min \frac{1}{2} \int_{0}^{1} (u^{2}(t) + 2x^{2}(t))dt \quad \text{subject to} \\ \dot{x}(t) = \frac{1}{2}x(t) + u(t), \quad t \in [0, 1], \quad x(0) = 1. \end{gathered} \qquad (4.42) $$

The optimal solution $(u^*, x^*)$ is given by $u^*(t) = \frac{2(e^{3t}-e^3)}{e^{3t/2}(2+e^3)}$, $x^*(t) = \frac{2e^{3t}+e^3}{e^{3t/2}(2+e^3)}$. As studied in [30] we modify problem (4.42) into a singularly perturbed (stiff) problem to illustrate the good stability properties of our new method. For a fixed $\varepsilon > 0$, we consider the following stiff optimal control problem,

$$ \begin{gathered} \min c(1) \quad \text{subject to} \\ \begin{aligned} \dot{c}(t) &= \frac{1}{2}(u^2(t) + x^2(t) + 4z^2(t)), \quad c(0) = 0, \\ \dot{x}(t) &= z(t) + u(t), \quad x(0) = 1, \\ \dot{z}(t) &= \frac{1}{\varepsilon}\left(\frac{1}{2}x(t) - z(t)\right), \quad z(0) = \frac{1}{2}, \end{aligned} \end{gathered} \qquad (4.43) $$

Figure 4.5 shows the convergence behavior, using the new RKC method (4.31)-(4.32), of the error in infinity norm between the solutions of the stiff problem (4.43) for $\varepsilon = 10^{-1}$ and $\varepsilon = 10^{-3}$ and different sizes of the time step $h_i = 2^{-i}$, $i = 0, \dots, 5$ and the reference solution is obtained with $h = 2^{-7}$. We observe lines of slope 2 which confirms the theoretical order two of accuracy of the scheme (Theorem 4.3.9). In the stiff case ($\varepsilon = 10^{-3}$), the method uses $s = 4$ to calculate the reference solution and $s = 40, 28, 20, 14, 10, 7$ respectively for the different time steps used to illustrate the convergence, these values coincide with the theoretical values that can be obtained using (4.18). Analogously to the case of stiff ODEs, the cost of scheme (4.31)-(4.32) is $O(\varepsilon^{-\frac{1}{2}})$ function evaluations of $f$, while using Euler method with its double adjoint would cost $O(\varepsilon^{-1})$.
<figure><img src="image_43.png" /><figcaption>Figure 4.5: Convergence plot of RKC (4.31)-(4.32) applied to problem (4.43).</figcaption></figure>

### 4.4.2 Optimal control of Burgers equation

To illustrate the performance of the new method, we consider the following optimal control problem of a nonlinear advection-diffusion PDE corresponding to the Burgers equation

$$
\begin{array}{ll}
\displaystyle \min_{u \in L^2([0,T]; L^2(\Omega))} & J(u) = \frac{1}{2} \|y(T) - y^{\text{target}}\|_{L^2(\Omega)}^2 + \frac{\alpha}{2} \int_0^T \|u(t)\|_{L^2(\Omega)}^2 \\
& \text{subject to} \\
& \partial_t y(t,x) = \mu \Delta y(t,x) - \frac{\nu}{2} \partial_x (y^2(t,x)) + u(t,x) \quad \text{in } (0,T) \times \Omega, \\
& y(0,x) = g(x) \quad \text{in } \Omega, \\
& y(t,x) = 0 \quad \text{on } \partial\Omega,
\end{array}
\qquad (4.44)
$$

where $\mu, \nu > 0$, in dimension $d = 1$ with domain $\Omega = (0,1)$ and the final time is given by $T = 2.5$. Here the control $u$ is a part of the source that we want to adjust in order to achieve a given final state $y^{\text{target}}: \Omega \to \mathbb{R}$.

We use a standard central finite difference space discretization for the state equation, and the trapezoid rule to discretize in space the norm $L^2(\Omega)$. We consider $M + 2$ points in space $x_m = m\Delta x$, with grid mesh size $\Delta x = \frac{1}{M+1}$, and we denote by $y_m(t)$ the approximation to $y(t, x_m)$, and define the vector $Y(t) = (y_0(t), y_1(t), \dots, y_{M+1}(t)) \in \mathbb{R}^{M+2}$. Similar notations are used for $U$ and $P$. We obtain the following optimal control problem
semi-discretized in space,

$$ \min_{U(t) \in \mathbb{R}^{M+2}} \Psi(c(T), Y(T)) = \frac{1}{2(M+1)} \sum_{m=0}^{M+1} (y_m(T) - y^{\text{target}}(x_m))^2 + \alpha c(T) $$

subject to

$$ \dot{c}(t) = \frac{1}{2(M+1)} \sum_{m=0}^{M+1}' u_m^2(t), \quad c(0) = 0, \qquad (4.45) $$

$$ \dot{y}_m(t) = F_m(U(t), Y(t)) := \frac{\mu}{\Delta x^2} (y_{m+1} - 2y_m + y_{m-1}) - \frac{\nu}{4\Delta x^2} (y_{m+1}^2 - y_{m-1}^2) + u_m, $$

$$ y_m(0) = g(x_m), \quad m = 0, \dots, M+1, $$

where $m = 0, \dots, M+1$ and the primed sum denotes a normal sum where the first and the last term are divided by 2 and we define $y_0 = y_{M+1} = 0$ to take into account the homogeneous Dirichlet boundary conditions. The function $F: \mathbb{R}^{M+2} \to \mathbb{R}^{M+2}$ with components $F_m: \mathbb{R}^{M+2} \to \mathbb{R}$ is obtained from the standard central finite difference discretization of the right hand side of the state equation (4.44), and adapted to the boundary conditions. The corresponding adjoint system is

$$ \begin{aligned} \dot{p}_c(t) &= 0, & p_c(T) &= \nabla_c \Psi = \alpha, \\ \dot{P}(t) &= -\nabla_Y F(U(t), Y(t))P, \end{aligned} \qquad (4.46) $$

where $P$ is a vector of length $M+2$ containing the costate values $p_m$, $m = 0, \dots, M+1$. In all our experiments we take $\mu = 0.1$, $\nu = 0.02$, $g(x) = \frac{3}{2}x(1-x)^2$, and $y^{\text{target}}(x) = \frac{1}{2}\sin(10x)(1-x)$.

In Figure 4.6 we plot the optimal control function (Fig. 4.6b) and the corresponding state function (Fig. 4.6a) obtained using scheme (4.31)-(4.32). When we use a small value for $\alpha$ in the model, we allow larger control values and thus a final state very close to the target (Fig. 4.6c), otherwise the control will be more limited and then the final state will not be that close to the target (Fig. 4.6d). Note that the method required $s = 24$ stages for $\Delta x = 1/100$ and $\Delta t = 2.5/30$, while using an Euler method with its double adjoint would require $\Delta t \le \Delta t_{\text{max,Euler}} := \Delta x^2/2$ at most. Hence, the standard Euler method would be $\Delta t/(s\Delta t_{\text{max,Euler}}) \approx 70$ times more expensive in terms of number of function evaluations for $\Delta x = 1/100$, a factor that grows arbitrarily as $\Delta x \to 0$.

In Figure 4.7, we plot the convergence curves for the state and control functions of the new RKC method (4.31)-(4.32) applied to the diffusion problem discretized in space (4.45), where the number of stages $s$ is computed adaptively using (4.18). We recover again lines of slope two, which corroborate the order two of the method. Although our convergence analysis is valid only in finite dimensions (Theorem 4.3.9), this suggests that the convergence of order two persists in the PDE case. For comparison, we also included the results for the following standard diagonally implicit Runge-Kutta method of order two, inspired from [30, Table 5.1] in the context of stiff optimal control problems, and given by the following Butcher tableau where $\gamma = 1 - \sqrt{2}/2$ (making the method L-stable),

$$ \begin{array}{c|cc} \gamma & & \\ 1 - 2\gamma & \gamma & \\ \hline 1/2 & 1/2 & \end{array} \qquad (4.47) $$
<figure><img src="image_44.png" /><figcaption>Figure 4.6: State, final state, and control, of problem (4.44). Figures (a), (b), and (c) are obtained using $\Delta x = 1/100$, $\Delta t = T/30$, and $s = 24$ stages, and $\alpha = 0.01$. Figure (d) uses the same $\Delta x$, $g$ and $y^{\text{target}}$ but $\alpha = 0.02$.</figcaption></figure>

Although for a fixed timestep, the second order implicit method (4.47) appears about two times more accurate than the RKC method (4.31)-(4.32) for the control and almost of the same accuracy for the state, we emphasize that these convergence plots do not take into account the extra cost of the implicitness of method (4.47). Indeed, the cost and difficulty of the implementation of the implicit methods (nonlinear iterations, preconditioners, etc.) would typically deteriorate in larger dimensions and for a nonlinear diffusion operator, as it is already the case for initial value PDEs [4], while as an explicit stabilized scheme, the RKC method (4.31)-(4.32) can be conveniently implemented in the spirit of the simplest explicit Euler method.
<figure><img src="image_45.png" /><figcaption>Figure 4.7: Convergence plot of the RKC method (4.31)-(4.32) and the implicit method (4.47) applied to problem (4.45) for many time steps $\Delta t_i = T/2^i$, $i = 3, \dots, 8$, $\Delta x = 1/100$, and $\alpha = 0.02$. The reference solution is obtained using $\Delta t = T/2^{12}$, $s = 3$.</figcaption></figure>
Chapter 5

Conclusion and outlook

We have designed new and explicit stabilized schemes for different types of problems, and we think that we now have enough material to look at many new interesting ideas related to our work. Before concluding we would like to present one idea that we looked at, and we think it is very interesting and promising.

5.1 Towards explicit implementation of implicit methods using optimization techniques and explicit stabilized integrators

We aim to implement excellent well established higher order implicit Runge-Kutta methods explicitly by combining some optimization techniques with an explicit stabilized method of order 1. We will explain the methodology in the rest of the current section.

In the recent paper [24], the authors proposed an explicit stabilized version of gradient descent (GD) to replace the standard GD method (5.3) for stiff optimization problems in large dimension to avoid step size restriction. Here we mean by stiff optimization problems the case where the gradient of the objective function is stiff. In order to solve the problem

$$
\min_{x \in \mathbb{R}^d} f(x), \qquad (5.1)
$$

where $d \in \mathbb{N} \setminus \{0\}$ and $f: \mathbb{R}^d \to \mathbb{R}$ a continuously differentiable function, with $L$-Lipschitz gradient, and such that the real parts of the eigenvalues of its Hessian matrix $\mathcal{H}f(x)$ are strictly positive for all $x \in \mathbb{R}^d$, the authors of [24] consider its gradient flow

$$
\dot{x} = -\nabla f(x), \quad x(0) = x_0 \in \mathbb{R}^d, \qquad (5.2)
$$

where $x_0$ is an initialization. Instead of applying an explicit Euler method to (5.2), which is equivalent to GD, they apply an explicit stabilized method. Indeed, if we discretize (5.2) using an explicit Euler method with step size $h$, we get the famous gradient descent (GD) method

$$
x_{n+1} = x_n - h \nabla f(x_n), \quad n = 0, 1, 2, \dots \tag{5.3}
$$
where $x_n$ is the numerical approximation of $x(t_n)$. For stiff gradients, this method faces severe restriction on $h$ in order to be stable. Another approach is to discretize (5.2) using an implicit Euler method with a step size $h$

$$x_{n+1} = x_n - h\nabla f(x_{n+1}), \quad n = 0, 1, 2, \dots \qquad (5.4)$$

It is easy to verify that $x_{n+1}$ in (5.4) is the unique minimizer of the problem

$$\min_{x \in \mathbb{R}^d} hf(x) + \frac{1}{2}\|x - x_n\|_2^2. \qquad (5.5)$$

It is reasonable to think of applying a GD method with step size $\tau$ to the gradient flow of (5.5) to avoid implicitness, and reduce the step size restriction, but this will not reduce enough the restriction on $\tau$ especially if $h$ is large. This pushes us to think of applying an explicit stabilized method to (5.5) with a step size $\tau = \mathcal{O}(1)$ and a suitable number of stages $s$, and iterate until convergence. The obtained algorithm is in fact explicit. We look at our stiff ODE as a gradient flow of some optimization problem, for which, the objective function needs not to be known. We aim to apply this idea to a powerful implicit method of higher order and perfect stability properties such as RADAU methods. If it works with the fifth order RADAU IIA for example, we obtain an explicit implementation of an excellent stiffly accurate method of order 5.

In a general framework, let $d$ be a positive integer and consider the ODE

$$\dot{y} = f(y), \quad y(0) = y_0, \qquad (5.6)$$

where $y_0 \in \mathbb{R}^d$, and $f: \mathbb{R}^d \to \mathbb{R}^d$ is an $L$-Lipschitz vector field. Suppose that the eigenvalues of the Jacobian matrix of $f$ have strictly negative real parts for all $y \in \mathbb{R}^d$. Consider a stiffly accurate implicit method (for example RADAU IIA of order 5) with $m$ stages, and with matrix of coefficients $\mathcal{A}$. After some calculations we arrive at the minimization problem

$$\min_{\mathcal{Y} \in \mathbb{R}^{d \times m}} hF(\mathcal{Y}) + \frac{1}{2} \|\mathcal{Y} - \mathcal{Y}_0\|_2^2, \qquad (5.7)$$

where $\mathcal{Y}$ is a vector containing all the internal stages of the considered method. The function $F$ needs not to be known and it is such that $-\nabla F(\mathcal{Y}) = \mathcal{A} \otimes I_d \begin{pmatrix} f(Y_1) \\ f(Y_2) \\ \vdots \\ f(Y_m) \end{pmatrix}$. The corresponding gradient flow is

$$\dot{\mathcal{Y}} = h\mathcal{A} \otimes I_d \begin{pmatrix} f(Y_1) \\ f(Y_2) \\ \vdots \\ f(Y_m) \end{pmatrix} - \mathcal{Y} + \mathcal{Y}_0, \quad \mathcal{Y}(0) = \mathcal{Y}_0. \qquad (5.8)$$

The idea now is to apply an explicit stabilized integrator to (5.8) with a time step $\tau = \mathcal{O}(1)$, and iterate.

Consider the test problem $\dot{y} = -\lambda y$, $y(0) = y_0$, where $\Re\lambda < 0$, which leads to

$$\dot{\mathcal{Y}} = -h\lambda\mathcal{A}\mathcal{Y} - \mathcal{Y} + \mathcal{Y}_0, \quad \mathcal{Y}(0) = \mathcal{Y}_0. \qquad (5.9)$$
In the case where the eigenvalues of A are real positive, we can apply the standard Chebyshev method (2.14). For some implicit methods (RADAU IIA for example) the eigenvalues of A are complex, which makes the standard Chebyshev method useless.

By diagonalizing A, we compute $A^r$ and $A^i$ such that $A = A^r + iA^i$ and the eigenvalues of $A^r$ are the real parts of the eigenvalues of A, and the eigenvalues of $A^i$ are the imaginary parts of the eigenvalues of A.

The above splitting of the matrix A requires an explicit stabilized method that is stable enough in the imaginary direction, and this can be done by stabilizing the part of the vector field with almost pure imaginary eigenvalues (advection-like) using a new polynomial. The overall stability function will look like

$$Q_s(h\tau\lambda A^r) + B_s(h\tau\lambda A^r)(h\tau\lambda iA^i). \qquad (5.10)$$

For optimization purposes, we ask that the above partition preserves the steady state of the system in both linear and nonlinear vector fields cases. We consider

$$Q_s(p) = \frac{T_s(\omega_0 + \omega_1' p)}{2T_s(\omega_0)} + \frac{1}{2}$$

and $B_s(p) = (Q_s(p) - 1)/p$, where $\omega_0 = 1 + \eta/s^2$, $\omega_1' = \frac{2T_s(\omega_0)}{T_s'(\omega_0)}$ and $\eta = 0.1$, and the overall stability function

$$R_s(p, q) = \frac{T_s(\omega_0 + \omega_1' p)}{2T_s(\omega_0)} + \frac{1}{2} + \frac{1}{p} \left( \frac{T_s(\omega_0 + \omega_1' p)}{2T_s(\omega_0)} - \frac{1}{2} \right) iq, \qquad (5.11)$$

where $p = h\tau\lambda A^r$ and $q = h\tau\lambda A^i$. To motivate the above choice, consider the test equation $\dot{y} = Uy + Vy$, $y(0) = y_0$, such that $Uy_0 + Vy_0 = 0$, we want $y_1 = y_0$, i.e.

$$\begin{aligned} Q_s(hU)y_0 + B_s(hU)hVy_0 &= y_0, \\ -B_s(hU)hVy_0 &= (Q_s(hU) - I)y_0, \\ B(hU)hUy_0 &= (Q_s(hU) - I)y_0, \quad \text{since } -Vy_0 = Uy_0, \\ B_s(p) &= (Q_s(p) - 1)/p, \end{aligned}$$

where $p = hU$.

For problems of the form

$$\dot{y} = f(y) + g(y), \quad y(0) = y_0, \qquad (5.12)$$

where $g$ contains the imaginary parts of the eigenvalues appearing in the case of RADAU IIA for example, we consider the following integrator

$$\begin{aligned} K_0 &= y_0, & K_1 &= y_0 + \frac{\omega_1'}{2\omega_0}h(f(y_0) + g(y_0)), \\ K_j &= \mu_j h(f(K_{j-1}) - \frac{1}{2}f(y_0) + \frac{1}{2}g(y_0)) + \nu_j K_{j-1} + (1 - \nu_j)K_{j-2}, \quad j = 2, \dots, s, \\ y_1 &= K_s, \end{aligned} \qquad (5.13)$$

where,

$$\mu_j = \frac{2\omega_1' T_{j-1}(\omega_0)}{T_j(\omega_0)}, \quad \nu_j = \frac{2\omega_0 T_{j-1}(\omega_0)}{T_j(\omega_0)}.$$
**Lemma 5.1.1.** Applied to the linear test problem $\dot{y} = \lambda y + i\sigma y$ the scheme (5.13) yields the following stability function

$$R_s(p, q) = \frac{T_s(\omega_0 + \omega_1' p)}{2T_s(\omega_0)} + \frac{1}{2} + \frac{1}{p} \left( \frac{T_s(\omega_0 + \omega_1' p)}{2T_s(\omega_0)} - \frac{1}{2} \right) iq, \quad (5.14)$$

where $p = h\lambda$ and $q = h\sigma$.

**Theorem 5.1.2.** The method (5.13) is of order 1 of accuracy for ODEs of the form (5.12), with $f$ and $g$ are Lipschitz continuous.

**Proposition 5.1.3** (steady state conservation). Let $\dot{y} = f(y) + g(y)$ and let $y(0) = y_0$ such that $f(y_0) + g(y_0) = 0$. Then, $y_1 = y_0$, with $y_1$ defined as in (5.13).

Proofs of Lemma 5.1.1 and Theorem 5.1.2 are similar to the proofs of Lemma 3.3.1 and Theorem 3.4.1 respectively. The proof of Proposition 5.1.3 is straightforward by induction.

## 5.2 Conclusion

In this thesis, we have constructed novel and efficient explicit stabilized methods for different types of stiff problems by combining techniques from stiff integration and geometric numerical integration. Part of the work is presented in two research articles [5, 14].

A remarkable feature of the explicit stabilized schemes proposed in Chapter 3 (stiff stochastic problems) and Chapter 4 (stiff optimal control problems) is that their extended size of stability domain is not only optimally large in many situations, but it is also based on rigorous estimates and it does not rely on empirical estimates as proposed in the past literature. This asset makes the new schemes promising in the context of variance reduction techniques (for instance combined with the multilevel Monte-Carlo method for stiff stochastic problems [6]) or parallel computing (for instance to be combined with the parareal algorithm for dissipative problems, in particular for optimal control [43]), where the ability of the integrators to apply large time steps with reliable stability is a key ingredient.

The SK-ROCK method with optimal stability domain presented in Chapter 3 is used in the article [47] to construct a highly efficient proximal Markov chain Monte Carlo methodology to perform Bayesian computation in imaging problems. Instead of the conventional Euler-Maruyama approximation that underpins existing proximal Monte Carlo methods, the authors use SK-ROCK method to significantly accelerate the convergence speed, similarly to accelerated gradient optimization methods. A multirate version of SK-ROCK is also used in the thesis [48] to solve stiff multirate stochastic differential equations.

Recently, the idea of stabilized Chebyshev methods was used in [21] in the context of ODEs arising from the discretization of wave equations. Based on Chebyshev polynomials, and in the spirit of explicit stabilized methods from the literature, the authors construct a stabilized version of leapfrog method for linear and semilinear second-order differential equations.

To conclude, we see that explicit stabilized methods are very useful in many contexts and not only for stiff dissipative problems, and we believe that they could still be used for many other classes of stiff problems and hence would give rise to many new interesting research ideas.
Bibliography

[1] A. Abdulle. On roots and error constants of optimal stability polynomials. *BIT Numerical Mathematics*, 40(1):177–182, 2000.

[2] A. Abdulle. *Chebyshev methods based on orthogonal polynomials*. PhD Thesis, University of Geneva, Department of Mathematics. University of Geneva, 2001.

[3] A. Abdulle. Fourth order Chebyshev methods with recurrence relation. *SIAM J. Sci. Comput.*, 23(6):2041–2054, 2002.

[4] A. Abdulle. *Explicit Stabilized Runge-Kutta Methods*, pages 460–468. Encyclopedia of Applied and Computational Mathematics, Springer Berlin Heidelberg, 2015.

[5] A. Abdulle, I. Almuslimani, and G. Vilmart. Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations. *SIAM/ASA J. Uncertain. Quantif.*, 6(2):937–964, 2018.

[6] A. Abdulle and A. Blumenthal. Stabilized multilevel Monte Carlo method for stiff stochastic differential equations. *J. Comput. Phys.*, 251:445–460, 2013.

[7] A. Abdulle and S. Cirilli. S-ROCK: Chebyshev methods for stiff stochastic differential equations. *SIAM J. Sci. Comput.*, 30(2):997–1014, 2008.

[8] A. Abdulle and T. Li. S-ROCK methods for stiff Ito SDEs. *Commun. Math. Sci.*, 6(4):845–868, 2008.

[9] A. Abdulle and A. Medovikov. Second order chebyshev methods based on orthogonal polynomials. *Numer. Math.*, 90(1):1–18, 2001.

[10] A. Abdulle and G. Vilmart. PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise. *J. Comput. Phys.*, 242:869–888, 2013.

[11] A. Abdulle, G. Vilmart, and K. C. Zygalakis. Mean-square A-stable diagonally drift-implicit integrators of weak second order for stiff Itô stochastic differential equations. *BIT*, 53(4):827–840, 2013.
[12] A. Abdulle, G. Vilmart, and K. C. Zygalakis. Weak second order explicit stabilized methods for stiff stochastic differential equations. *SIAM J. Sci. Comput.*, 35(4):A1792–A1814, 2013.

[13] G. Albi, M. Herty, and L. Pareschi. Linear multistep methods for optimal control problems and applications to hyperbolic relaxation systems. *Appl. Math. Comput.*, 354:460–477, 2019.

[14] I. Almuslimani and G. Vilmart. Explicit stabilized integrators for stiff optimal control problems. *Under revision for SIAM J. Sci. Comput.*, arXiv: 1910.10584, 2020.

[15] M. Bakker. Analytical aspects of a minimax problem. 1971. Technical Note TN 62 (in Dutch), Mathematical centre, Amsterdam.

[16] J. F. Bonnans and J. Laurent-Varin. Computation of order conditions for symplectic partitioned Runge-Kutta schemes with application to optimal control. *Numer. Math.*, 103(1):1–10, 2006.

[17] C.-E. Bréhier and G. Vilmart. High Order Integrator for Sampling the Invariant Distribution of a Class of Parabolic Stochastic PDEs with Additive Space-Time Noise. *SIAM J. Sci. Comput.*, 38(4):A2283–A2306, 2016.

[18] E. Buckwar and C. Kelly. Towards a systematic linear stability analysis of numerical methods for systems of stochastic differential equations. *SIAM Journal on Numerical Analysis*, 48(1):298–321, 2010.

[19] K. Burrage, P. Burrage, and T. Tian. Numerical methods for strong solutions of stochastic differential equations: an overview. *Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci.*, 460(2041):373–402, 2004.

[20] J. C. Butcher. The effective order of Runge-Kutta methods. In J. L. Morris, editor, *Proceedings of Conference on the Numerical Solution of Differential Equations*, volume 109 of *Lecture Notes in Math.*, pages 133–139, 1969.

[21] C. Carle, M. Hochbruck, and A. Sturm. On Leapfrog-Chebyshev Schemes. *SIAM J. Numer. Anal.*, 58(4):2404–2433, 2020.

[22] Y. Chong and J. B. Walsh. The roughness and smoothness of numerical solutions to the stochastic heat equation. *Potential Anal.*, 37(4):303–332, 2012.

[23] A. M. Davie and J. G. Gaines. Convergence of numerical schemes for the solution of parabolic stochastic partial differential equations. *MATH. COMP*, 70:121–134, 2000.

[24] A. Eftekhari, B. Vandereycken, G. Vilmart, and K. C. Zygalakis. Explicit stabilised gradient descent for faster strongly convex optimisation. *To appear in BIT Numerical Mathematics (2020)*.

[25] T. Gard. *Introduction to stochastic differential equations*. Marcel Dekker, New York, 1988.

[26] M. B. Giles. Multilevel Monte Carlo path simulation. *Oper. Res.*, 56(3):607–617, 2008.
[27] W. W. Hager. Runge-Kutta methods in optimal control and the transformed adjoint system. *Numer. Math.*, 87(2):247-282, 2000.

[28] E. Hairer, C. Lubich, and G. Wanner. *Geometric numerical integration*, volume 31 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, second edition, 2006. Structure-preserving algorithms for ordinary differential equations.

[29] E. Hairer and G. Wanner. *Solving ordinary differential equations II. Stiff and differential-algebraic problems*. Springer-Verlag, Berlin and Heidelberg, 1996.

[30] M. Herty, L. Pareschi, and S. Steffensen. Implicit-explicit Runge-Kutta schemes for numerical discretization of optimal control problems. *SIAM J. Numer. Anal.*, 51(4):1875-1899, 2013.

[31] M. Herty and V. Schleper. Time discretizations for numerical optimisation of hyperbolic problems. *Appl. Math. Comput.*, 218(1):183-194, 2011.

[32] D. Higham. An algorithmic introduction to numerical simulation of stochastic differential equations. *SIAM Review*, 43(3):525-546, 2001.

[33] D. J. Higham. Mean-square and asymptotic stability of the stochastic theta method. *SIAM J. Numer. Anal.*, 38(3):753-769, 2000.

[34] W. Hundsdorfer and J. Verwer. *Numerical solution of time-dependent advection-diffusion-reaction equations*, volume 33 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 2003.

[35] C. Y. Kaya. Inexact restoration for Runge-Kutta discretization of optimal control problems. *SIAM J. Numer. Anal.*, 48(4):1492-1517, 2010.

[36] P. Kloeden and E. Platen. *Numerical solution of stochastic differential equations*. Springer-Verlag, Berlin and New York, 1992.

[37] J. Lang and J. G. Verwer. W-methods in optimal control. *Numer. Math.*, 124(2):337-360, 2013.

[38] A. Laurent and G. Vilmart. Exotic aromatic B-series for the study of long time integrators for a class of ergodic SDEs. *Math. Comp.*, 89(321):169-202, 2020.

[39] B. Leimkuhler and C. Matthews. Rational construction of stochastic numerical methods for molecular sampling. *Appl. Math. Res. Express.*, 2013(1):34-56, 2013.

[40] B. Leimkuhler, C. Matthews, and M. V. Tretyakov. On the long-time integration of stochastic gradient systems. *Proc. R. Soc. A*, 470(2170), 2014.

[41] Q. Li, L. Chen, C. Tai, and W. E. Maximum principle based algorithms for deep learning. *Journal of Machine Learning Research*, 18(165):1-29, 2018.

[42] X. Liu and J. Frank. Symplectic runge-kutta discretization of a regularized forward-backward sweep iteration for optimal control problems. *Journal of Computational and Applied Mathematics*, 383:113133, 2021.
[43] Y. Maday, M.-K. Riahi, and J. Salomon. Parareal in time intermediate targets methods for optimal control problems. In *Control and optimization with PDE constraints*, volume 164 of *Internat. Ser. Numer. Math.*, pages 79–92. Birkhäuser/Springer Basel AG, Basel, 2013.

[44] G. N. Mil'shteǐn. A theorem on the order of convergence of mean-square approximations of solutions of systems of stochastic differential equations. *Teor. Veroyatnost. i Primenen.*, 32(4):809–811, 1987.

[45] G. Milstein. Weak approximation of solutions of systems of stochastic differential equations. *Theory Probab. Appl.*, 30(4):750–766, 1986.

[46] G. Milstein and M. Tretyakov. *Stochastic numerics for mathematical physics*. Scientific Computing. Springer-Verlag, Berlin and New York, 2004.

[47] M. Pereyra, L. Vargas Mieles, and K. C. Zygalakis. Accelerating proximal Markov chain Monte Carlo by using an explicit stabilized method. *SIAM J. Imaging Sci.*, 13(2):905–935, 2020.

[48] G. Rosilho De Souza. Numerical methods for deterministic and stochastic differential equations with multiple scales and high contrasts. page 216, 2020.

[49] Y. Saito and T. Mitsui. Stability analysis of numerical schemes for stochastic differential equations. *SIAM J. Numer. Anal.*, 33:2254–2267, 1996.

[50] J. M. Sanz-Serna. Symplectic Runge-Kutta schemes for adjoint equations, automatic differentiation, optimal control, and more. *SIAM Rev.*, 58(1):3–33, 2016.

[51] B. Sommeijer, L. Shampine, and J. Verwer. RKC: an explicit solver for parabolic PDEs. *J. Comput. Appl. Math.*, 88:316–326, 1998.

[52] B. P. Sommeijer and J. G. Verwer. *A performance evaluation of a class of Runge-Kutta-Chebyshev methods for solving semidiscrete parabolic differential equations*. Afdeling Numerieke Wiskunde [Department of Numerical Mathematics], 91. Mathematisch Centrum, Amsterdam, 1980.

[53] A. Tocino. Mean-square stability of second-order Runge-Kutta methods for stochastic differential equations. *J. Comput. Appl. Math.*, 175(2):355–367, 2005.

[54] P. J. van der Houwen and B. P. Sommeijer. On the internal stability of explicit, *m*-stage Runge-Kutta methods for large *m*-values. *Z. Angew. Math. Mech.*, 60(10):479–485, 1980.

[55] J. Verwer. Explicit Runge-Kutta methods for parabolic partial differential equations. *Special issue of Appl. Num. Math.*, 22:359–379, 1996.

[56] G. Vilmart. Postprocessed integrators for the high order integration of ergodic SDEs. *SIAM J. Sci. Comput.*, 37(1):A201–A220, 2015.

[57] A. Walther. Automatic differentiation of explicit Runge-Kutta methods for optimal control. *Comput. Optim. Appl.*, 36(1):83–108, 2007.
[58] C. J. Zbinden. Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems. *SIAM J. Sci. Comput.*, 33(4):1707–1725, 2011.
List of Figures

<table>
  <tbody>
    <tr>
      <td>1.1</td>
      <td>Different parts of the second moment of the stability function of SK-ROCK</td>
      <td>4</td>
    </tr>
    <tr>
      <td>1.2</td>
      <td>Deterministic complex stability domain for different damping parameters</td>
      <td>4</td>
    </tr>
    <tr>
      <td>1.3</td>
      <td>Stochastic mean-square stability domains</td>
      <td>5</td>
    </tr>
    <tr>
      <td>1.4</td>
      <td>Stability domains of the RKC method and the Heun method</td>
      <td>6</td>
    </tr>
    <tr>
      <td>1.5</td>
      <td>Internal stages and stability polynomial of the RKC method and its double adjoint</td>
      <td>9</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Graphical illustration of the approximation by Heun method after one step.</td>
      <td>13</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Stability domains of the explicit and the implicit Euler methods.</td>
      <td>15</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Stability domains and stability functions of the Chebyshev method for different damping values</td>
      <td>19</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Internal stages and stability polynomials of the Chebyshev method with and without damping</td>
      <td>19</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>Internal stages and stability polynomials of the RKC method and the ROCK2 method</td>
      <td>21</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>Stability domains of the RKC method and the Heun method</td>
      <td>21</td>
    </tr>
    <tr>
      <td>2.7</td>
      <td>A plot of 100 Brownian paths, their average, and their variance</td>
      <td>23</td>
    </tr>
    <tr>
      <td>2.8</td>
      <td>Exact solution of (2.33) and its approximation using Euler-Maruyama.</td>
      <td>26</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Stability domains and stability functions of the deterministic Chebyshev method for different damping values</td>
      <td>34</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Mean-square stability domains of the standard and new stochastic Chebyshev methods</td>
      <td>36</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Stability function of the new SK-ROCK method as a function of p with drift and diffusion contributions</td>
      <td>39</td>
    </tr>
    <tr>
      <td>3.4</td>
      <td>Strong and weak convergence plots using SK-ROCK of the nonlinear problem (3.52)</td>
      <td>50</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Weak convergence plots using SK-ROCK of the stiff nonlinear problem(3.53)</td>
      <td>51</td>
    </tr>
    <tr>
      <td>3.6</td>
      <td>Second moment error for the linear additive problem (3.55)</td>
      <td>52</td>
    </tr>
    <tr>
      <td>3.7</td>
      <td>PSK-ROCK without damping (η = 0)</td>
      <td>53</td>
    </tr>
    <tr>
      <td>3.8</td>
      <td>Second moment errors versus the average number of drift function evaluations for problem (3.56)</td>
      <td>54</td>
    </tr>
    <tr>
      <td>3.9</td>
      <td>SPDE problem (3.57) using the space discretization stepsize Δx = 1/100.</td>
      <td>56</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>3.10</td>
      <td>SPDE problem (3.57) with the initial condition $u(0, x) = 1$ and different damping parameters</td>
      <td>56</td>
    </tr>
    <tr>
      <td>3.11</td>
      <td>Stability domain of the new DA-ROCK method</td>
      <td>59</td>
    </tr>
    <tr>
      <td>3.12</td>
      <td>Solution and convergence plot of the DA-ROCK method applied to problem (3.65).</td>
      <td>60</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Internal stages and stability polynomials of the Chebyshev method with and without damping</td>
      <td>69</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Internal stages and stability polynomials of the classical and the new RKC implementations</td>
      <td>70</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Internal stages and stability polynomials of the ROCK2 method and its double adjoint</td>
      <td>71</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Internal stages and stability polynomials of the double adjoints of Chebyshev and RKC methods</td>
      <td>76</td>
    </tr>
    <tr>
      <td>4.5</td>
      <td>Convergence plot of RKC (4.31)-(4.32) applied to problem (4.43).</td>
      <td>81</td>
    </tr>
    <tr>
      <td>4.6</td>
      <td>State, final state, and control, of Burgers equation.</td>
      <td>83</td>
    </tr>
    <tr>
      <td>4.7</td>
      <td>Convergence plot of the RKC method (4.31)-(4.32) and the implicit method (4.47) applied to problem (4.45).</td>
      <td>84</td>
    </tr>
  </tbody>
</table>