Fast Numerical Solvers for Parameter Identification Problems in Mathematical Biology

TL;DR

This paper introduces efficient numerical methods for solving complex optimization problems in mathematical biology.

Contribution

The paper proposes commutative discretization and optimization techniques for PDE-constrained optimization in biology.

Findings

The proposed discretization ensures commutativity with optimization operations.

Efficient preconditioned iterative methods are applied to solve large-scale systems.

Numerical experiments confirm the approach's viability and efficiency.

Abstract

In this paper, we consider effective discretization strategies and iterative solvers for nonlinear PDE-constrained optimization models of pattern evolution within biological processes. Upon a Sequential Quadratic Programming linearization of the optimization problem, we devise appropriate time-stepping schemes and discrete approximations of the cost functionals such that the discretization and optimization operations are commutative, a highly desirable property of a discretization of such problems. We formulate the large-scale, coupled linear systems in such a way that efficient preconditioned iterative methods can be applied within a Krylov subspace solver. Numerical experiments demonstrate the viability and efficiency of our approach.

Figures

Snapshot of the solution at $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.7$$\end{document}$ for the desired states $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{u}$$\end{document}$ , $\documentclass[12pt]{minimal} \usepa

Full text

Introduction

To deepen our understanding of complex biological systems, it is crucial to develop mathematical models to uncover underlying mechanisms and analyze such systems. One of the most fundamental challenges in biology is to investigate how structures and shapes of organisms arise [1], or how an organism develops from, e.g., a single fertilized egg [2], which involves modelling and studying developmental processes. A crucial mechanism in embryonic development is the formation of spatial patterns that then evolve into formation of different tissues and organs. Modelling the mechanisms that give rise to these patterns is an active research area and a variety of structures have been studied, ranging from patterns on the skin that give rise to hair follicles [3], to the formation of organs such as the villi of the intestine [4] or the development of digits of a limb [5].

The most well-known theory for spatial pattern formation was first suggested by Alan Turing in his seminal paper [6]. Turing observed that asymmetric organisms develop from symmetric embryos, and the symmetry is broken during morphogenesis, the process of tissue and organ formation in an embryo. He proposed that this “symmetry-breaking” is caused by reactions between the chemicals involved, termed the morphogens, as their interactions can amplify small disturbances in the nearly homogeneous tissue. By adding diffusion to each morphogen at different rates, instability is added to the system, which gives rise to so-called Turing patterns. While there exist other theories for pattern formation, Turing’s reaction–diffusion theory is still a widely used modelling tool [7, 8].

The dynamics of morphogens can be modelled using a system of partial differential equations (PDEs). To analyse and simulate the system numerically, we first need to determine the parameters and source functions included in the equations. In an ideal world, all of these can be measured directly, however, in practice we are frequently required to fit them using indirect measurements. A parameter identification problem consists of solving an inverse problem. Assuming we have one or more unknown parameters (or functions) within a PDE model describing a biological system, along with one or more observations such as images, a valuable approach for determining or approximating the missing parameters is to formulate an optimization problem with these PDEs serving as constraints and the unknown parameters as controls. This is an example of a PDE-constrained optimization (PDECO) problem; see [9, 10] for comprehensive overviews of such problems. The optimization aspect consists of minimizing the deviation of solutions of the PDEs from the observation in combination with the amount of control required to do so.

In order to solve such a PDECO problem numerically, we derive a system of discrete optimality conditions in the form of discretized PDEs. While there are a variety of approaches for solving PDE systems, a key challenge when tackling time-dependent problems is that the coupled system involves PDEs evolving forward and backward in time. In order to achieve good convergence properties, we elect to solve the system with the all-at-once method, that is we simultaneously solve the linearized and discretized system of optimality conditions for all variables, across the entire time interval and within the whole spatial domain. Clearly, the dimensions of the linearized systems for fine mesh and time interval discretizations can be extremely large, which necessitates a highly efficient and robust solution strategy, with minimal computational storage requirements. We therefore solve the fully discretized system iteratively, and devise a suitable preconditioning strategy to accelerate convergence and allow us to obtain the solution fast and robustly. We highlight that the all-at-once approach for solving reaction–diffusion models, employing backward Euler discretization in time, was demonstrated in [11]. Another class of methods for addressing such problems are gradient-based techniques. Notable examples of image-driven PDECO problems for reaction–diffusion equations include [12], and [13–15] for spatial patterns in particular. Additionally, similar approaches have been applied to reaction–diffusion–advection equations [16].

The process of transforming the optimization problem into a fully discretized PDE system to be solved numerically can be done by performing the optimization and discretization steps in either order. Developing adjoint-consistent approaches that enable discretization and optimization to commute, thus making the system matrices resulting from the Optimize-then-Discretize (OTD) and Discretize-then-Optimize (DTO) approaches equivalent, is of great interest in the field. It is also of importance when one wishes to quantify the discretization error properties of such problems. In general, the DTO approach can be applied such that one obtains a symmetric system matrix, while the symmetry of the OTD matrix depends on the discretization scheme of choice, including the time discretization if applicable. We refer to [17] and [18, Ch. 3] for discussions of the two approaches, and [19, Ch. 2] for a derivation for Poisson control problems showing adjoint-consistency given suitable choices of finite element discretizations within each approach.

However, obtaining adjoint-consistent methods can be a challenging task even for relatively simple (linear) PDEs, e.g., applying the two approaches to the (steady) advection–diffusion equation with streamline upwind/Petrov Galerkin (SUPG) stabilization was shown to lead to different results [20], and this would extend to the time-dependent setting if the same spatial stabilization were applied. We refer to [21, Sec. 2.1] for a discussion of systems obtained for steady and time-dependent problems, including a demonstration that heat control problems with backward Euler discretization will generally lead to different systems for OTD and DTO. This may be mitigated if specific choices are taken to approximate the time-dependent cost functional, but doing so may affect the order of convergence. See also [22, 23] for derivation of an adjoint-consistent strategy for heat control problems, which is utilized within our methodology. Descriptions of the systems obtained for problems constrained by advection equations are given in [24]. In this work, we wish to develop a strategy that enables commutation between the discretization and optimization steps for challenging nonlinear (in particular semilinear) PDE systems arising from modelling pattern formation. Inspired by the Störmer–Verlet integrator (see [25] for a review), we also aim to obtain satisfactory convergence properties upon refinement in the time discretization.

The aim of this paper is hence to propose a numerical scheme to solve PDECO problems with reaction–diffusion equations, which allows the commutation of the discretization and optimization steps. The discretization scheme that results in a symmetric system matrix also allows us to exploit the robust and theoretically-predictable convergence properties of suitable iterative methods for symmetric systems. In particular, we wish to derive a fast and efficient iterative numerical solver for the fully discretized problem by exploiting the algebraic structure of the system matrix and designing a suitable preconditioner. We highlight that the framework we present is general and can be applied to a wide range of nonlinear evolution PDEs. It also opens up the possibility of applying such a strategy to other time-stepping methods for a range of PDEs, which one may have reason to apply in order to achieve desirable stability properties or make use of implicit–explicit (IMEX) methods1.

The paper is organized as follows. In Section 2, we formulate a parameter identification problem for a system of reaction–diffusion equations in the PDECO context, which is linearized using the Sequential Quadratic Programming method in Section 3. The linear–quadratic PDECO problem that emerges is then transformed into a fully discretized system in the following two sections: in Section 4, we apply a Störmer–Verlet time discretization method to the optimality conditions in the OTD approach, while in Section 5 we first approximate the Lagrangian with midpoint and trapezoidal methods and then derive the optimality system. The system is then written in matrix–vector form in Section 6. To be able to solve the system of linear algebraic equations fast, in Section 7 we derive a suitable preconditioner. In Section 8 we introduce the backward Euler method for our PDECO problem and in Section 9 we compare numerical results for the Störmer–Verlet and classical backward Euler methods. The concluding remarks are presented in Section 10.

Notation

The matrices are denoted by uppercase letters in bold font. As the work in this paper involves block diagonal, sub- and super-diagonal matrices, we introduce a concise notation for their representation to enhance brevity in mathematical expressions. We denote block diagonal matrices, with N diagonal blocks that contain matrices \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{P}^i$$\end{document} or \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{R}^i$$\end{document} , by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\text {blkdiag}} \left( \left\{ \textbf{P}^i\right\} _{i=1}^N \right) \hspace{3mm} {\text {or}} \hspace{3mm} {\text {blkdiag}} \left( \left\{ \textbf{P}^i\right\} _{i=1}^{N_1}, \left\{ \textbf{R}^i\right\} _{i=1}^{N_2} \right) , \end{aligned}$$\end{document}

where the superscripts indicate the time-step of the elements involved in the block matrix and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_1+N_2=N$$\end{document} . Similarly, we define the operators for block sub- and super-diagonal matrices, with sub-diagonals \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{U}^i$$\end{document} and super-diagonals \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{V}^i$$\end{document} , as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\text {blksubdiag}} \left( \left\{ \textbf{U}^i\right\} _{i=1}^{N-1} \right) \hspace{3mm} {\text {and}} \hspace{3mm} {\text {blksupdiag}} \left( \left\{ \textbf{V}^i\right\} _{i=1}^{N-1} \right) , \end{aligned}$$\end{document}

respectively. Calligraphic typeface denotes the main system matrices and their preconditioners.

The vectors are represented by lowercase letters in bold font. To differentiate between various types of vectors utilized in different contexts, continuous variable vectors and vector functions are represented using bold Roman typeface, e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {u}} = [u,v]^\top $$\end{document} , whereas vectors containing fully discretized variables are indicated by bold italic, e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}^i$$\end{document} represents the variable u discretized over the mesh at the i-th time-step. Underlined bold italic vectors combine two variables discretized over the mesh and across some time interval, e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\underline{\boldsymbol{u}}:= \left[ (\boldsymbol{u}^1)^{\top },..., (\boldsymbol{u}^{N_t})^{\top }, (\boldsymbol{v}^1)^{\top },..., (\boldsymbol{v}^{N_t})^{\top }\right] ^{\top }$$\end{document} for time-steps \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1,..., N_t$$\end{document} .

Parameter Identification for Reaction–Diffusion Equations

To determine or approximately identify missing parameters or functions in a system of PDEs, we may formulate a constrained optimization problem with the PDEs as constraints. For the problems considered here, we think of the missing parameters as control variables, and aim to find optimal controls which minimize the misfit between the solutions of the PDEs (state variables) and observed data in the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} -norm, which is a frequent choice of norm for these problems. To ensure the existence of a minimum in the optimal control problem, we add to the cost functional defining the optimization problem a regularization term, minimizing the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2$$\end{document} -norm of controls. We are interested in finding solutions in a bounded Lipschitz domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega \subset \mathbb {R}^2$$\end{document} throughout the time interval (0, T) and we solve the optimization problem in the time–space cylinder \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Q:=(0,T)\times \varOmega $$\end{document} . The boundary of the spatial domain is denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\partial \varOmega $$\end{document} and the time–space boundary is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varGamma := (0,T)\times \partial \varOmega $$\end{document} .

For a PDECO problem with two state variables u, v and controls a, b, we consider the cost functional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {J} = \mathcal {J}(u,v,a,b)$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \mathcal {J} :={}&\frac{\alpha _1}{2} \Vert u-\hat{u}\Vert ^2_{L^2(Q)} +\frac{\alpha _2}{2} \Vert v-\hat{v}\Vert ^2_{L^2(Q)} + \frac{\beta _1}{2} \Vert a\Vert ^2_{L^2(Q)} + \frac{\beta _2}{2} \Vert b\Vert ^2_{L^2{(Q)}}, \end{aligned} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{u}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{v}$$\end{document} represent the observations, known as desired states. For the weights \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _i$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _i$$\end{document} ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i=1,2$$\end{document} ) we typically choose \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _i>\beta _i$$\end{document} to prioritize minimizing the misfit. In the case of large  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _i$$\end{document} , one would be penalizing the cost of implementing the controls, which may result in the state variables failing to match the desired states closely.

We formulate the PDECO problem for reaction–diffusion systems

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \min _{u,v,a,b} \mathcal {J}(u,v,a,b) \end{aligned}$$\end{document}

subject to:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} u_t - D_u \nabla ^2u + \varPhi (u,v)&= \gamma a \hspace{5mm} {\text { in }} Q, \\ v_t - D_v \nabla ^2v + \varPsi (u,v)&= \gamma b \hspace{5.3mm} {\text { in }} Q, \end{aligned} \end{aligned}$$\end{document}

with initial and boundary conditions

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} u(0,{\textbf {x}}) = u_0({\textbf {x}}), \hspace{2mm} v(0,{\textbf {x}}) = v_0({\textbf {x}})\hspace{5mm}&{\text { in }} \varOmega , \\ \frac{\partial u}{\partial \boldsymbol{\nu }} = \frac{\partial v}{\partial \boldsymbol{\nu }} = 0 \hspace{11.4mm}&{\text { on }} \varGamma , \end{aligned} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varPhi (u,v)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varPsi (u,v)$$\end{document} are nonlinear reaction terms and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} is a positive constant. The parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_u$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_v$$\end{document} denote the positive diffusion coefficients and, for the occurrence of Turing patterns, we require \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_u<D_v$$\end{document} (see e.g., [26, Sec. 2.3]). The control variables a and b represent the production rates for u and v, respectively. In this paper, we only consider the case where the control variables are acting on the whole domain, i.e., the controls are distributed, and where they act as source terms in the model equations. The initial conditions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0$$\end{document} are assumed to be known. While we consider zero-flux boundary conditions with the vector normal to the boundary denoted by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\nu }$$\end{document} , it is also possible to use Dirichlet boundary conditions instead.

The numerical scheme that we present in this paper is built upon the Schnakenberg model, a system of reaction–diffusion equations commonly used to model pattern formation, for which \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varPhi (u,v)= \gamma (u-u^2v)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varPsi (u,v)= \gamma u^2v$$\end{document} . Here, the given positive parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} can be interpreted as the relative strength of the reaction (and also control) terms. The equations were initially proposed by Schnakenberg to study limit cycles arising in chemical reactions [27]. They have since also been used to study self-organization phenomena in various biological systems, for example, the emergence of skin patterns [28] or initiation of plant root hair [29].

Since the reaction terms are nonlinear, the optimization problem (1)–(3) is also nonlinear. To resolve this, we need to utilize some nonlinear programming algorithm. In the following section, we derive the first-order optimality conditions using the Lagrange method and apply the Sequential Quadratic Programming (SQP) method; see [10, Sec. 5.9.2] and [30, Sec. 3.2].

Nonlinear Programming

Below, we perform linearization for the nonlinear problem (1)–(3) with the SQP method, making use of a Newton linearization for the PDE constraints. We commence by following the OTD approach, and via our specific choice of time-stepping and approximation of certain terms, we find that the same optimality conditions are obtained from both OTD and DTO formulations, thus achieving that the discretization and optimization operations are commutative.

Lagrange method

Let us start by formulating the continuous Lagrangian \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {L} = \mathcal {L}(u,v,a,b,p,q)$$\end{document} . We multiply the PDE constraints by the Lagrange multipliers and add them to the cost functional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {J}$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \mathcal {L} (u,v,a,b,p,q):= \mathcal {J} + \int _Q p \left( u_t - D_u \nabla ^2u + \varPhi (u,v) - \gamma a \right) dQ \\ + \int _Q q \left( v_t - D_v \nabla ^2v + \varPsi (u,v) - \gamma b \right) dQ, \end{aligned} \end{aligned}$$\end{document}

where p and q denote the multipliers, also called the adjoint variables. We omit including the boundary conditions within the Lagrangian (4) and the following sections, for the sake of brevity of presentation. The following optimality conditions are a result of finding the stationary points of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {L}$$\end{document} , which we obtain by taking Fréchet derivatives; for a comprehensive description of this process for a range of problems we refer to [10].

Firstly, from the partial derivatives with respect to the state variables u and v, we derive the adjoint equations:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} -p_t - D_u \nabla ^2 p + \varPhi _u(u,v) p + \varPsi _u(u,v) q + \alpha _1(u-\hat{u})&= 0 \hspace{5mm} {\text {in }} Q, \\ -q_t - D_v \nabla ^2 q + \varPhi _v(u,v) p + \varPsi _v(u,v) q + \alpha _2(v-\hat{v})&= 0 \hspace{5mm} {\text {in }} Q, \end{aligned} \end{aligned}$$\end{document}

equipped with final-time conditions and zero-flux boundary conditions:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} p(T,{\textbf {x}}) = q(T,{\textbf {x}})&= 0 \hspace{5mm} {\text {in }} \varOmega , \\ \frac{\partial p}{\partial \boldsymbol{\nu }} = \frac{\partial q}{\partial \boldsymbol{\nu }}&= 0 \hspace{5mm} {\text {on }} \varGamma . \end{aligned} \end{aligned}$$\end{document}

Next, the derivatives with respect to the controls give us the gradient or control equations:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \beta _1 a - \gamma p&= 0 \hspace{5mm} {\text {in }} Q,\\ \beta _2 b - \gamma q&= 0 \hspace{5mm} {\text {in }} Q. \end{aligned} \end{aligned}$$\end{document}

We note that the gradient equations are algebraic for this problem, and thus can be used to reduce the complexity of the optimization problem in later sections by substituting for the variables a and b. Finally, the set of optimality conditions is completed by studying optimality with respect to the Lagrange multipliers: the Fréchet derivatives will yield the original state equations and initial and boundary conditions as in (2)–(3).

Lagrange–SQP method

After obtaining the optimality conditions, we proceed to their linearization. For our optimality system, assuming previously calculated Newton iterates at the k-th iteration, marked in the subscripts of the variables, we need to solve for the variables u, v, p, q at the current \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(k+1)$$\end{document} -th Newton iteration. Using the more concise notation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\varPhi }:= \varPhi (u_k,v_k)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\varPsi }:= \varPsi (u_k,v_k)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\varPhi }_{u}:= \partial _{u}\varPhi (u_k,v_k)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{\varPsi }_{uv}:= \partial ^2_{uv} \varPsi (u_k,v_k)$$\end{document} , and similarly for other derivatives (where the subscripts of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varPhi $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varPsi $$\end{document} denote partial derivatives), we thus need to solve the linearized state equations:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} u_t - D_u \nabla ^2 u + {\bar{\varPhi }} + {\bar{\varPhi }}_u(u-u_k) + {\bar{\varPhi }}_v(v-v_k) - \frac{\gamma ^2}{\beta _1} p&= 0, \\ v_t - D_v\nabla ^2v + {\bar{\varPsi }} + {\bar{\varPsi }}_u(u-u_k) + {\bar{\varPsi }}_v(v-v_k)- \frac{\gamma ^2}{\beta _2} q&= 0, \\ \end{aligned} \end{aligned}$$\end{document}

and, upon an application of the SQP method, the linearized adjoint equations of the form:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&-p_t- D_u\nabla ^2 p + {\bar{\varPhi }}_u p + {\bar{\varPsi }}_u q + (\alpha _1 + {\bar{\varPhi }}_{uu} p_k + {\bar{\varPsi }}_{uu} q_k) u + ({\bar{\varPhi }}_{uv} p_k + {\bar{\varPsi }}_{uv} q_k) v \\&\hspace{5mm} = \alpha _1 \hat{u} + ({\bar{\varPhi }}_{uu} p_k + {\bar{\varPsi }}_{uu} q_k)u_k + ({\bar{\varPhi }}_{uv} p_k + {\bar{\varPsi }}_{uv} q_k) v_k, \\&- q_t - D_v \nabla ^2 q + {\bar{\varPhi }}_v p + {\bar{\varPsi }}_v q + ({\bar{\varPhi }}_{vu} p_k + {\bar{\varPsi }}_{vu} q_k) u + (\alpha _2 + {\bar{\varPhi }}_{vv} p_k + {\bar{\varPsi }}_{vv} q_k) v \\&\hspace{5mm} = \alpha _2 \hat{v} + ({\bar{\varPhi }}_{vu} p_k + {\bar{\varPsi }}_{vu} q_k) u_k + ({\bar{\varPhi }}_{vv} p_k + {\bar{\varPsi }}_{vv} q_k) v_k. \end{aligned} \end{aligned}$$\end{document}

The linearized equations (7)–(8) are equivalent to the first-order optimality conditions of a linear–quadratic optimization problem, obtained via an iteration of the SQP method, as described in [10]. The SQP cost functional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\mathcal {J}}$$\end{document} is related to the cost functional and Lagrangian of the original problem (1)–(3) as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \tilde{\mathcal {J}}(u,v,a,b):= \mathcal {J}^\prime (u_k,v_k,a_k,b_k){\textbf {w}} + \frac{1}{2} {\textbf {w}} ^\top \mathcal {L}^{\prime \prime }(u_k,v_k,a_k,b_k,p_k,q_k){\textbf {w}}, \end{aligned} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {w}}=[u-u_k, v-v_k, a-a_k,b-b_k]^\top $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {L}^{\prime \prime }$$\end{document} is the Hessian of (4) with respect to u, v, a, and b. At each SQP iteration, we therefore solve the quadratic optimization problem

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \min _{u,v,a,b} \mathcal {\tilde{J}}(u,v,a,b) \end{aligned}$$\end{document}

subject to the system of linear reaction–diffusion equations:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} u_t - D_u \nabla ^2 u + {\bar{\varPhi }} + {\bar{\varPhi }}_u(u-u_k) + {\bar{\varPhi }}_v(v-v_k) - \gamma a&= 0 \hspace{11.4mm} {\text { in }} Q, \\ v_t - D_v\nabla ^2v + {\bar{\varPsi }} + {\bar{\varPsi }}_u(u-u_k) + {\bar{\varPsi }}_v(v-v_k)- \gamma b&= 0 \hspace{11.4mm} {\text { in }} Q, \\ u(0,{\textbf {x}}) = u_0({\textbf {x}}), \hspace{2mm} v(0,{\textbf {x}})&= v_0({\textbf {x}})\hspace{6mm} {\text { in }} \varOmega , \\ \frac{\partial u}{\partial \boldsymbol{\nu }} = \frac{\partial v}{\partial \boldsymbol{\nu }}&= 0 \hspace{11.4mm} {\text { on }} \varGamma . \end{aligned} \end{aligned}$$\end{document}

We now evaluate all the terms in (9). The first term is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \mathcal {J}^\prime (u_k,v_k,a_k,b_k){\textbf {w}} ={}&\alpha _1 \int _Q (u_k-\hat{u})(u-u_k) \, dQ + \alpha _2 \int _Q (v_k-\hat{v})(v-v_k) \, dQ \\&+ \beta _1 \int _Q a_k(a-a_k) \, dQ + \beta _2 \int _Q b_k(b-b_k) \, dQ, \end{aligned} \end{aligned}$$\end{document}

and the second term is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \frac{1}{2} {\textbf {w}}^\top \mathcal {L}^{\prime \prime }(u_k,v_k,a_k,b_k{,p_k,q_k}){\textbf {w}} = \frac{1}{2} \int _Q {\textbf {w}}^\top \textbf{D}_{\alpha , \beta } {\textbf {w}} \, dQ + \frac{1}{2} \int _Q {\textbf {w}}^\top \textbf{A}_{\varPhi ,\varPsi } {\textbf {w}} \, dQ, \end{aligned} \end{aligned}$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \textbf{D}_{\alpha , \beta } = \begin{bmatrix} \alpha _1 {\text {id}} & 0 & 0 & 0 \\ 0 & \alpha _2 {\text {id}} & 0 & 0 \\ 0 & 0 & \beta _1 {\text {id}} & 0 \\ 0 & 0 & 0 & \beta _2 {\text {id}} \end{bmatrix},\quad ~ \textbf{A}_{\varPhi ,\varPsi } = \begin{bmatrix} p_k {\bar{\varPhi }}_{uu} + q_k {\bar{\varPsi }}_{uu} & p_k {\bar{\varPhi }}_{uv} + q_k {\bar{\varPsi }}_{uv} & 0 & 0 \\ p_k {\bar{\varPhi }}_{vu} + q_k {\bar{\varPsi }}_{vu} & p_k {\bar{\varPhi }}_{vv}+ q_k {\bar{\varPsi }}_{vv}& 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}, $$\end{document}

and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {id}}$$\end{document} denotes the identity operator.

In the remainder of the work, we consider reaction–diffusion equations specific to the Schnakenberg model, but the approach considered here can be applied to any system in the form (2). The cost functional (9) can be split into two terms:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathcal {\tilde{J}}(u,v,a,b) = \int _Q (\mathcal {\tilde{J}}_1 + \mathcal {\tilde{J}}_2) \, dQ, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {\tilde{J}}_1$$\end{document} contains the terms from (12) and the first integral in (13):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \mathcal {\tilde{J}}_1 ={}&\alpha _1 (u_k-\hat{u})(u-u_k) + \alpha _2 (v_k-\hat{v})(v-v_k) + \beta _1 a_k(a-a_k) + \beta _2 b_k(b-b_k) \\&+ \frac{\alpha _1}{2}(u-u_k)^2 + \frac{\alpha _2}{2}(v-v_k)^2 + \frac{\beta _1}{2}(a-a_k)^2 + \frac{\beta _2}{2}(b-b_k)^2 \end{aligned} \end{aligned}$$\end{document}

and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {\tilde{J}}_2$$\end{document} denotes the second term in (13):

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathcal {\tilde{J}}_2 = \frac{1}{2} {\textbf {w}}^{\top } \textbf{A}_{\varPhi ,\varPsi } {\textbf {w}}, \quad {\text { where }} \; \; \textbf{A}_{\varPhi ,\varPsi } = \begin{bmatrix} 2\gamma (q_k-p_k) v_k & 2\gamma (q_k-p_k) u_k & 0 & 0 \\ 2\gamma (q_k-p_k) u_k & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{bmatrix}. \end{aligned}$$\end{document}

As \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {\tilde{J}}$$\end{document} is a cost functional, and it will later be part of a Lagrangian corresponding to the problem (10)–(11), the terms independent of u, v, a, b within \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {\tilde{J}}_1$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {\tilde{J}}_2$$\end{document} will not appear in the optimality conditions. We can therefore conveniently add the constant terms \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\alpha _1}{2}\hat{u}^2$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\alpha _2}{2}\hat{v}^2$$\end{document} to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {\tilde{J}}_1$$\end{document} which allows us to reconstruct the original cost functional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {\tilde{J}}$$\end{document} , while the other terms in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {\tilde{J}}_1$$\end{document} will be constant. Specifically,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathcal {\tilde{J}}_1 = \mathcal {J} - \frac{\alpha _1}{2}\hat{u}^2 - \frac{\alpha _2}{2}\hat{v}^2 - \frac{\alpha _1}{2} u_k^2 + \alpha _1 \hat{u}u_k - \frac{\alpha _2}{2} v_k^2 + \alpha _2 \hat{v}v_k - \frac{\beta _1}{2}a_k^2 - \frac{\beta _2}{2}b_k^2. \end{aligned}$$\end{document}

Omitting the constant terms, the SQP method consists of minimizing at each iteration the cost functional

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \mathcal {\tilde{J}} ={}&\frac{\alpha _1}{2} \int _Q (u-\hat{u})^2 \, dQ + \frac{\alpha _2}{2} \int _Q (v-\hat{v})^2 \, dQ + \frac{\beta _1}{2} \int _Q a^2 \, dQ + \frac{\beta _2}{2} \int _Q b^2 \, dQ \\&+\int _Q \gamma (q_k-p_k)v_k(u-u_k)^2 \, dQ + \int _Q 2\gamma (q_k-p_k)u_k(u-u_k)(v-v_k) \, dQ, \end{aligned} \end{aligned}$$\end{document}

and we can rearrange the linearized state equations to obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} u_t - D_u \nabla ^2u + \gamma (1-2u_kv_k)u - \gamma u_k^2 v - \gamma a + 2\gamma u_k^2v_k&= 0 \hspace{5mm} {\text {in }} Q, \\ v_t - D_v \nabla ^2v + 2\gamma u_kv_k u + \gamma u_k^2v - \gamma b - 2\gamma u_k^2v_k&= 0 \hspace{5mm} {\text {in }} Q, \end{aligned} \end{aligned}$$\end{document}

with the initial and boundary conditions (3). We therefore reformulated the nonlinear optimization problem (1)–(3) as the quadratic optimization problem (14), (15), (3). The next step is to provide a method for which the discretization and optimization steps commute for the SQP problem.

Optimize-then-Discretize Approach

To solve the SQP problem (14), (15), (3) using the Optimize-then-Discretize approach, we apply the Lagrange multiplier technique described in Section 3.2 to derive the optimality conditions. With the Schnakenberg equations as constraints, the first two optimality conditions recover the linearized state equations (15), and the other two are the adjoint equations:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} -p_t - D_u\nabla ^2 p + \alpha _1 u + 2\gamma v_k \left( q_k-p_k\right) u + 2\gamma u_k \left( q_k-p_k\right) v + \gamma (1- 2 u_k v_k) p&\\ {}+2\gamma u_k v_k q = \alpha _1 \hat{u} + 4\gamma u_k v_k\left( q_k - p_k \right) {}&, \\ - q_t - D_v \nabla ^2 q + 2\gamma u_k \left( q_k-p_k\right) u + \alpha _2 v + \gamma u_k^2 \left( q - p\right) = \alpha _2 \hat{v} + 2\gamma u_k^2 \left( q_k - p_k\right) {}&, \end{aligned} \end{aligned}$$\end{document}

the same as in (8). We also recover the same linear gradient equations as in (6), and initial and boundary conditions in (3) and (5), respectively.

Störmer–Verlet scheme

Störmer–Verlet type schemes for time-dependent optimization problems with linear PDEs as constraints, in particular the heat equation, were considered in [22, 23]; see also [31, 32]. We now introduce such a scheme in the context of the nonlinear problems of this paper, denoting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {u}} = [u,v]^\top $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {p}} = [p,q]^\top $$\end{document} and considering the state and adjoint equations in the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\textbf {u}}_t={\textbf {f}}({\textbf {u}},{\textbf {p}}), \hspace{5mm} {\textbf {p}}_t={\textbf {g}}({\textbf {u}},{\textbf {p}}). \end{aligned}$$\end{document}

Notice that using the gradient equations (6), the right-hand side functions in both state and adjoint equations can be written as functions of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {u}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {p}}$$\end{document} . We therefore have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\textbf {f}}({\textbf {u}},{\textbf {p}}) = \begin{bmatrix} f_1({\textbf {u}},{\textbf {p}})\\ f_2({\textbf {u}},{\textbf {p}}) \end{bmatrix} \quad {\text { and }} \quad {\textbf {g}}({\textbf {u}},{\textbf {p}}) = \begin{bmatrix} g_1({\textbf {u}},{\textbf {p}})\\ g_2({\textbf {u}},{\textbf {p}}) \end{bmatrix}, $$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} f_1({\textbf {u}},{\textbf {p}}) ={}&D_u \nabla ^2u - \gamma (1-2u_k v_k)u + \gamma u_k^2 v + \frac{\gamma ^2}{\beta _1} p - 2\gamma u_k^2v_k, \\ f_2({\textbf {u}},{\textbf {p}}) ={}&D_v \nabla ^2v - 2\gamma u_k v_k u - \gamma u_k^2v + \frac{\gamma ^2}{\beta _2} q + 2\gamma u_k^2 v_k, \\ g_1({\textbf {u}},{\textbf {p}}) ={}&-D_u \nabla ^2p + \gamma (1-2u _k v_k)p + 2\gamma u_k v_k q + 2\gamma v_k(q_k-p_k)u \\&+ 2\gamma u_k(q_k-p_k)v - \alpha _1(\hat{u}-u) - 4\gamma u_kv_k(q_k-p_k),\\ g_2({\textbf {u}},{\textbf {p}}) ={}&- D_v \nabla ^2q - \gamma u_k^2 p + \gamma u_k^2 q + 2\gamma u_k(q_k-p_k) u - \alpha _2(\hat{v}-v) \\&- 2\gamma u_k^2(q_k-p_k). \end{aligned} \end{aligned}$$\end{document}

Similarly to [25, Eq. (1.24)], we discretize equations (17) using the time-stepping scheme described as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} {\textbf {p}}^{n+\frac{1}{2}}&= {\textbf {p}}^n + \frac{\tau }{2} {\textbf {g}}({\textbf {u}}^n,{\textbf {p}}^{n+\frac{1}{2}}), \\ {\textbf {u}}^{n+1}&= {\textbf {u}}^n + \frac{\tau }{2}\left( {\textbf {f}}({\textbf {u}}^n,{\textbf {p}}^{n+\frac{1}{2}}) + {\textbf {f}}({\textbf {u}}^{n+1},{\textbf {p}}^{n+\frac{1}{2}})\right) , \\ {\textbf {p}}^{n+1}&= {\textbf {p}}^{n+\frac{1}{2}} + \frac{\tau }{2} {\textbf {g}}({\textbf {u}}^{n+1},{\textbf {p}}^{n+\frac{1}{2}}), \end{aligned} \end{aligned}$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=0,1,...,N_t-1$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_t$$\end{document} is the number of time-steps of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} , and with superscripts denoting the time indices. To maintain time stepping of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} in each variable, we reverse the last equation in (18) by one step:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} {\textbf {p}}^n = {\textbf {p}}^{n-\frac{1}{2}} + \frac{\tau }{2} {\textbf {g}}({\textbf {u}}^n,{\textbf {p}}^{n-\frac{1}{2}}), \end{aligned}$$\end{document}

Substituting (19) into the first equation in (18), we obtain the scheme

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} {\textbf {p}}^{n+\frac{1}{2}}&= {\textbf {p}}^{n-\frac{1}{2}} + \frac{\tau }{2}\left( {\textbf {g}}({\textbf {u}}^n,{\textbf {p}}^{n-\frac{1}{2}}) + {\textbf {g}}({\textbf {u}}^n,{\textbf {p}}^{n+\frac{1}{2}})\right) ,&{\text { for }} \; n=1,2,...,N_t-1,\\ {\textbf {u}}^{n+1}&= {\textbf {u}}^n + \frac{\tau }{2}\left( {\textbf {f}}({\textbf {u}}^n,{\textbf {p}}^{n+\frac{1}{2}}) + {\textbf {f}}({\textbf {u}}^{n+1},{\textbf {p}}^{n+\frac{1}{2}})\right) ,&{\text { for }} \; n=0,1,...,N_t-1, \\ \end{aligned} \end{aligned}$$\end{document}

with the values at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=N_t - \frac{1}{2}$$\end{document} in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {p}}$$\end{document} calculated using

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} {\textbf {p}}^{\frac{1}{2}}&= {\textbf {p}}^0 + \frac{\tau }{2} {\textbf {g}}({\textbf {u}}^0,{\textbf {p}}^{\frac{1}{2}}), \\ {\textbf {p}}^{N_t}&= {\textbf {p}}^{N_t-\frac{1}{2}} + \frac{\tau }{2} {\textbf {g}}({\textbf {u}}^{N_t},{\textbf {p}}^{N_t-\frac{1}{2}}). \end{aligned} \end{aligned}$$\end{document}

We note that the expressions for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {u}}^{n+1}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {p}}^{n+\frac{1}{2}}$$\end{document} in (20) resemble trapezoidal rule type schemes, and this fact is used in Section 5 when constructing the scheme using the DTO approach. Due to the final-time condition for the adjoint variables, the values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {p}}^{N_t}$$\end{document} in (21) will be equal to zero. The same time-stepping will be followed when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {g}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {f}}$$\end{document} also involve the previous iterates of the variables, e.g., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {g}}({\textbf {u}}^n,{\textbf {p}}^{n-\frac{1}{2}}) = {\textbf {g}}({\textbf {u}}^n,{\textbf {u}}^n_k,{\textbf {p}}^{n-\frac{1}{2}}, {\textbf {p}}^{n-\frac{1}{2}}_k)$$\end{document} .

Applying the scheme over the time interval (0, T) yields the semi-discretized state equations

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&-u^n - \frac{\tau D_u}{2} \nabla ^2 u^n + \frac{\tau \gamma }{2} (1-2u_k^n v_k^n) u^n + u^{n+1} - \frac{\tau D_u}{2} \nabla ^2 u^{n+1} \\&{}+ \frac{\tau \gamma }{2}( 1-2u_k^{n+1} v_k^{n+1}) u^{n+1} - \frac{\tau \gamma }{2} ( u_k^n) ^2 v^n - \frac{\tau \gamma }{2} ( u_k^{n+1}) ^2 v^{n+1} \\&{}-\tau \gamma a^{n+\frac{1}{2}} = -\tau \gamma \left[ \left( u_k^n\right) ^2 v_k^n + ( u_k^{n+1}) ^2 v_k^{n+1} \right] \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\tau \gamma u_k^n v_k^n u^n + \tau \gamma u_k^{n+1} v_k^{n+1} u^{n+1} -v^n - \frac{\tau D_v }{2}\nabla ^2 v^n + \frac{\tau \gamma }{2} \left( u_k^n\right) ^2 v^n + v^{n+1} \\&{}- \frac{\tau D_v}{2} \nabla ^2 v^{n+1} + \frac{\tau \gamma }{2} ( u_k^{n+1}) ^2 v^{n+1} - \tau \gamma b^{n+\frac{1}{2}} = \tau \gamma \left[ \left( u_k^n\right) ^2 v_k^n + (u_k^{n+1})^2 v_k^{n+1} \right] , \end{aligned} \end{aligned}$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=0,1,...,N_t-1$$\end{document} , with the initial conditions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^0 = u(0,{\textbf {x}})$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^0 = v(0,{\textbf {x}})$$\end{document} in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega $$\end{document} . For the semi-discretized adjoint equations, starting with the first half step, we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\frac{\tau \alpha _1}{2} u^0 + \tau \gamma v_k^0 (q_k^{\frac{1}{2}} - p_k^{\frac{1}{2}}) u^0 + \tau \gamma u_k^0 (q_k^{\frac{1}{2}} - p_k^{\frac{1}{2}}) v^0 + p^0 - p^{\frac{1}{2}} - \frac{\tau D_u}{2}\nabla ^2 p^{\frac{1}{2}} \\&{} + \frac{\tau \gamma }{2} (1- 2 u_k^0 v_k^0) p^{\frac{1}{2}} + \tau \gamma u_k^0 v_k^0 q^{\frac{1}{2}} = \frac{\tau \alpha _1}{2} \hat{u}^0 + 2 \tau \gamma u_k^0 v_k^0 (q_k^{\frac{1}{2}} - p_k^{\frac{1}{2}}) \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\tau \gamma u_k^0 (q_k^{\frac{1}{2}} - p_k^{\frac{1}{2}}) u^0 + \frac{\tau \alpha _2}{2} v^0 -\frac{\tau \gamma }{2} ( u_k^0 )^2 p^{\frac{1}{2}} + q^0 - q^{\frac{1}{2}} - \frac{\tau D_v}{2}\nabla ^2 q^{\frac{1}{2}} \\&{}+\frac{\tau \gamma }{2} \left( u_k^0 \right) ^2 q^{\frac{1}{2}} = \frac{\tau \alpha _2}{2} \hat{v}^0 + \tau \gamma \left( u_k^0 \right) ^2 (q_k^{\frac{1}{2}} - p_k^{\frac{1}{2}}). \end{aligned} \end{aligned}$$\end{document}

The adjoint equations for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=1,...,N_t-1$$\end{document} are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\tau \alpha _1 u^n + \tau \gamma v_k^n \left[ (q_k^{n-\frac{1}{2}} - p_k^{n-\frac{1}{2}}) + (q_k^{n+\frac{1}{2}} - p_k^{n+\frac{1}{2}}) \right] u^n + \tau \gamma u_k^n \left[ (q_k^{n-\frac{1}{2}} - p_k^{n-\frac{1}{2}}) \right. \\&\left. {}\!\!+ (q_k^{n+\frac{1}{2}} - p_k^{n+\frac{1}{2}}) \right] v^n + p^{n-\frac{1}{2}} - \frac{\tau D_u}{2}\nabla ^2 p^{n-\frac{1}{2}} + \frac{\tau \gamma }{2} (1- 2 u_k^n v_k^n) p^{n-\frac{1}{2}} \\&{}- p^{n+\frac{1}{2}} - \frac{\tau D_u}{2}\nabla ^2 p^{n+\frac{1}{2}} + \frac{\tau \gamma }{2} (1- 2 u_k^n v_k^n) p^{n+\frac{1}{2}} + \tau \gamma u_k^n v_k^n q^{n-\frac{1}{2}} \\&{}+ \tau \gamma u_k^n v_k^n q^{n+\frac{1}{2}} = \tau \alpha _1 \hat{u}^n + 2 \tau \gamma u_k^n v_k^n \left[ (q_k^{n-\frac{1}{2}} - p_k^{n-\frac{1}{2}}) + (q_k^{n+\frac{1}{2}} - p_k^{n+\frac{1}{2}}) \right] \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\tau \gamma u_k^n \left[ (q_k^{n-\frac{1}{2}} - p_k^{n-\frac{1}{2}}) + (q_k^{n+\frac{1}{2}} - p_k^{n+\frac{1}{2}}) \right] u^n + \tau \alpha _2 v^n - \frac{\tau \gamma }{2} ( u_k^n )^2 p^{n-\frac{1}{2}} \\&{} -\frac{\tau \gamma }{2} ( u_k^n )^2 p^{n+\frac{1}{2}} + q^{n-\frac{1}{2}} - \frac{\tau D_v}{2}\nabla ^2 q^{n-\frac{1}{2}} + \frac{\tau \gamma }{2} \left( u_k^n \right) ^2 q^{n-\frac{1}{2}} - q^{n+\frac{1}{2}} \\&{}- \frac{\tau D_v}{2}\nabla ^2 q^{n+\frac{1}{2}} +\frac{\tau \gamma }{2} \left( u_k^n \right) ^2 q^{n+\frac{1}{2}} = \tau \alpha _2 \hat{v}^n + \tau \gamma \left( u_k^n \right) ^2 \left[ (q_k^{n-\frac{1}{2}} - p_k^{n-\frac{1}{2}}) \right. \\&{}+\left. (q_k^{n+\frac{1}{2}} - p_k^{n+\frac{1}{2}}) \right] . \end{aligned} \end{aligned}$$\end{document}

At the final half-time-step, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\frac{\tau \alpha _1}{2} u^{N_t} + \tau \gamma v_k^{N_t} ( q_k^{N_t-\frac{1}{2}} - p_k^{N_t-\frac{1}{2}}) u^{N_t} + \tau \gamma u_k^{N_t} (q_k^{N_t-\frac{1}{2}} - p_k^{N_t-\frac{1}{2}}) v^{N_t} + p^{N_t-\frac{1}{2}} \\&{} - \frac{\tau D_u}{2}\nabla ^2 p^{N_t-\frac{1}{2}} + \frac{\tau \gamma }{2} (1- 2 u_k^{N_t} v_k^{N_t}) p^{N_t-\frac{1}{2}} + \tau \gamma u_k^{N_t} v_k^{N_t} q^{N_t-\frac{1}{2}} = \frac{\tau \alpha _1}{2} \hat{u}^{N_t} \\&{}+ 2 \tau \gamma u_k^{N_t} v_k^{N_t} (q_k^{N_t-\frac{1}{2}} - p_k^{N_t-\frac{1}{2}}) \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\tau \gamma u_k^{N_t} (q_k^{N_t-\frac{1}{2}} - p_k^{N_t-\frac{1}{2}}) u^{N_t} + \frac{\tau \alpha _2}{2} v^{N_t} - \frac{\tau \gamma }{2} ( u_k^{N_t} )^2 p^{N_t-\frac{1}{2}} + q^{N_t-\frac{1}{2}} \\&{}- \frac{\tau D_v}{2}\nabla ^2 q^{N_t-\frac{1}{2}} + \frac{\tau \gamma }{2} (u_k^{N_t} )^2 q^{N_t-\frac{1}{2}} = \frac{\tau \alpha _2}{2} \hat{v}^{N_t} + \tau \gamma (u_k^{N_t} )^2 (q_k^{N_t-\frac{1}{2}} - p_k^{N_t-\frac{1}{2}}). \end{aligned} \end{aligned}$$\end{document}

Finite element discretization

For the spatial discretization of (22)–(29) we apply a finite element method and use P1 elements and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_x$$\end{document} nodes in the discretization of the domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega $$\end{document} . Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{M}$$\end{document} denote the mass matrix and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{K}$$\end{document} the stiffness matrix, with entries

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left[ \textbf{M}\right] _{r,s} = \int _\varOmega \varphi _r \varphi _s \, dx \; \; {\text { and }} \; \; \left[ \textbf{K}\right] _{r,s} = \int _\varOmega \nabla \varphi _r \cdot \nabla \varphi _s \, dx, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _r$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varphi _s$$\end{document} denote the finite element basis functions. Both these matrices are symmetric, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{M}$$\end{document} is positive definite, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{K}$$\end{document} is positive semi-definite (positive definite when considering Dirichlet boundary conditions). For the terms arising from linearization, we define the pseudo-mass matrices \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{M}_{u^i v^i}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{M}_{u^i(q^j - p^j)}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{M}_{(u^i)^2}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \textbf{M}_{v^i(q^j - p^j)}$$\end{document} as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \left[ \textbf{M}_{u^i v^i}\right] _{r,s} := \int _\varOmega u_k^i v_k^i \varphi _r \varphi _s \, dx, \hspace{5mm}&\left[ \textbf{M}_{u^i(q^j - p^j)}\right] _{r,s} := \int _\varOmega u_k^i(q_k^j - p_k^j) \varphi _r \varphi _s \, dx, \\ \left[ \textbf{M}_{(u^i)^2}\right] _{r,s} := \int _\varOmega (u_k^i)^2 \varphi _r \varphi _s \, dx, \hspace{5mm}&\left[ \textbf{M}_{v^i(q^j - p^j)}\right] _{r,s} := \int _\varOmega v_k^i(q_k^j - p_k^j) \varphi _r \varphi _s \, dx, \end{aligned} \end{aligned}$$\end{document}

which are also symmetric. The superscripts i and j of the variables involved denote the i-th and j-th time points ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^i=i \tau $$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^j=j \tau $$\end{document} ). For the right-hand sides of the discretized equations, we define the vectors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{d}^i$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{c}^{i,j}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{h}^{i,j}$$\end{document} , given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \left[ \boldsymbol{d}^i \right] _{r}&:= \int _\varOmega 2 \gamma (u_k^i)^2 v_k^i \varphi _r \, dx, \\ \left[ \boldsymbol{c}^{i,j} \right] _{r}&:= \int _\varOmega \left[ \alpha _1 \hat{u}^i + 4\gamma u_k^i v_k^i (q_k^j - p_k^j ) \right] \varphi _r \, dx, \\ \left[ \boldsymbol{h}^{i,j} \right] _{r}&:= \int _\varOmega \left[ \alpha _2 \hat{v}^i + 2\gamma ( u_k^i)^2 (q_k^j - p_k^j )\right] \varphi _r \, dx, \end{aligned} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{u}^i({\textbf {x}})=\hat{u}(t^i,{\textbf {x}})$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{v}^i({\textbf {x}})=\hat{v}(t^i,{\textbf {x}})$$\end{document} . We also define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\hat{u}}^i$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[\boldsymbol{\hat{u}}^i]_r = \int _\varOmega \hat{u}^i \psi _r \, dx$$\end{document} , similarly for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\hat{v}}^i$$\end{document} . Incorporating the above matrices and vectors yields the fully discretized state equations

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\left[ -\textbf{M} + \frac{\tau D_u}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M} - \tau \gamma \textbf{M}_{u^n v^n} \right] \boldsymbol{u}^n + \left[ \textbf{M} + \frac{\tau D_u}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M} \right. \\&{}-\left. \tau \gamma \textbf{M}_{u^{n+1} v^{n+1}} \right] \boldsymbol{u}^{n+1} - \frac{\tau \gamma }{2} \textbf{M}_{(u^n)^2 } \boldsymbol{v}^n - \frac{\tau \gamma }{2} \textbf{M}_{(u^{n+1})^2 } \boldsymbol{v}^{n+1} \\&{}- { \tau \gamma } \textbf{M} \boldsymbol{a}^{n+\frac{1}{2}} = - \frac{\tau }{2} \left( \boldsymbol{d}^n + \boldsymbol{d}^{n+1} \right) \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\tau \gamma \textbf{M}_{u^n v^n} \boldsymbol{u}^n + \tau \gamma \textbf{M}_{u^{n+1} v^{n+1}} \boldsymbol{u}^{n+1} + \left[ -\textbf{M} + \frac{\tau D_v}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M}_{(u^n)^2 } \right] \boldsymbol{v}^n \\&{}+ \left[ \textbf{M} + \frac{\tau D_v}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M}_{(u^{n+1})^2}\right] \boldsymbol{v}^{n+1} - \tau \gamma \textbf{M} \boldsymbol{b}^{n+\frac{1}{2}} = \frac{\tau }{2} \left( \boldsymbol{d}^n + \boldsymbol{d}^{n+1} \right) , \end{aligned} \end{aligned}$$\end{document}

with initial conditions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{M} \boldsymbol{u}^0 = \boldsymbol{u}_0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{M} \boldsymbol{v}^0 = \boldsymbol{v}_0$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[\boldsymbol{u}_0]_r = \int _\varOmega u_0 \psi _r \, dx$$\end{document} , similarly for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{v}_0$$\end{document} , using the initial conditions given in (3). For the first half step, the discretized adjoint equations are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\frac{\tau \alpha _1}{2} \textbf{M} \boldsymbol{u}^0 + \tau \gamma \textbf{M}_{v^0(q^{\frac{1}{2}} - p^{\frac{1}{2}})} \boldsymbol{u}^0 + \tau \gamma \textbf{M}_{u^0(q^{\frac{1}{2}} - p^{\frac{1}{2}})} \boldsymbol{v}^0 + \textbf{M} \boldsymbol{p}^0 \\&{}+ \left[ - \textbf{M} + \frac{\tau D_u}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M} - \tau \gamma \textbf{M}_{u^0v^0} \right] \boldsymbol{p}^{\frac{1}{2}} + \tau \gamma \textbf{M}_{u^0 v^0} \boldsymbol{q}^{\frac{1}{2}} = \frac{\tau }{2}\boldsymbol{c}^{0,\frac{1}{2}} \\ \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\tau \gamma \textbf{M}_{u^0(q^{\frac{1}{2}} - p^{\frac{1}{2}})} \boldsymbol{u}^0 + \frac{\tau \alpha _2}{2} \textbf{M} \boldsymbol{v}^0 - \frac{\tau \gamma }{2} \textbf{M}_{(u^0)^2} \boldsymbol{p}^{\frac{1}{2}} + \textbf{M} \boldsymbol{q}^0 \\&{} + \left[ -\textbf{M} + \frac{\tau D_v}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M}_{(u^0)^2} \right] \boldsymbol{q}^{\frac{1}{2}} = \frac{\tau }{2} \boldsymbol{h}^{0,\frac{1}{2}}. \end{aligned} \end{aligned}$$\end{document}

For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=1,...,N_t-1$$\end{document} , we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\tau \alpha _1 \textbf{M} \boldsymbol{u}^n + \tau \gamma \left[ \textbf{M}_{v^n(q^{n-\frac{1}{2}} - p^{n-\frac{1}{2}})} + \textbf{M}_{v^n(q^{n+\frac{1}{2}} - p^{n+\frac{1}{2}})} \right] \boldsymbol{u}^n \\&{}+ \tau \gamma \left[ \textbf{M}_{u^n(q^{n-\frac{1}{2}} - p^{n-\frac{1}{2}})} + \textbf{M}_{u^n(q^{n+\frac{1}{2}} - p^{n+\frac{1}{2}})} \right] \boldsymbol{v}^n + \left[ \textbf{M} + \frac{\tau D_u}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M} \right. \\&\left. {}\!\!- \tau \gamma \textbf{M}_{u^nv^n} \right] \boldsymbol{p}^{n-\frac{1}{2}} + \left[ - \textbf{M} + \frac{\tau D_u}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M} - \tau \gamma \textbf{M}_{u^nv^n} \right] \boldsymbol{p}^{n+\frac{1}{2}} \\&{}+ \tau \gamma \textbf{M}_{u^n v^n} \boldsymbol{q}^{n-\frac{1}{2}} + \tau \gamma \textbf{M}_{u^n v^n} \boldsymbol{q}^{n+\frac{1}{2}} = \frac{\tau }{2} \left( \boldsymbol{c}^{n, n-\frac{1}{2}} + \boldsymbol{c}^{n, n+\frac{1}{2}} \right) \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\tau \gamma \left[ \textbf{M}_{u^n(q^{n-\frac{1}{2}} - p^{n-\frac{1}{2}})} + \textbf{M}_{u^n(q^{n+\frac{1}{2}} - p^{n+\frac{1}{2}})} \right] \boldsymbol{u}^n +\tau \alpha _2 \textbf{M} \boldsymbol{v}^n - \frac{\tau \gamma }{2} \textbf{M}_{(u^n)^2} \boldsymbol{p}^{n-\frac{1}{2}} \\&- \frac{\tau \gamma }{2} \textbf{M}_{(u^n)^2} \boldsymbol{p}^{n+\frac{1}{2}} + \left[ \textbf{M} + \frac{\tau D_v}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M}_{(u^n)^2} \right] \boldsymbol{q}^{n-\frac{1}{2}} \\&+ \left[ -\textbf{M} + \frac{\tau D_v}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M}_{(u^n)^2} \right] \boldsymbol{q}^{n+\frac{1}{2}} = \frac{\tau }{2} \left( \boldsymbol{h}^{n, n-\frac{1}{2}} + \boldsymbol{h}^{n, n+\frac{1}{2}} \right) , \end{aligned} \end{aligned}$$\end{document}

and at the final half step, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\frac{\tau \alpha _1}{2} \textbf{M} \boldsymbol{u}^{N_t} + \tau \gamma \textbf{M}_{v^{N_t}(q^{N_t-\frac{1}{2}} - p^{N_t-\frac{1}{2}})} \boldsymbol{u}^{N_t} + \tau \gamma \textbf{M}_{u^{N_t}(q^{N_t-\frac{1}{2}} - p^{N_t-\frac{1}{2}})} \boldsymbol{v}^{N_t} + \left[ \textbf{M} \right. \\&\left. {}+ \frac{\tau D_u}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M} - \tau \gamma \textbf{M}_{u^{N_t}v^{N_t}} \right] \boldsymbol{p}^{N_t-\frac{1}{2}} + \tau \gamma \textbf{M}_{u^{N_t} v^{N_t}} \boldsymbol{q}^{N_t-\frac{1}{2}} = \frac{\tau }{2}\boldsymbol{c}^{N_t,N_t-\frac{1}{2}} \\ \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\tau \gamma \textbf{M}_{u^{N_t}(q^{N_t-\frac{1}{2}} - p^{N_t-\frac{1}{2}})} \boldsymbol{u}^{N_t} + \frac{\tau \alpha _2}{2} \textbf{M} \boldsymbol{v}^{N_t} - \frac{\tau \gamma }{2} \textbf{M}_{(u^{N_t})^2} \boldsymbol{p}^{N_t-\frac{1}{2}} \\&{}+ \left[ \textbf{M} + \frac{\tau D_v}{2} \textbf{K} + \frac{\tau \gamma }{2} \textbf{M}_{(u^{N_t})^2} \right] \boldsymbol{q}^{N_t-\frac{1}{2}} = \frac{\tau }{2} \boldsymbol{h}^{N_t,N_t-\frac{1}{2}}. \end{aligned} \end{aligned}$$\end{document}

Finally, according to the gradient equations (6), the controls depend on the adjoint variables which we evaluate at half time-steps:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \beta _1 \boldsymbol{a}^{n+\frac{1}{2}} - \gamma \boldsymbol{p}^{n+\frac{1}{2}} = \boldsymbol{0}, \hspace{5mm} \beta _2 \boldsymbol{b}^{n+\frac{1}{2}} - \gamma \boldsymbol{q}^{n+\frac{1}{2}} = \boldsymbol{0}, \quad {\text { for }} \; n=0,1,...,N_t-1. \end{aligned} \end{aligned}$$\end{document}

We have thus arrived at a fully discretized system following the Optimize-then-Discretize approach. To solve the system, we will form a system matrix for the all-at-once problem. First, however, we show that the Störmer–Verlet method leads to the same discrete optimality conditions when considering the Discretize-then-Optimize approach.

Discretize-then-Optimize Approach

In this section, we approximate the Lagrangian for the SQP problem by introducing a time-stepping strategy treating each term individually, to achieve a match with the time discretization in Section 4. The difficulty arises from the linearized terms that couple the state and adjoint variables, thus combining integer and half time-steps within the same term.

We note that the Störmer–Verlet scheme (20) is similar to a trapezoidal rule approximation, equivalent to the implicit midpoint rule for linear equations. We therefore use these quadrature rules, with the notation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^n \approx u(t^n)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^{n+\frac{1}{2}} \approx u(t^{n+\frac{1}{2}})$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t^{n+\frac{1}{2}}:= t^n + \frac{\tau }{2}$$\end{document} . The latter term can in turn be approximated to second order by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{2} \left( u^n + u^{n+1}\right) $$\end{document} .

We first examine the SQP cost functional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {\tilde{J}}$$\end{document} defined in (14). Starting with the integrals involving \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( u-\hat{u}\right) ^2$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( v-\hat{v}\right) ^2$$\end{document} , we can apply the trapezoidal quadrature rule to obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\alpha _1}{2} \int _Q \left( u-\hat{u}\right) ^2 \, d Q \approx \frac{\tau \alpha _1}{4} \int _\varOmega \Big [(u^0 - \hat{u}^0)^2 + 2\sum _{i=1}^{N_t-1}(u^i - \hat{u}^i)^2 + (u^{N_t} - \hat{u}^{N_t})^2 \Big ] \, dx, \end{aligned}$$\end{document}

and similarly for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( v-\hat{v}\right) ^2$$\end{document} . The control and adjoint variables need to be discretized at half steps, i.e., at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{2},...,N_t-\frac{1}{2}$$\end{document} . For the control variables, we apply the midpoint rule to the integrals involving \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a^2$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b^2$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{\beta _1}{2} \int _Q a^2 \, d Q \approx \frac{\tau \beta _1}{2} \int _\varOmega \sum _{i=0}^{N_t-1} (a^{i+\frac{1}{2}})^2 \, dx, \end{aligned}$$\end{document}

and similarly for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b^2$$\end{document} . The rest of the terms are integrals of polynomials involving adjoint variables coupled with state variables. Here, we also apply the midpoint rule, and to achieve an analogous time-stepping to the OTD approach, we apply a further second-order approximation in time. We illustrate the application of these steps on the following term, with other calculations following similarly:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \int _Q&\gamma \left( q_k - p_k\right) v_k \left( u-u_k\right) ^2 \, d Q \\&\approx \int _\varOmega \tau \gamma \sum _{i=0}^{N_t-1} (q_k^{i+\frac{1}{2}} - p_k^{i+\frac{1}{2}})v_k^{i+\frac{1}{2}} (u^{i+\frac{1}{2}}-u_k^{i+\frac{1}{2}} )^2 \, dx \\&\approx \int _\varOmega \frac{\tau \gamma }{2} \left( (q_k^{\frac{1}{2}} - p_k^{\frac{1}{2}}) v_k^0 ( u^0-u_k^0 ) ^2 + \sum _{i=1}^{N_t-1} \left[ (q_k^{i-\frac{1}{2}} - p_k^{i-\frac{1}{2}}) \right. \right. \\&\hspace{3mm} {}+ \left. \left. (q_k^{i+\frac{1}{2}} - p_k^{i+\frac{1}{2}}) \right] v_k^i (u^i-u_k^i )^2 + (q_k^{N_t-\frac{1}{2}} - p_k^{N_t-\frac{1}{2}}) v_k^{N_t} ( u^{N_t}-u_k^{N_t} ) ^2 \right) \, dx. \end{aligned} \end{aligned}$$\end{document}

Next, we approximate the constraints using a similar two-step scheme as the one above. For instance, approximating the time derivative within the first state equation in (15) yields

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&u^{n+1} - u^n - \tau D_u \nabla ^2 u^{n+\frac{1}{2}} + \tau \gamma \big ( 1-2u_k^{n+\frac{1}{2}} v_k^{n+\frac{1}{2}}\big ) u^{n+\frac{1}{2}} - \tau \gamma (u_k^{n+\frac{1}{2}} )^2 v^{n+\frac{1}{2}}\\&\quad -\tau \gamma a^{n+\frac{1}{2}} + 2 \tau \gamma (u_k^{n+\frac{1}{2}} )^2 v_k^{n+\frac{1}{2}} \approx 0, \end{aligned}$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=0,1,...,N_t-1$$\end{document} . Then, approximating the terms that involve state variables at half-time-steps using the trapezoidal rule, we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&u^{n+1} - u^n - \frac{\tau D_u}{2} \nabla ^2 \left( u^n + u^{n+1}\right) + \frac{\tau \gamma }{2} \left[ (1-2u_k^n v_k^n)u^n \right. \\&+\left. (1-2u_k^{n+1} v_k^{n+1})u^{n+1}\right] -\frac{\tau \gamma }{2} \left[ \left( u_k^n \right) ^2 v^n + (u_k^{n+1})^2 v^{n+1} \right] \\&-\tau \gamma a^{n+\frac{1}{2}} + \tau \gamma \left[ \left( u_k^n \right) ^2 v_k^n + (u_k^{n+1})^2 v_k^{n+1}\right] \approx 0. \end{aligned} \end{aligned}$$\end{document}

This leads to integer time-steps for the state variables and half time-steps for the control variables. For each n we multiply the discretized equation by the multiplier \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^{n+\frac{1}{2}}$$\end{document} , and add the integral of this term to the cost functional. We repeat the same process with the second equation in (15) with the corresponding multiplier \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q^{n+\frac{1}{2}}$$\end{document} . Then we discretize the Lagrangian using the finite element method and obtain the fully discrete problem with the mass, stiffness, and pseudo-mass matrices as in Section 4.2.

To obtain the optimality conditions from the discrete Lagrangian, we compute its partial derivatives. Differentiating with respect to the adjoint variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{p}^{\frac{1}{2}},..., \boldsymbol{p}^{N_t-\frac{1}{2}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{q}^{\frac{1}{2}},..., \boldsymbol{q}^{N_t-\frac{1}{2}}$$\end{document} gives the same system of discretized state equations as (30) and (31), respectively. Differentiation with respect to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}^1,..., {\boldsymbol{u}^{N_t-1}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{v}^1,..., {\boldsymbol{v}^{N_t-1}}$$\end{document} yields adjoint equations (32)–(33), while differentiating with respect to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}^0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{v}^0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}^{N_t}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{v}^{N_t}$$\end{document} gives the corresponding adjoint equations at first and last half-step, respectively. We also recover the gradient equations as in (34). Not obtaining terms with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{p}^{N_t},\boldsymbol{q}^{N_t}$$\end{document} is equivalent to applying the zero final-time condition in the Optimize-then-Discretize approach.

All-at-Once System

We have constructed a methodology by which the linearized systems arising from the OTD and DTO approaches coincide. To solve for the state, control, and adjoint variables, we can now construct the all-at-once system matrix. For a more convenient notation, we introduce the matrices

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \textbf{L}_{(1)}^i := \frac{D_u}{2} \textbf{K} + \frac{\gamma }{2} \textbf{M} - \gamma \textbf{M}_{u^iv^i}, \hspace{2mm} \textbf{L}_{(2)}^i := \frac{D_v}{2} \textbf{K} + \frac{\gamma }{2} \textbf{M}_{(u^i)^2}, \hspace{2mm} {\text { for }} \; i=0,...,N_t, \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \textbf{A}_{(1)}^i := {\left\{ \begin{array}{ll} \gamma \textbf{M}_{v^0(q^\frac{1}{2} - p^\frac{1}{2})}, & i=0, \\ \gamma \left( \textbf{M}_{v^i(q^{i-\frac{1}{2}} - p^{i-\frac{1}{2}})} + \textbf{M}_{v^i(q^{i+\frac{1}{2}} - p^{i+\frac{1}{2}})} \right) , & i=1,...,N_t-1, \\ \gamma \textbf{M}_{v^{N_t}(q^{N_t-\frac{1}{2}} - p^{N_t-\frac{1}{2}})}, & i=N_t, \end{array}\right. } \\ \textbf{A}_{(12)}^i := {\left\{ \begin{array}{ll} \gamma \textbf{M}_{u^0(q^\frac{1}{2} - p^\frac{1}{2})}, & i=0, \\ \gamma \left( \textbf{M}_{u^i(q^{i-\frac{1}{2}} - p^{i-\frac{1}{2}})} + \textbf{M}_{u^i(q^{i+\frac{1}{2}} - p^{i+\frac{1}{2}})} \right) , & i=1,...,N_t-1, \\ \gamma \textbf{M}_{u^{N_t}(q^{N_t-\frac{1}{2}} - p^{N_t-\frac{1}{2}})}, & i=N_t. \end{array}\right. } \end{aligned} \end{aligned}$$\end{document}

To devise an efficient iterative solver, we transform the system matrix into a simpler and more convenient form. Firstly, using the initial conditions for the state equations and making use of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{A}_{(1)}^0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{A}_{(12)}^0$$\end{document} , we can separate the adjoint equations at the first half time-step and obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \textbf{M} \boldsymbol{p}^0&= \frac{\tau \alpha _1}{2} \left( \boldsymbol{\hat{u}}^0 - \boldsymbol{u}_0 \right) - \left[ -\textbf{M} + \tau \textbf{L}_{(1)}^0 \right] \boldsymbol{p}^\frac{1}{2} - \tau \gamma \textbf{M}_{u^0v^0} \boldsymbol{q}^\frac{1}{2}, \\ \textbf{M} \boldsymbol{q}^0&= \frac{\tau \alpha _2}{2} \left( \boldsymbol{\hat{v}}^0 - \boldsymbol{v}_0 \right) + \frac{\tau \gamma }{2} \textbf{M}_{( u^0)^2}\boldsymbol{p}^\frac{1}{2} - \left[ -\textbf{M} + \tau \textbf{L}_{(2)}^0 \right] \boldsymbol{q}^\frac{1}{2}, \end{aligned} \end{aligned}$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{p}^0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{q}^0$$\end{document} . In block notation, the remaining all-at-once system can be written as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{bmatrix} \textbf{A}_1 & \textbf{0} & \textbf{B}_1^\top \\ \textbf{0} & \textbf{A}_2 & \textbf{B}_2^\top \\ \textbf{B}_1 & \textbf{B}_2 & \textbf{0} \end{bmatrix} \begin{bmatrix} \underline{\boldsymbol{u}} \\ \underline{\boldsymbol{a}} \\ \underline{\boldsymbol{p}} \end{bmatrix} = \begin{bmatrix} \underline{\boldsymbol{c}} \\ \underline{\boldsymbol{0}} \\ \underline{\boldsymbol{d}} \end{bmatrix}, \end{aligned}$$\end{document}

where the vectors for the fully discretized state, adjoint and control variables are defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \underline{\boldsymbol{u}}&:= \left[ (\boldsymbol{u}^1)^{\top }, ..., (\boldsymbol{u}^{N_t})^{\top }, (\boldsymbol{v}^1)^{\top }, ..., (\boldsymbol{v}^{N_t})^{\top }\right] ^{\top }, \\ \underline{\boldsymbol{p}}&:= \left[ (\boldsymbol{p}^{\frac{1}{2}})^{\top }, ..., (\boldsymbol{p}^{N_t-\frac{1}{2}})^{\top }, (\boldsymbol{q}^{\frac{1}{2}})^{\top }, ..., (\boldsymbol{q}^{N_t-\frac{1}{2}})^{\top }\right] ^{\top },\\ \underline{\boldsymbol{a}}&:= \left[ (\boldsymbol{a}^{\frac{1}{2}})^{\top }, ..., (\boldsymbol{a}^{N_t-\frac{1}{2}})^{\top }, (\boldsymbol{b}^{\frac{1}{2}})^{\top }, ..., (\boldsymbol{b}^{N_t-\frac{1}{2}})^{\top }\right] ^{\top }, \end{aligned} \end{aligned}$$\end{document}

and the right-hand side vectors are defined by the right-hand side terms derived in Section 4.2, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\underline{\boldsymbol{c}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\underline{\boldsymbol{d}}$$\end{document} corresponding to the adjoint and state equations, respectively. The block-matrices in (36) are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \textbf{A}_1 :=&\begin{bmatrix} \textbf{A}_{11} & \textbf{A}_{12} \\ \textbf{A}_{21} & \textbf{A}_{22} \end{bmatrix}, \hspace{5mm} \textbf{A}_2:={\text {blkdiag}}\big ( \left\{ \tau \beta _1 \textbf{M}\right\} _{i=1}^{N_t} , \left\{ \tau \beta _2 \textbf{M}\right\} _{i=1}^{N_t} \big ) , \\ \textbf{B}_1 :=&\begin{bmatrix} \textbf{B}_{11} & \textbf{B}_{12} \\ \textbf{B}_{21} & \textbf{B}_{22} \end{bmatrix}, \hspace{5mm} \textbf{B}_2:={\text {blkdiag}}\big ( \left\{ -\tau \gamma \textbf{M}\right\} _{i=1}^{2N_t}\big ). \end{aligned} \end{aligned}$$\end{document}

The sub-blocks of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{A}_1$$\end{document} are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \textbf{A}_{11}&:= {\text {blkdiag}} \left( \left\{ \tau \alpha _1 \textbf{M} + \tau \textbf{A}_{(1)}^i \right\} _{i=1}^{N_t-1}, \frac{\tau \alpha _1}{2} \textbf{M} + \tau \textbf{A}_{(1)}^{N_t} \right) , \\ \textbf{A}_{12} \equiv \textbf{A}_{21}&:= {\text {blkdiag}} \left( \left\{ \tau \textbf{A}_{(12)}^i \right\} _{i=1}^{N_t} \right) , \\ \textbf{A}_{22}&:= {\text {blkdiag}} \left( \left\{ \tau \alpha _2 \textbf{M} \right\} _{i=1}^{N_t-1}, \frac{\tau \alpha _2}{2} \textbf{M} \right) , \end{aligned} \end{aligned}$$\end{document}

while the (lower-triangular) sub-blocks of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{B}_1$$\end{document} are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \textbf{B}_{11}&:= {\text {blkdiag}} \left( \left\{ \textbf{M} + \tau \textbf{L}_{(1)}^i \right\} _{i=1}^{N_t} \right) + {\text {blksubdiag}} \left( \left\{ -\textbf{M} + \tau \textbf{L}_{(1)}^i\right\} _{i=1}^{N_t-1} \right) , \\ \textbf{B}_{12}&:= {\text {blkdiag}} \left( \left\{ -\frac{\tau \gamma }{2} \textbf{M}_{(u^i)^2} \right\} _{i=1}^{N_t} \right) + {\text {blksubdiag}} \left( \left\{ -\frac{\tau \gamma }{2} \textbf{M}_{(u^i)^2} \right\} _{i=1}^{N_t-1} \right) , \\ \textbf{B}_{21}&:= {\text {blkdiag}} \left( \left\{ \tau \gamma \textbf{M}_{u^iv^i} \right\} _{i=1}^{N_t} \right) + {\text {blksubdiag}} \left( \left\{ \tau \gamma \textbf{M}_{u^iv^i} \right\} _{i=1}^{N_t-1} \right) , \\ \textbf{B}_{22}&:= {\text {blkdiag}} \left( \left\{ \textbf{M} + \tau \textbf{L}_{(2)}^i \right\} _{i=1}^{N_t} \right) + {\text {blksubdiag}} \left( \left\{ -\textbf{M} + \tau \textbf{L}_{(2)}^i \right\} _{i=1}^{N_t-1} \right) . \end{aligned} \end{aligned}$$\end{document}

We can now eliminate the discretized gradient equations and reorder the equations and variables. With two state and two adjoint PDEs, we can think of the system matrix as a matrix composed of four blocks. For computational convenience, we wish for the (1, 1)-block of the matrix to be positive semi-definite. This can be implemented by absorbing the negative sign into part of the vector of unknowns and by negating the (2, 2)-block of the matrix and part of the right-hand side. With these modifications, we obtain the system

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \underbrace{\begin{bmatrix} \textbf{A} & \textbf{B}^\top \\ \textbf{B} & -\textbf{C} \end{bmatrix}}_{\boldsymbol{\mathcal {A}}} \underbrace{\begin{bmatrix} -\underline{\boldsymbol{p}} \\ \underline{\boldsymbol{u}} \end{bmatrix}}_{\underline{\boldsymbol{w}}} := \begin{bmatrix} \textbf{A}_3 & \textbf{B}_1 \\ \textbf{B}_1^\top & -\textbf{A}_1 \end{bmatrix} \begin{bmatrix} -\underline{\boldsymbol{p}} \\ \underline{\boldsymbol{u}} \end{bmatrix} = \underbrace{\begin{bmatrix} \underline{\boldsymbol{d}} \\ -\underline{\boldsymbol{c}} \end{bmatrix}}_{\underline{\boldsymbol{b}}}, \end{aligned}$$\end{document}

where the block \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{A}_3:= \textbf{B}_2 \textbf{A}_2^{-1} \textbf{B}_2^\top $$\end{document} arises by incorporating the gradient equations:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{A}_3:={\text {blkdiag}}\left( \left\{ \frac{\tau \gamma ^2}{\beta _1} \textbf{M}\right\} _{i=1}^{N_t} , \left\{ \frac{\tau \gamma ^2}{\beta _2}\textbf{M}\right\} _{i=1}^{N_t} \right) . \end{aligned}$$\end{document}

We highlight that the system matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal {A}}$$\end{document} is symmetric indefinite and the main challenge is to find a solution to the large-scale (generalized) saddle-point problem (37). After solving the system above and simplifying the terms, we solve (35) for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{p}^0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{q}^0$$\end{document} using Chebyshev semi-iteration [33, 34] with Jacobi splitting.

Preconditioning

Solving large (generalized) saddle-point problems of the form (37) is a topic of wide interest in the numerical analysis community. Since finding an “exact” solution to this problem by using a direct solver is computationally too expensive, we make use of iterative solvers known to be suitable for this type of problem structure, specifically Krylov subspace solvers. We need to ensure that we find the solutions of the systems robustly and efficiently, that the iterative solver will converge fast, and we take advantage of the problem structures within the algorithms used for their solution.

To address the first issue, we note that the convergence rate of iterative methods often depends on the distribution of eigenvalues of matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal {A}}$$\end{document} . As \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal {A}}$$\end{document} is considered to be indefinite, and its eigenvalues can be distributed across many orders of magnitude, the convergence of an iterative method can be slow. By constructing a preconditioner \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\boldsymbol{\mathcal {P}}}$$\end{document} , we can transform the linear system (37) into the left-preconditioned system \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal {P}}^{-1} \boldsymbol{\mathcal {A}}\underline{\boldsymbol{w}} = \boldsymbol{\mathcal {P}}^{-1}\underline{\boldsymbol{b}}$$\end{document} or the right-preconditioned system \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\boldsymbol{\mathcal {A}}\boldsymbol{\mathcal {P}}^{-1}) (\boldsymbol{\mathcal {P}}\underline{\boldsymbol{w}}) = \underline{\boldsymbol{b}}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal {P}}^{-1} \boldsymbol{\mathcal {A}}$$\end{document} (equivalently \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal {A}}\boldsymbol{\mathcal {P}}^{-1}$$\end{document} ) has improved properties that ensure faster convergence of iterative methods. An “ideal” preconditioner for (generalized) saddle-point systems (37), given certain properties of matrices \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{A}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{B}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{C}$$\end{document} , is of the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \boldsymbol{\mathcal {P}}= \begin{bmatrix} \textbf{A} & \textbf{0} \\ \textbf{0} & \textbf{S} \end{bmatrix}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}$$\end{document} denotes the (negative) Schur complement \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}:= \textbf{C} + \textbf{B} \textbf{A}^{-1}\textbf{B}^\top $$\end{document} . Specifically, we assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{A}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}$$\end{document} are symmetric positive definite so that the preconditioner may be applied within the MINRES algorithm [35]. If \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{C}$$\end{document} is also symmetric positive semi-definite, we would have that the eigenvalues of the preconditioned matrix are contained within \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[-1, \frac{1}{2}(1-\sqrt{5})] \cup [1, \frac{1}{2}(1+\sqrt{5})]$$\end{document} regardless of problem dimension; see e.g., [36, Cor. 1], [21, Thm. 4], and [37, Lem. 2.2]. Now, due to the structure of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal {A}}$$\end{document} , applying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}^{-1}$$\end{document} can be very expensive due to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}$$\end{document} typically being a dense matrix. To offer a more favourable alternative to applying \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}^{-1}$$\end{document} explicitly, one option is to approximate the Schur complement using a suitable product of matrices, the so-called “matching strategy” which was devised for simpler PDE-constrained optimization problems in [38, 39]. We start by approximating the Schur complement with

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{S} \approx \hat{\textbf{S}} := (\textbf{B}+\textbf{D}) \textbf{A}^{-1} (\textbf{B}+\textbf{D})^\top , \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D}$$\end{document} is some real matrix to be determined. Since we can rewrite the approximated Schur complement as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\textbf{S}} = [ \left( \textbf{B}+\textbf{D}\right) \textbf{A}^{-\frac{1}{2}} ] [ \left( \textbf{B}+\textbf{D}\right) \textbf{A}^{-\frac{1}{2}} ]^\top $$\end{document} , we find that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\textbf{S}}$$\end{document} must be symmetric positive semi-definite, and provided \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D}$$\end{document} is suitably chosen (as below) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\textbf{S}}$$\end{document} will be symmetric positive definite. As this property also holds for the block \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{A}$$\end{document} , the preconditioner is then symmetric positive definite, which allows us to consider more options for iterative solvers shown later in this section. Expanding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\textbf{S}}$$\end{document} and matching the terms with those in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}$$\end{document} , we determine \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D}$$\end{document} by specifying that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D} \textbf{A}^{-1} \textbf{D}^\top $$\end{document} should approximate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{C}$$\end{document} . We can again rewrite this quantity as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$[ \textbf{D} \textbf{A}^{-\frac{1}{2}} ] [ \textbf{D} \textbf{A}^{-\frac{1}{2}} ]^\top $$\end{document} , so the approximation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{C}$$\end{document} must be a symmetric and positive semi-definite matrix. A suitable approximation of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{C}$$\end{document} can be a block diagonal matrix, incorporating components that include mass matrices, so that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{C} \approx \hat{\textbf{C}} = {\text {blkdiag}} \left( \left\{ \tau \alpha _1 \textbf{M} \right\} _{i=1}^{N_t-1}, \frac{\tau \alpha _1}{2} \textbf{M}, \left\{ \tau \alpha _2 \textbf{M} \right\} _{i=1}^{N_t-1}, \frac{\tau \alpha _2}{2} \textbf{M} \right) . \end{aligned}$$\end{document}

This is a convenient choice for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\textbf{C}}$$\end{document} because it captures many of the (linear) terms within the PDEs, and is also positive definite and hence approximable by a term of the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D} \textbf{A}^{-1} \textbf{D}^\top $$\end{document} . Now we need to find the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D}$$\end{document} . Since the block \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{A}$$\end{document} , and consequently, its inverse, is block diagonal, we obtain that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D}$$\end{document} is a block diagonal matrix of the form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{D}:={\text {blkdiag}}\left( \left\{ \textbf{D}_{(1)}^i \right\} _{i=1}^{N_t} , \left\{ \textbf{D}_{(2)}^i\right\} _{i=1}^{N_t} \right) , \end{aligned}$$\end{document}

and we calculate that the diagonal blocks of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D}$$\end{document} must be

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \textbf{D}_{(1)}^i&:= \tau \gamma \sqrt{\frac{\alpha _1}{\beta _1}} \textbf{M} \quad {\text { for}} \; i=1,..., N_t-1, \hspace{5mm}&\textbf{D}_{(1)}^{N_t}&:= \tau \gamma \sqrt{\frac{\alpha _1}{2\beta _1}} \textbf{M}, \\ \textbf{D}_{(2)}^i&:= \tau \gamma \sqrt{\frac{\alpha _2}{\beta _2}} \textbf{M} \quad {\text { for}} \; i=1,..., N_t-1, \hspace{5mm}&\textbf{D}_{(2)}^{N_t}&:= \tau \gamma \sqrt{\frac{\alpha _2}{2\beta _2}} \textbf{M}. \end{aligned} \end{aligned}$$\end{document}

With an approximation of the ideal preconditioner \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\hat{\mathcal {P}}}$$\end{document} , for computational efficiency we must define the solution of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\hat{\mathcal {P}}} \underline{\boldsymbol{y}} = \underline{\boldsymbol{z}}$$\end{document} for a given \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\underline{\boldsymbol{z}}$$\end{document} without forming the inverse of the preconditioner explicitly. As \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\hat{\mathcal {P}}}$$\end{document} is block diagonal, we can split this computation into two parts. First, dealing with the (1, 1)-block of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\boldsymbol{\mathcal {P}}}$$\end{document} comes down to solving a system with scaled mass matrices on the diagonal. Instead of using a direct solver, we can use Chebyshev semi-iteration with Jacobi splitting, which provides an efficient approximation of mass matrices. For 2D problems discretized using P1 finite elements, as we later consider, we know that the eigenvalues of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{M}$$\end{document} preconditioned with its diagonal must lie between 0.5 and 2; see e.g., [40]. We use this information within our implementation of Chebyshev semi-iteration and let the method run for 20 iterations. Second, considering the approximated Schur complement within the (2, 2)-block of the preconditioner, and with vector superscripts denoting the second halves of the vectors, we now aim to solve the system \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\textbf{B}+\textbf{D}) \textbf{A}^{-1} (\textbf{B}+\textbf{D})^\top \underline{\boldsymbol{y}}^2 = \underline{\boldsymbol{z}}^2$$\end{document} . As the matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{B}+\textbf{D}$$\end{document} is block upper-triangular, we first solve \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\textbf{B}+\textbf{D}) \underline{\boldsymbol{m}} = \underline{\boldsymbol{z}}^2$$\end{document} by block-backward substitution, after which we solve \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ (\textbf{B}+\textbf{D})^\top \underline{\boldsymbol{y}}^2 = \textbf{A} \underline{\boldsymbol{m}}$$\end{document} by block-forward substitution. Each of the sub-systems can be solved using an algebraic multigrid solver (see [41], for instance), and we apply 6 V-cycles to ensure a sufficiently accurate approximation.

Finally, we highlight that there are a variety of preconditioned Krylov subspace iterative solvers for different classes of matrices and preconditioners. In particular, for symmetric matrices with symmetric and positive definite preconditioners of the structure (38), the MINRES solver is the preferred choice. One could also apply the Bramble–Pasciak Conjugate Gradient Method [42] with a modified preconditioner, and GMRES [43], BICG [44], or its stabilized variant BiCGSTAB [45] for problems involving non-symmetric systems or preconditioners, for example block triangular variants of (38), which may also be readily applied. A key advantage of using such Krylov subspace solvers in general is their computational cost: they require multiplication with the system matrix, but this can be performed without having to form the matrix in full. That is, we are only required to define the linear operator that computes \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal {A}}\underline{\boldsymbol{v}}$$\end{document} for an arbitrary vector \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\underline{\boldsymbol{v}}$$\end{document} of consistent dimension. We may therefore utilize the specific block structure of the matrix under consideration, and define a matrix–vector product by simply using block-forward substitution for lower-triangular blocks and block-backward substitution for upper-triangular blocks.

Remark 1

To further motivate the above approximation of the Schur complement, consider the distribution of eigenvalues of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\textbf{S}}^{-1} \textbf{S}$$\end{document} in a “linearized” setting, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{u} = \bar{v} = \bar{p} = \bar{q} = 0$$\end{document} . Then, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{C}$$\end{document} is positive definite, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D} \equiv \textbf{C}^{\frac{1}{2}} \textbf{A}^{\frac{1}{2}}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D} \textbf{A}^{-1} \textbf{D}^\top \equiv \textbf{C}$$\end{document} , and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \hat{\textbf{S}} = \textbf{S} + \textbf{C}^{\frac{1}{2}} \textbf{A}^{-\frac{1}{2}} \textbf{B}^\top + \textbf{B} \textbf{A}^{-\frac{1}{2}} \textbf{C}^{\frac{1}{2}}. \end{aligned}$$\end{document}

By construction, both \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\textbf{S}}$$\end{document} are symmetric positive definite. Hence, following the same reasoning as in [39, Thm. 4.1], it may be shown that the eigenvalues of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{S}^{-1} S$$\end{document} are bounded below by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{1}{2}$$\end{document} , irrespective of problem dimension and parameters involved. Further, it holds that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{C}^{\frac{1}{2}} \textbf{A}^{-\frac{1}{2}} \textbf{B}^\top + \textbf{B} \textbf{A}^{-\frac{1}{2}} \textbf{C}^{\frac{1}{2}} ={}&\frac{\sqrt{\alpha _1 \beta _1}}{\gamma } \left( \begin{bmatrix} \textbf{T}_1 & \textbf{0} \\ \textbf{0} & \textbf{0} \end{bmatrix} \otimes \textbf{M} + \begin{bmatrix} \textbf{T}_2 & \textbf{0} \\ \textbf{0} & \textbf{0} \end{bmatrix} \otimes \tau \textbf{L}_{(1)} \right) \\&+ \frac{\sqrt{\alpha _2 \beta _2}}{\gamma } \left( \begin{bmatrix} \textbf{0} & \textbf{0} \\ \textbf{0} & \textbf{T}_1 \end{bmatrix} \otimes \textbf{M} + \begin{bmatrix} \textbf{0} & \textbf{0} \\ \textbf{0} & \textbf{T}_2 \end{bmatrix} \otimes \tau \textbf{L}_{(2)} \right) , \end{aligned}$$\end{document}

where for the linearized problem \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{L}_{(1)} = \frac{D_u}{2} \textbf{K} + \frac{\gamma }{2} \textbf{M}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{L}_{(2)} = \frac{D_v}{2} \textbf{K}$$\end{document} , and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{T}_1 = \begin{bmatrix} 2 & -1 \\ -1 & 2 & \ddots \\ & \ddots & \ddots & -1 \\ & & -1 & 2 & -\frac{1}{\sqrt{2}} \\ & & & -\frac{1}{\sqrt{2}} & \sqrt{2} \\ \end{bmatrix}, \hspace{5mm} \textbf{T}_2 = \begin{bmatrix} 2 & 1 \\ 1 & 2 & \ddots \\ & \ddots & \ddots & 1 \\ & & 1 & 2 & \frac{1}{\sqrt{2}} \\ & & & \frac{1}{\sqrt{2}} & \sqrt{2} \\ \end{bmatrix}. \end{aligned}$$\end{document}

As \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{L}_{(1)}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{L}_{(2)}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{T}_1$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{T}_2$$\end{document} are all at least positive semi-definite, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{C}^{\frac{1}{2}} \textbf{A}^{-\frac{1}{2}} \textbf{B}^\top + \textbf{B} \textbf{A}^{-\frac{1}{2}} \textbf{C}^{\frac{1}{2}}$$\end{document} is positive semi-definite, hence \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S} \preceq \hat{\textbf{S}}$$\end{document} and the eigenvalues of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\textbf{S}}^{-1} \textbf{S}$$\end{document} are bounded above by 1.

Comparison with Backward Euler Method

To evaluate the performance of our new numerical scheme, we compare it to the problem obtained by applying the backward Euler method for the time discretization of (7)–(8); see [11]. For the Schnakenberg model, the linearized equations used in this section are (15) and (16). To reduce the number of variables, we can substitute for the control variables using the gradient equations (6).

Applying the backward Euler method over the time interval (0, T), with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_t$$\end{document} time-steps of size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} , we begin by discretizing the initial and final-time conditions,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{M}\boldsymbol{u}^0 = \boldsymbol{u}_0, \quad \textbf{M}\boldsymbol{v}^0 = \boldsymbol{v}_0, \quad \textbf{M}\boldsymbol{p}^{N_t} = \textbf{M}\boldsymbol{q}^{N_t} = \boldsymbol{0}, \end{aligned}$$\end{document}

with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}_0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \boldsymbol{v}_0$$\end{document} defined as in Section 4.2. At each Newton iteration, discretization of the linearized PDEs leads to the following system:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} - \textbf{M}\boldsymbol{u}^n + [\textbf{M} + 2\tau \textbf{L}_{(1)}^{n+1}] \boldsymbol{u}^{n+1} - \tau \gamma \textbf{M}_{(u^{n+1})^2}\boldsymbol{v}^{n+1} - \frac{\tau \gamma ^2}{\beta _1}\textbf{M}\boldsymbol{p}^{n+1} = -\tau \boldsymbol{d}^{n+1}{}&, \\ 2\tau \gamma \textbf{M}_{u^{n+1} v^{n+1}} \boldsymbol{u}^{n+1} - \textbf{M}\boldsymbol{v}^n + [\textbf{M} + 2\tau \textbf{L}_{(2)}^{n+1}] \boldsymbol{v}^{n+1} - \frac{\tau \gamma ^2}{\beta _2}\textbf{M}\boldsymbol{q}^{n+1} = \tau \boldsymbol{d}^{n+1}{}&, \\ [\tau \alpha _1 \textbf{M} + \tau \textbf{A}_{(1)}^{E,n}] \boldsymbol{u}^n +\tau \textbf{A}_{(12)}^{E,n} \boldsymbol{v}^n + [\textbf{M} + 2\tau \textbf{L}_{(1)}^n]\boldsymbol{p}^n - \textbf{M}\boldsymbol{p}^{n+1} + 2\tau \gamma \textbf{M}_{u^n v^n} \boldsymbol{q}^n&\\ = \tau \boldsymbol{c}^{n,n}{}&, \\ \tau \textbf{A}_{(12)}^{E,n}\boldsymbol{u}^n + \tau \alpha _2 \textbf{M} \boldsymbol{v}^n - \tau \gamma \textbf{M}_{(u^n)^2}\boldsymbol{p}^n + [\textbf{M} + 2\tau \textbf{L}_{(2)}^n]\boldsymbol{q}^n - \textbf{M}\boldsymbol{q}^{n+1} = \tau \boldsymbol{h}^{n,n}{}&, \end{aligned} \end{aligned}$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n=0,1,...,N_t-1$$\end{document} . Here we re-use the right-hand side vectors and some of the matrices introduced in Sections 4.2 and 6. The matrices specific to the backward Euler scheme are

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{A}_{(1)}^{E,i}:= 2 \gamma \textbf{M}_{v^i(q^i - p^i)}, \quad \textbf{A}_{(12)}^{E,i}:=2 \gamma \textbf{M}_{u^i(q^i - p^i)}, \quad {\text { for }} \; i=0,1,...,N_t-1. \end{aligned}$$\end{document}

As this system contains equations evolving forward and backward in time, we again solve for all the discretized variables at once until a stopping criterion for the Newton iterates is achieved. That is, we need to find solutions to the state variables u and v at time-steps \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1,...,N_t$$\end{document} and to the adjoint variables p and q at time-steps \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0,1,...,N_t-1$$\end{document} .

We construct the system matrix for (39), involving all four variables at time-steps \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0,1,...,N_t$$\end{document} . Notice that we can exclude the variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{p}^0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{q}^0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}^{N_t}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{v}^{N_t}$$\end{document} from the all-at-once system, which can be done by decoupling the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_t$$\end{document} -th equation of each of the state equations and the first equation of each of the adjoint equations. To obtain a symmetric system matrix, we reorder the variables, and using similar arguments as in Section 7, we solve for the vectors of unknowns \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left( -\underline{\boldsymbol{p}}, \underline{\boldsymbol{u}} \right) $$\end{document} , where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \underline{\boldsymbol{u}}&:= \left[ (\boldsymbol{u}^1)^{\top }, ..., (\boldsymbol{u}^{N_t-1})^{\top }, (\boldsymbol{v}^1)^{\top }, ..., (\boldsymbol{v}^{N_t-1})^{\top }\right] ^{\top }, \\ \underline{\boldsymbol{p}}&:= \left[ (\boldsymbol{p}^1)^{\top }, ..., (\boldsymbol{p}^{N_t-1})^{\top }, (\boldsymbol{q}^1)^{\top }, ..., (\boldsymbol{q}^{N_t-1})^{\top }\right] ^{\top }, \end{aligned} \end{aligned}$$\end{document}

and with the system matrix

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \boldsymbol{\mathcal {A}} := \begin{bmatrix} \textbf{A}^E & \begin{bmatrix} (\textbf{B}_{11}^E)^\top & \textbf{B}_{21}^E \\ \textbf{B}_{12}^E & (\textbf{B}_{22}^E)^\top \end{bmatrix} \\ \begin{bmatrix} \textbf{B}_{11}^E & \textbf{B}_{12}^E \\ \textbf{B}_{21}^E & \textbf{B}_{22}^E \end{bmatrix} & \begin{bmatrix} -\textbf{C}_{11}^E & -\textbf{C}_{12}^E \\ -\textbf{C}_{12}^E & -\textbf{C}_{22}^E \end{bmatrix} \end{bmatrix} = \begin{bmatrix} \textbf{A}^E & \left( \textbf{B}^E\right) ^\top \\ \textbf{B}^E & -\textbf{C}^E \end{bmatrix} . \end{aligned}$$\end{document}

The (1, 1)-block and most of the sub-blocks are block diagonal:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \textbf{A}^E&:= {\text {blkdiag}} \left( \hspace{-0.1 cm} \left\{ \frac{\tau \gamma ^2}{\beta _1} \textbf{M}\right\} _{i=1}^{N_t-1} , \left\{ \frac{\tau \gamma ^2}{\beta _2} \textbf{M}\right\} _{i=1}^{N_t-1} \right) , \\ \textbf{C}_{11}^E&:= {\text {blkdiag}} \left( \hspace{-0.1 cm} \left\{ \tau \alpha _1 \textbf{M} + \tau \textbf{A}_{(1)}^{E,i}\right\} _{i=1}^{N_t-1}\right) , \end{aligned} \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \textbf{B}_{12}^E&:= {\text {blkdiag}} \left( \left\{ 2\tau \gamma \textbf{M}_{u^iv^i}\right\} _{i=1}^{N_t-1}\right) ,&\textbf{C}_{12}^E&:= {\text {blkdiag}} \left( \left\{ \tau \textbf{A}_{(12)}^{E,i}\right\} _{i=1}^{N_t-1}\right) ,\\ \textbf{B}_{21}^E&:= {\text {blkdiag}} \left( \left\{ -\tau \gamma \textbf{M}_{(u^i)^2}\right\} _{i=1}^{N_t-1} \right) ,&\textbf{C}_{22}^E&:= {\text {blkdiag}} \left( \left\{ \tau \alpha _2 \textbf{M}\right\} _{i=1}^{N_t-1}\right) . \end{aligned} \end{aligned}$$\end{document}

Two of the sub-blocks appearing in the (1, 2)- and (2, 1)-blocks are upper-triangular:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \textbf{B}_{11}^E&:= {\text {blkdiag}} \left( \left\{ \textbf{M} + 2\tau \textbf{L}_{(1)}^i \right\} _{i=1}^{N_t-1} \right) + {\text {blksupdiag}} \left( \left\{ -\textbf{M} \right\} _{i=1}^{N_t-2} \right) , \\ \textbf{B}_{22}^E&:= {\text {blkdiag}} \left( \left\{ \textbf{M} + 2\tau \textbf{L}_{(2)}^i \right\} _{i=1}^{N_t-1} \right) + {\text {blksupdiag}} \left( \left\{ -\textbf{M} \right\} _{i=1}^{N_t-2} \right) . \end{aligned} \end{aligned}$$\end{document}

Additionally, we need to solve the equations that are excluded from the system above. We solve the coupled system for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{p}^0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{q}^0$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned}&\begin{bmatrix} \textbf{M} + 2\tau \textbf{L}_{(1)}^0 & 2\tau \gamma \textbf{M}_{u^0v^0} \\ -\tau \gamma \textbf{M}_{(u^0)^2} & \textbf{M} + 2\tau \textbf{L}_{(2)}^0 \end{bmatrix} \begin{bmatrix} \boldsymbol{p}^0 \\ \boldsymbol{q}^0 \end{bmatrix}\\&\quad = \begin{bmatrix} \tau \boldsymbol{c}^{0,0} + \textbf{M}\boldsymbol{p}^1 - \tau [ \alpha _1 \textbf{M} + \textbf{A}_{(1)}^{E,0} ] \boldsymbol{u}_0 - \tau \textbf{A}_{(12)}^{E,0}\boldsymbol{v}_0 \\ \tau \boldsymbol{h}^{0,0} + \textbf{M}\boldsymbol{q}^1 - \tau \textbf{A}_{(12)}^{E,0} \boldsymbol{u}_0 - \tau \alpha _2 \textbf{M} \boldsymbol{v}_0 \end{bmatrix}, \end{aligned}$$\end{document}

and another system for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}^{N_t}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{v}^{N_t}$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{bmatrix} \textbf{M} + 2\tau \textbf{L}_{(1)}^{N_t} & -\tau \gamma \textbf{M}_{(u^{N_t})^2} \\ 2\tau \gamma \textbf{M}_{u^{N_t}v^{N_t}} & \textbf{M} + 2\tau \textbf{L}_{(2)}^{N_t} \end{bmatrix} \begin{bmatrix} \boldsymbol{u}^{N_t} \\ \boldsymbol{v}^{N_t} \end{bmatrix} =\begin{bmatrix} -\tau \boldsymbol{d}^{N_t} + \textbf{M}\boldsymbol{u}^{N_t-1} \\ \tau \boldsymbol{d}^{N_t} + \textbf{M}\boldsymbol{v}^{N_t-1} \end{bmatrix}, \end{aligned}$$\end{document}

using the vectors \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{p}^1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{q}^1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}^{N_t-1}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{v}^{N_t-1}$$\end{document} previously calculated.

The matrix \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{\mathcal {A}}$$\end{document} is again of (generalized) saddle-point form, and we proceed using a similar approach as in Section 7 and construct a block-diagonal preconditioner. In particular, the (1, 1)-block of the preconditioner should approximate the (1, 1)-block of the system matrix, which only has factors of mass matrices on the block diagonal. In the numerical implementation, we approximate the mass matrices using Chebyshev semi-iteration, as in Section 7. The (2, 2)-block of the preconditioner is the approximation of the (negative) Schur complement \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}^E:= \textbf{C}^E + \textbf{B}^E \left( \textbf{A}^E \right) ^{-1} \left( \textbf{B}^E\right) ^\top $$\end{document} . The block \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{C}^E$$\end{document} can be approximated by dropping the matrices emerging from nonlinear terms, leaving \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\textbf{C}}^E = {\text {blkdiag}} \left( \left\{ \tau \alpha _1 \textbf{M} \right\} _{i=1}^{N_t-1} \left\{ \tau \alpha _2 \textbf{M} \right\} _{i=1}^{N_t-1} \right) $$\end{document} . We can thus approximate the Schur complement by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{S}^E \approx \hat{\textbf{S}}^E = (\textbf{D}^E + \textbf{B}^E) \left( \textbf{A}^E \right) ^{-1} \left( \textbf{D}^E+ \textbf{B}^E\right) ^\top $$\end{document} , and we find \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D}^E$$\end{document} by matching \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D}^E \left( \textbf{A}^E \right) ^{-1} \left( \textbf{D}^E\right) ^\top = \hat{\textbf{C}}^E$$\end{document} and neglecting the mixed terms after factorization of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{\textbf{S}}^E$$\end{document} , to obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \textbf{D}^E := {\text {blkdiag}} \left( \left\{ \tau \gamma \sqrt{\frac{\alpha _1}{\beta _1}} \textbf{M}\right\} _{i=1}^{N_t-1} \left\{ \tau \gamma \sqrt{\frac{\alpha _2}{\beta _2}} \textbf{M}\right\} _{i=1}^{N_t-1} \right) . \end{aligned}$$\end{document}

Since the system matrix is symmetric and positive definite, we use the preconditioned MINRES method to approximately solve the system. The equations for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{p}^0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{q}^0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}^{N_t}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{v}^{N_t}$$\end{document} can be solved with a direct solver, due to the relatively modest dimension of this system.

Numerical Results

To validate our numerical method, we now examine its performance on two relevant test problems. First, in Section 9.1, we benchmark our method and compare it with the backward Euler method described in Section 8, using a problem with an analytical solution. In Section 9.2, we validate our approach for parameter identification in the context of a data-driven problem, with target states generated using the Schnakenberg equations with chosen functions a and b.

In all the following experiments we consider \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha := \alpha _1 = \alpha _2$$\end{document} and, without loss of generality, set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha = 1$$\end{document} , and define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta :=\beta _1=\beta _2$$\end{document} . Further, we choose the diffusion parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_u=1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D_v=10$$\end{document} . All the simulations are performed on a unit square domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varOmega =(0,1)^2$$\end{document} . We consider zero-flux boundary conditions due to their more realistic physical significance, however, the Dirichlet boundary conditions can also be easily applied by adjusting the linear system.

To generate the matrices and the right-hand side vectors that emerge in the schemes, we use FEniCS [46] (version 2019.1.0), in particular the DOLFIN component. For other components of our implementation, we use Python [47] (v. 3.8.15) and its libraries: scipy [48] (v. 1.6.3), numpy [49] (v. 1.22.3), matplotlib [50] (v. 3.4.0), and pyamg [51] (v. 4.1.0).

Using DOLFIN, we implement a finite element method in two dimensions with P1 elements and adhere to the corresponding node ordering. Once the matrices and vectors are generated, we convert them to Python arrays. The scipy.sparse library allows us to efficiently store the sparse matrices, while dense arrays and vectors are stored as numpy arrays. We also use the MINRES implementation from this library, which we provide with routines applying the system matrix and inverse of the preconditioner, both implemented as linear operators for memory efficiency. To perform the multigrid calculations in the Schur complement, we use the smoothed_aggregation_solver from the pyamg library. The plots in the following sections are generated with the matplotlib.pyplot library.

Convergence scaling

Table 1. Störmer–Verlet method for the Schnakenberg model: weighted errors, mean number of MINRES iterations and total SQP iterations, and CPU times, for different mesh sizes and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-2}$$\end{document} . Convergence rates, computed by evaluating the ratio of successive error quantities, are given in bold i DoF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{{\text {it}}}_M$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {it}}_{{\text {SQP}}}$$\end{document} CPU124,2008.73e-28.55e-28.64e-36.70e-3256165––––2176,4002.10e-22.04e-22.10e-31.66e-3304542 4.16

4.20

4.12

4.04 31,344,8005.01e-34.96e-35.10e-44.12e-43241,202 4.19

4.10

4.11

4.03 410,497,6001.39e-31.37e-31.25e-41.03e-4109419,233 3.61

3.62

4.08

3.99

Table 2. Störmer–Verlet method for the Schnakenberg model: weighted errors, mean number of MINRES iterations and total SQP iterations, and CPU times, for different mesh sizes and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-3}$$\end{document} . Convergence rates are given in bold i DoF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{{\text {it}}}_M$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {it}}_{{\text {SQP}}}$$\end{document} CPU 124,2004.61e-12.04e-16.18e-32.92e-3306169––––2176,4001.69e-17.24e-21.73e-38.44e-4385410 2.73

2.82

3.58

3.46 31,344,8005.29e-22.25e-24.66e-42.37e-43851,745 3.19

3.22

3.70

3.56 410,497,6001.39e-25.99e-31.15e-46.06e-574516,605 3.81

3.76

4.06

3.91

Table 3. Backward Euler method for the Schnakenberg model: weighted errors, mean number of MINRES iterations and total SQP iterations, and CPU times, for different mesh sizes and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-2}$$\end{document} . Convergence rates are given in bold i DoF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{{\text {it}}}_M$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {it}}_{{\text {SQP}}}$$\end{document} CPU123,7161.03e-19.53e-28.13e-36.90e-330493––––2351,0362.47e-22.25e-21.96e-31.74e-3304500 4.14

4.23

4.14

3.96 35,372,4765.92e-35.47e-34.77e-44.30e-4254 3,672 4.18

4.11

4.11

4.05 483,954,5561.23e-31.29e-31.13e-49.90e-527440,770 4.81

4.25

4.21

4.35

Table 4. Backward Euler method for the Schnakenberg model: weighted errors, mean number of MINRES iterations and total SQP iterations, and CPU times, for different mesh sizes and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-3}$$\end{document} . Convergence rates are given in bold i DoF \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q_{{\text {error}}}$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{{\text {it}}}_M$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {it}}_{{\text {SQP}}}$$\end{document} CPU123,7165.09e-12.14e-15.95e-32.38e-3296131––––2351,0361.98e-17.97e-21.76e-38.40e-4365 761 2.57

2.68

3.37

2.84 35,372,4766.49e-22.61e-24.96e-42.50e-43656,147 3.06

3.05

3.56

3.36 483,954,5561.95e-27.78e-31.39e-47.25e-534568,518 3.33

3.36

3.56

3.45

To evaluate the performance of our method, we test it along with the backward Euler scheme on a problem with the following analytical solution:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} u^{\text {sol}}(t,{\textbf {x}})&= {\text {e}}^{0.1t} \left( \kappa ({\textbf {x}})+1\right) ,&p^{\text {sol}}(t,{\textbf {x}})&= \left( {\text {e}}^{0.1t} -{\text {e}}^{0.1T} \right) \left( \kappa ({\textbf {x}}) +1 \right) , \\ v^{\text {sol}}(t,{\textbf {x}})&= {\text {e}}^{0.15t} \left( \eta ({\textbf {x}}) +1 \right) , \hspace{5mm}&q^{\text {sol}}(t,{\textbf {x}})&= \left( {\text {e}}^{0.15t} -{\text {e}}^{0.15T} \right) \left( \eta ({\textbf {x}}) +1 \right) . \end{aligned} \end{aligned}$$\end{document}

Here, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa ({\textbf {x}}) = \cos (2 \pi x_1) \cos (2 \pi x_2)$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta ({\textbf {x}}) = \cos (\pi x_1) \cos (\pi x_2)$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\textbf {x}}=(x_1, x_2)$$\end{document} . The initial conditions for the state variables are constructed using the analytical solutions at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0$$\end{document} . For the desired states, we define the functions

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \hat{u}^{\text {sol}}(t,{\textbf {x}}) ={}&\frac{1}{\alpha }\left[ -0.1{\text {e}}^{0.1t}(\kappa ({\textbf {x}}) + 1) + 8D_u\pi ^2\left( {\text {e}}^{0.1t} -{\text {e}}^{0.1T} \right) \kappa ({\textbf {x}}) + \alpha u^{\text {sol}}(t,{\textbf {x}}) \right. \\&\left. {}+ 2\gamma u^{\text {sol}}(t,{\textbf {x}}) v^{\text {sol}}(t,{\textbf {x}}) (q^{\text {sol}}(t,{\textbf {x}})-p^{\text {sol}}(t,{\textbf {x}})) + \gamma p^{\text {sol}}(t,{\textbf {x}})) \right] , \\ \hat{v}^{\text {sol}}(t,{\textbf {x}}) ={}&\frac{1}{\alpha }\left[ -0.15{\text {e}}^{0.15t}(\eta ({\textbf {x}}) + 1) + 2D_v\pi ^2\left( {\text {e}}^{0.15t} -{\text {e}}^{0.15T} \right) \eta ({\textbf {x}}) \right. \\&\left. {}+ \alpha v^{\text {sol}}(t,{\textbf {x}}) + \gamma u^{\text {sol}}(t,{\textbf {x}})^2 (q^{\text {sol}}(t,{\textbf {x}})-p^{\text {sol}}(t,{\textbf {x}})) \right] . \end{aligned} \end{aligned}$$\end{document}

Finally, we need to introduce source functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f^{\text {sol}}(t,{\textbf {x}})$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g^{\text {sol}}(t,{\textbf {x}})$$\end{document} to the right-hand side of the first and second state equations in (2), respectively:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} f^{\text {sol}}(t,{\textbf {x}}) ={}&(0.1 + \gamma )u^{\text {sol}}(t,{\textbf {x}}) + 8D_u\pi ^2 {\text {e}}^{0.1t}\kappa ({\textbf {x}}) - \gamma |u^{\text {sol}}(t,{\textbf {x}})|^2 v^{\text {sol}}(t,{\textbf {x}}) \\&{}- \frac{\gamma ^2}{\beta } p^{\text {sol}}(t,{\textbf {x}}), \\ g^{\text {sol}}(t,{\textbf {x}}) ={}&0.15v^{\text {sol}}(t,{\textbf {x}}) + 2D_v\pi ^2 {\text {e}}^{0.15t}\eta ({\textbf {x}}) + \gamma |u^{\text {sol}}(t,{\textbf {x}})|^2 v^{\text {sol}}(t,{\textbf {x}}) \\ &{}- \frac{\gamma ^2}{\beta } q^{\text {sol}}(t,{\textbf {x}}). \end{aligned} \end{aligned}$$\end{document}

In the notation of (2), this corresponds to taking \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varPhi (u,v) = \gamma (u-u^2 v) - f^{\text {sol}}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varPsi (u,v) = \gamma u^2 v - g^{\text {sol}}$$\end{document} , which does not alter the adjoint or gradient equations derived thereafter.

As both schemes yield a symmetric system matrix, we use the MINRES solver. The preconditioner is built using the approach described in Sections 7 and 8. The MINRES solver terminates when the relative residual is below \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {tol}}_{{\text {M}}}$$\end{document} . We find that for our method, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {tol}}_{{\text {M}}} = 10^{-8}$$\end{document} is sufficient for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-3}$$\end{document} , however, we report results with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {tol}}_{{\text {M}}} = 10^{-9}$$\end{document} as this tolerance was required to achieve the correct scaling of the error for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-2}$$\end{document} . The outer iteration, i.e., the SQP method for the Störmer–Verlet scheme or Newton method for the backward Euler scheme, stops after the relative error with the solution from the previous iteration reaches \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {tol}}_{{\text {SQP}}}=10^{-5}$$\end{document} in all variables.

The first iteration of the SQP/Newton method is commenced using suitable initial guesses. In the case of the Störmer–Verlet scheme, we first solve the problem on the coarsest mesh where our initial guesses for the state variables are taken to be the desired states, while we take vectors of zeros as the initial guesses for the adjoint variables. After obtaining the solution, we proceed with solving on a finer mesh, using the previously calculated solution interpolated onto the current mesh and time interval and multiplied by 0.8 as an initial guess. This multiplication is to ensure that we are able to retrieve an accurate solution, without initializing in the immediate neighbourhood of the solution. We repeat this process for finer meshes. To accommodate the inbuilt MINRES solver, we use slightly different initial guesses for the backward Euler scheme on the coarsest grid level: specifically, we use the target states multiplied by 0.4 as initial guesses for u, v. For finer meshes, we proceed in the same way as outlined for the Störmer–Verlet scheme. We note that using this approach for the Störmer–Verlet scheme leads to the same numerical solution.

The results in Tables 1, 2, 3 and 4 were computed on a CentOS Linux HPC facility with Intel(R) Xeon(R) Gold 6348 CPU @ 2.60GHz. The meshes are defined by the spacing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h=\frac{2^{1-i}}{10}$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i\in \{1,2,3,4\}$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau =\frac{h}{5}$$\end{document} for the Störmer–Verlet scheme and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau =2h^2$$\end{document} for the backward Euler scheme, with the total number of time-steps \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_t^S$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_t^E$$\end{document} , respectively. The number of degrees of freedom (DoF) is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4N_t^SN_x$$\end{document} for the Störmer–Verlet and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$4(N_t^E-1)N_x$$\end{document} for the backward Euler scheme. We state the weighted errors with the exact solution for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma =2$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta \in \{10^{-2},10^{-3}\}$$\end{document} . We report the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell ^2$$\end{document} in space and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell ^\infty $$\end{document} in time error, given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} u_{{\text {error}}} := \max _{ 0 \le i \le N_t} (h |\boldsymbol{u}^i - u^{{\text {sol}}}(i\tau ,{\textbf {x}}) |_2 ), \end{aligned}$$\end{document}

u̧sing values of the exact solution on nodal points, which is used for all the variables. The CPU time measurements only include the SQP/Newton method and the calculation of the system that is solved with MINRES, i.e., we do not include the operations preceding the linear solve nor the calculation of the variables that are excluded from the large system, i.e., \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{p}^0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{q}^0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{u}^{N_t}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\boldsymbol{v}^{N_t}$$\end{document} , in the time measurement.

Comparing the accuracy between the two schemes in Tables 1–4, we notice that for the same mesh size, both schemes achieve errors of roughly the same order. Comparing the errors across the two values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} considered, our strategy generally outperforms the backward Euler method. Note in particular that for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-3}$$\end{document} at the finest mesh, the Störmer–Verlet scheme achieves smaller errors than the backward Euler scheme despite having eight times larger time-step size.

While the number of SQP iterations in Tables 1 and 2 is fairly consistent across the different meshes, the number of MINRES iterations increases about three-fold for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-2~}$$\end{document} and two-fold for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-3}$$\end{document} at the finest mesh level. This increase might be due to the differing structure of the system matrix resulting from our strategy applied to the PDECO problem, as opposed to the backward Euler scheme, and the different relations between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document} and h used to test each method, combined with the large dimension of the problem and the features of the nonlinearity of the problem affecting the linear algebraic properties of the system at the finest mesh level. Despite this, we can noticeably reduce the CPU time for both values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} , while achieving similar accuracy as with the backward Euler scheme. Therefore, using our strategy for solving the given PDECO problem at a fine resolution allows us to save computational resources since we can reduce the problem dimension by taking larger time-steps while maintaining good accuracy.

Data-driven problem

Table 5. Parameter identification problem: squared \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q)$$\end{document} norms of individual terms in cost functional, average number of MINRES iterations, and total SQP iterations, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-2}$$\end{document} i

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [u-\hat{u}]^2$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [v-\hat{v}]^2 $$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [a]^2$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [b]^2$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{{\text {it}}}_M$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {it}}_{{\text {SQP}}}$$\end{document} 22.89e-55.74e-50.17860.191114432.85e-54.16e-50.17890.187713442.89e-53.75e-50.17910.1890128 Table 6. Parameter identification problem: squared \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q)$$\end{document} norms of individual terms in cost functional, average number of MINRES iterations, and total SQP iterations, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-3}$$\end{document} i

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [u-\hat{u}]^2$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [v-\hat{v}]^2 $$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [a] ^2$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [b] ^2$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{{\text {it}}}_M$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {it}}_{{\text {SQP}}}$$\end{document} 23.03e-72.25e-60.18080.205711433.31e-71.04e-60.18100.199010444.72e-77.96e-70.18060.199796 Table 7. Parameter identification problem: squared \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q)$$\end{document} norms of individual terms in cost functional, average number of MINRES iterations, and total SQP iterations, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-4}$$\end{document} i

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [u-\hat{u}]^2$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [v-\hat{v}]^2 $$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [a] ^2$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q) [b] ^2$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{{\text {it}}}_M$$\end{document}

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {it}}_{{\text {SQP}}}$$\end{document} 23.42e-98.84e-80.18110.21108434.35e-99.00e-80.18120.20158446.30e-88.39e-80.18130.201464

Fig. 1. Snapshot of the solution at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.7$$\end{document} for the desired states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{u}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{v}$$\end{document} (left) and computed state variables u, v (middle), and at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.695$$\end{document} for the computed control variables a, b (right), with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-2}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau = 10^{-2}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h = 0.0125$$\end{document} . The colorbar takes account of values from the whole time interval Fig. 2. Snapshot of the solution at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=2$$\end{document} for the desired states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{u}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{v}$$\end{document} (left) and computed state variables u, v (middle), and at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=1.995$$\end{document} for the computed control variables a, b (right), with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-2}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau = 10^{-2}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h = 0.0125$$\end{document} . The colorbar takes account of values from the whole time interval Fig. 3. Evolution of the mean of the numerical solution for the control variables a, b across the time interval compared to those of the control variables \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_G$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_G$$\end{document} used to generate the target state at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=5$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta =10^{-2}$$\end{document}

For testing purposes, we generate the target states \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{u}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{v}$$\end{document} directly using the Schnakenberg model implemented in FEniCS. For the Turing patterns to emerge, we introduce small perturbations to the homogeneous steady states as in [13], with values for the control parameters \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_G=0.126779$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_G = 0.792366$$\end{document} taken to be constant across the whole domain. We take snapshots of the patterns once they are stationary, at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=5$$\end{document} . Since in our PDECO problem we assume the target states evolve in time, for simplicity we linearly interpolate from target states that are zero at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0$$\end{document} to acquire the snapshots at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=T$$\end{document} . Here we use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=2$$\end{document} , the fixed time-step size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau = 10^{-2}$$\end{document} , and tolerances \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {tol}}_{{\text {M}}} = 10^{-7}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\text {tol}}_{{\text {SQP}}} = 10^{-6}$$\end{document} .

The indices i correspond to mesh sizes h where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h = \frac{2^{1-i}}{10}$$\end{document} . The squared \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(Q)$$\end{document} -norms of the misfits between the desired and computed states remain very small and similar for different mesh sizes across different values of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} . The norms of the control parameters also remain similar in each table, which suggests consistency of the results regardless of the mesh resolution. The average number of MINRES iterations decreases for finer meshes; however, as the size of the problem increases, more SQP iterations are needed (See Tables 5, 6 and 7), which differs from the behaviour observed in the previous section when the time-step size was varied.

The outcome of the numerical simulations can be observed in Figures 1 and 2 with snapshots of the desired states, computed states and computed controls at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=0.7$$\end{document} and the final time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=2$$\end{document} , respectively. We recall that the state variables are evaluated at full time steps, while the controls are calculated via the adjoint variables evaluated at half time-steps as derived in Sections 4 and 5. Both the desired and computed states in each snapshot are in good agreement which is confirmed by the norms of the differences \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u - \hat{u}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v - \hat{v}$$\end{document} in Table 5. We note that since no precautions were taken for positivity preservation, the control variable b has reached negative values within the time interval. Importantly, it is obvious that to drive the state variables to the desired pattern, the values of the control variables cannot be homogeneous across the spatial domain. Therefore we compare the mean values of the control variables across the domain with the constants \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_G$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b_G$$\end{document} , used to generate the images at the final time. In Figure 3, we observe that at the final time the means of the computed controls do indeed roughly achieve the desired control values. Since the desired state to be achieved at final time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=2$$\end{document} was generated using the state equations at time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=5$$\end{document} , we can also conclude that we were able to drive the state to the desired state in a shorter time. This means that even though in reality we may not know the actual duration of the evolution into a particular pattern for a biological process, we can still obtain a good approximation for the unknown parameters.

Conclusions

In this paper, we have devised an SQP-based method for challenging and large-scale parameter identification problems constrained by reaction–diffusion equations, based on the Schnakenberg model for pattern formation in mathematical biology. Specifically, we have devised a paradigm for time-stepping and approximation of cost functional terms, used in conjunction with a suitable resolution in the spatial variables, such that optimization and discretization steps commute within our numerical algorithm. Further, in computational experiments our method demonstrates second-order convergence upon mesh refinement in time and space. We constructed a feasible and efficient preconditioned iterative solver for the linear systems arising from each iteration of the SQP solver. Problems with dimensions of the order of tens of millions have been successfully tackled using our solution strategy. Future work includes devising a theoretical stability and convergence analysis for our numerical scheme and applying it to even more demanding problems from mathematical biology.