# Precision and Adaptability of YOLOv5 and YOLOv8 in Dynamic Robotic
Environments

1st Victor A. Kich  Intelligent Robot Laboratory
University of Tsukuba
Tsukuba, Japan
<EMAIL_ADDRESS>2nd Muhammad A. Muttaqien  Intelligent Robot
Laboratory
University of Tsukuba
Tsukuba, Japan
<EMAIL_ADDRESS>3rd Junya Toyama  Intelligent Robot Laboratory
University of Tsukuba
Tsukuba, Japan
<EMAIL_ADDRESS>4rd Ryutaro Miyoshi  Intelligent Robot
Laboratory
University of Tsukuba
Tsukuba, Japan
<EMAIL_ADDRESS>5th Yosuke Ida  Intelligent Robot
Laboratory
University of Tsukuba
Tsukuba, Japan
<EMAIL_ADDRESS>6th Akihisa Ohya  Intelligent Robot
Laboratory
University of Tsukuba
Tsukuba, Japan
<EMAIL_ADDRESS>7th Hisashi Date  Intelligent Robot Laboratory
University of Tsukuba
Tsukuba, Japan
<EMAIL_ADDRESS>
###### Abstract

Recent advancements in real-time object detection frameworks have spurred
extensive research into their application in robotic systems. This study
provides a comparative analysis of YOLOv5 and YOLOv8 models, challenging the
prevailing assumption of the latter’s superiority in performance metrics.
Contrary to initial expectations, YOLOv5 models demonstrated comparable, and
in some cases superior, precision in object detection tasks. Our analysis
delves into the underlying factors contributing to these findings, examining
aspects such as model architecture complexity, training dataset variances, and
real-world applicability. Through rigorous testing and an ablation study, we
present a nuanced understanding of each model’s capabilities, offering
insights into the selection and optimization of object detection frameworks
for robotic applications. Implications of this research extend to the design
of more efficient and contextually adaptive systems, emphasizing the necessity
for a holistic approach to evaluating model performance.

###### Index Terms:

YOLO, Robotic Mapping, Object Detection, Real-Time Processing

## Suplementary Material

The code and pre-trained models used in this research are available in the
link: https://github.com/victorkich/BGBox-TsukubaChallenge. The dataset for
box classification can be found in: https://universe.roboflow.com/victor-
augusto-kich/bluebox-and-greenbox. While the dataset for the letter
classification can be found by the link:
https://universe.roboflow.com/tsukuba-challenge/letter-classification.

##  I Introduction

The Tsukuba Challenge, an annual event since 2007 in Tsukuba, Ibaraki, Japan,
is a technology challenge that has been held every year. Its objective is to
enhance real-world autonomous driving technology, operating in the same
environments as people. Participants from various organizations strive to
advance technology by developing mobile robots, experimenting, and sharing
outcomes.

The utility of robotic systems is evident across various domains such as
autonomous navigation, surveillance, and environmental monitoring. Among
these, terrestrial robots play a crucial role in navigating and understanding
complex environments. The Challenge serves as a testament to advancements in
real-world autonomous navigation, focusing on the seamless integration of
mobile robots within human environments. This challenge underscores the
importance of robust object detection capabilities, which are addressed in
this study.

In response to the 2023 Tsukuba Challenge 111Tsukuba Challenge website:
https://tsukubachallenge.jp, which emphasizes the detection and classification
of specific objects like traffic signs and boxes, this study contributes by
offering a comparative analysis of several YOLO models, including the latest
iterations like YOLOv8 [1]. Despite YOLOv8’s advancements, such as improved
feature fusion and spatial attention, this paper maintains a critical
perspective on each model’s suitability for robotic applications, particularly
in terms of real-time processing capabilities and adaptability to diverse
settings [2, 3, 4].

Figure 1: Image samples from the proposed dataset designed to train various
YOLO models for real-time object detection

This work seeks to evaluate these models comprehensively across a range of
scenarios pertinent to terrestrial robots, thereby assessing their performance
in terms of detection accuracy, processing efficiency, and practical
deployment in outdoor navigation tasks. The contributions of this study are
listed as follows:

  1. 1.

A broad evaluation of YOLO models, highlighting the practical integration of
these frameworks in robotic mapping to tackle real-world environmental
complexities.

  2. 2.

Development of a novel 1700-image dataset tailored to the specific
requirements of the proposed task, addressing the challenges outlined in the
study. Various data augmentation techniques were applied to enhance the
dataset’s diversity and robustness.

  3. 3.

Detailed performance comparisons across various YOLO versions, with a focus on
enhancing the robustness and efficiency of object detection in dynamic
settings.

By examining the intersection of advanced object detection technologies and
their application in robotic systems, this paper aspires to contribute to the
broader dialogue on enhancing the operational capabilities of autonomous
robots in diverse environments.

##  II Related Works

The domain of object detection has witnessed a remarkable evolution,
predominantly influenced by the advancements in deep learning architectures.
Seminal works such as those by Redmon et al. [5], who introduced the concept
of YOLO, revolutionized real-time object detection by framing it as a
regression problem, enabling the simultaneous prediction of object classes and
bounding boxes. Successive iterations of YOLO, including YOLOv4 presented by
Bochkovskiy et al. [6], have progressively enhanced the balance between
detection speed and accuracy, setting new benchmarks in the field. These
models have been pivotal in robotic applications where rapid environmental
perception is critical.

Recent advancements in object detection have been instrumental in addressing a
wide range of applications and challenges. A notable comparative study by
Agrawal et al. [7] delved into various recognition methods, including YOLO,
Faster R-CNN, and R-CNN, with a particular emphasis on real-time processing.
They evaluated key metrics such as mean Average Precision (mAP) and Frames Per
Second (FPS), highlighting YOLO’s effectiveness in object localization and
tracking, a critical factor in high-stakes scenarios requiring precise object
detection. In service automation, Ge et al. [8] introduced YOLO-GG, a
customized model for empty-dish recycling robots in the food service sector.
Designed for embedded systems, YOLO-GG achieved an impressive 99.34% mAP and
71.2 FPS, showcasing YOLO’s practicality and efficiency in real-world settings
and its capability to meet specific industry needs.

Further to these foundational models, more recent iterations such as YOLOv5 by
Jocher et al. [9], and the advancements encapsulated in YOLOv8, have pushed
the frontiers of object detection performance. Comparative analyses, such as
the work by Tan et al. [10], often emphasize key metrics like mean Average
Precision (mAP) and Frames Per Second (FPS), to gauge the improvements offered
by these newer models. Notably, YOLOv5 has emerged as a robust framework that
not only accelerates inference times but also offers a level of precision that
challenges its successors. This has led to a rich body of literature that not
only compares these models but also scrutinizes their performance across a
spectrum of real-world applications.

Figure 2: The vision system of the Kerberos Robot utilizes three cameras to
identify and categorize the designated target box. The procedure includes
cropping the area of interest, implementing a mask over the captured box to
isolate it, and finally, engaging a specialized neural network for letter
classification. This streamlined workflow ensures precise detection and
analysis.

Furthermore, in the realm of robotics, Putra et al. [11] successfully
implemented YOLOv4-tiny for the Indonesian Search and Rescue Robot Contest,
achieving an impressive 98.25% mAP and 17.8 FPS even with limited hardware
capabilities. This underscores YOLO’s ability to adapt to varying light
conditions and distances, solidifying its suitability for dynamic search and
rescue operations. In industrial contexts, Zuo et al. [12] demonstrated YOLO’s
versatility by employing it to detect weld surface defects, achieving a
significant 75.5% mAP and highlighting its adaptability across various
industrial scenarios. Additionally, Khokhlov et al. [13] further expanded
YOLO’s application by enhancing Tiny-YOLO object detection with geometrical
data, showcasing the continuous evolution and wide-ranging utility of YOLO-
based systems in diverse fields.

Despite the comprehensive benchmarks available, there remains a gap in
literature regarding the in-depth analysis of why certain models perform
better under specific conditions, particularly in robotics where environmental
variables are unpredictable. Ablation studies, such as those by Howard et al.
[14], provide insight into the contribution of individual model components to
overall performance. However, such studies are seldom applied to comparative
analyses between different versions of object detection frameworks. Our study
aims to fill this void by dissecting the performance of YOLOv5 and YOLOv8
within the context of robotic applications, offering a granular view of their
operational strengths and weaknesses.

##  III Methodology

Our study adopted an empirical approach to scrutinize the performance nuances
of the YOLOv5 and YOLOv8 frameworks. The methodology was twofold: firstly, to
implement and train the models using a standardized dataset designed for
robotic object detection; secondly, to analyze the trained models’ performance
through a series of metrics, primarily focusing on mAP and precision-recall
values. A brief overview of our approach to solve the proposed task can be
seen in the Fig.2

The experimentation involved the training of YOLOv5 and YOLOv8 models on a
curated dataset comprising images annotated for robotic vision tasks. Training
was executed over 200 epochs with a batch size of 16, using a constant
learning rate of 0.01 and employing early stopping based on a patience
parameter set to 20 epochs. Both, YOLOv8 and YOLOv5 models were fine-tuned
from weights pre-trained on the ImageNet dataset [15].

Figure 3: The overview of YOLOv8 architecture.

###  III-A Optimization Techniques

The implementation of class balancing has evidently supported the models in
achieving uniformity in detection capabilities across different classes. This
is particularly visible in the balanced precision and recall rates across
models, ensuring that no single class dominates the detection process, which
is vital for the varied requirements of robotic applications.

By focusing on challenging samples where misclassification was frequent, such
as differentiating between blue boxes and the blue sky, hard negative mining
significantly reduced false positives. This is critical in robotics, where
distinguishing between similar objects can mean the difference between correct
navigation and a potential collision.

The improvement in mAP, especially under varied lighting and environmental
conditions, underscores the effectiveness of domain adaptation. This technique
has enhanced the robustness of the models, allowing them to perform reliably
in different settings—an essential attribute for robots operating in real-
world conditions.

The application of advanced post-processing techniques like non-maximum
suppression (NMS) and varying IoU thresholds has refined the detection
precision. This is evidenced by the sustained high performance in precision
metrics, which is crucial for ensuring that robots make decisions based on
reliable data.

Adjusting the depth and width of the models to fit specific use cases has
shown to optimize the computational efficiency without compromising the
accuracy. Smaller models have demonstrated surprisingly robust performance,
particularly in environments with limited computational resources, which is a
common scenario in mobile robotics.

###  III-B Data Preparation and Model Training

The dataset used for training was assembled to represent the variety of
scenarios a robotic system would encounter. It included a diversity of
lighting conditions, occlusions, and angles to challenge the models’
generalizability. The images were evenly distributed across the classes to
ensure class balance and were augmented to simulate the conditions of a
robot’s operational environment. Both YOLOv8 and YOLOv5 models underwent
rigorous training, with specific attention to the fine-tuning of
hyperparameters that could influence their real-world performance.

###  III-C Ablation Insights

An ablation study was conducted (Section IV) to understand the influence of
various factors on the models’ performance. We examined the impact of class
balancing, hard negative mining, domain adaptation, post-processing
techniques, and model size adjustment on the precision and recall values.

Initial findings suggested that while YOLOv8 models benefit from architectural
optimizations (Fig. 3), these do not necessarily translate to higher precision
in complex, real-world scenarios often encountered in robotic applications.
Conversely, the matured training procedures and optimizations of YOLOv5 models
could be contributing factors to their robust precision performance. Detailed
analysis into the ablation results is presented in the subsequent section,
shedding light on the intricate dynamics between model architecture, training
paradigms, and dataset characteristics.

##  IV Ablation Study

The ablation study conducted provides a deep insight into the performance
enhancements and limitations of the YOLOv8 and YOLOv5 models through various
optimization techniques including class balancing, hard negative mining,
domain adaptation, post-processing techniques, and model size adjustments.
This study critically evaluates how each factor contributes to the overall
performance in robotic applications.

###  IV-A Performance Evaluation

Performance evaluation focused on mAP across two different IoU thresholds:
0.5, a standard for object detection, and 0.5:0.95, a more stringent metric
that considers a range of IoU thresholds. Precision and recall rates were
calculated to reflect the models’ accuracy and completeness in identifying
objects.

Table I: Performance metrics across all models. Architecture | mAP (0.5) | mAP (0.5:0.95) | Precision | Recall
---|---|---|---|---
YOLOv8m | 0.964 | 0.821 | 0.885 | 0.878
YOLOv8l | 0.966 | 0.826 | 0.886 | 0.878
YOLOv8x | 0.951 | 0.772 | 0.865 | 0.866
YOLOv5m | 0.978 | 0.829 | 0.918 | 0.917
YOLOv5l | 0.976 | 0.768 | 0.867 | 0.852
YOLOv5x | 0.983 | 0.834 | 0.917 | 0.913

As Table I indicates, YOLOv5 models demonstrated a competitive edge in
precision, particularly the ‘m’ and ‘x’ configurations, despite the
anticipated superiority of YOLOv8 in terms of architectural advancements. This
unexpected turn of results initiated a deeper dive into the specific aspects
of model training and data handling that may have contributed to this outcome.

###  IV-B Analysis of Performance Metrics

(a) mAP increase.

(b) mAP 0.5:0.95 increase.

(c) Class loss decrease.

(d) YOLOv8 object loss decrease.

(e) YOLOv5 object loss decrease.

Figure 4: Performance metrics for YOLOv8 and YOLOv5 during training and
validation over several epochs. The legends for each graph include variations
of the model: ‘m’, ‘l’, and ‘x’, representing different model sizes or
configurations, with ‘train’ indicating training data and ‘val’ indicating
validation data.

As illustrated in Figures 4(a) and 4(b), the Mean Average Precision (mAP) for
both YOLOv8 and YOLOv5 models showed significant improvements over epochs,
with YOLOv5 models generally maintaining a slight edge over YOLOv8,
particularly in the more stringent mAP 0.5-0.95 metric. This could be
attributed to the more mature training processes and architectural
optimizations that are specifically tuned in YOLOv5 models for handling
diverse and complex object classes.

The class loss and object loss curves (Figures 4(c), 4(d), and 4(e)])
demonstrate a consistent decrease, indicating effective learning and
adaptation by both model architectures. However, YOLOv5 models exhibit a more
rapid and stable decline in object loss, suggesting better generalization
across varied object detections which is crucial for dynamic environments in
robotic navigation.

###  IV-C Implications for Robotic Applications

The outcomes of this ablation study illustrate that while newer model
architectures like YOLOv8 offer significant advancements, the specific tuning
and optimizations inherent in YOLOv5 models provide a competitive edge in
terms of precision and reliability in object detection tasks. These insights
direct a clear path towards optimizing deep learning models for practical and
efficient use in robotics, where accuracy, speed, and adaptability to diverse
environments are paramount.

The refined understanding of how different optimizations impact model
performance not only aids in selecting appropriate models for specific robotic
tasks but also guides future developments in robotic vision systems. The
results encourage a more nuanced approach to model selection and training,
advocating for a balanced consideration of various factors beyond mere
architectural advancements.

###  IV-D Insights from Comparative Performance

The higher precision of YOLOv5 models in certain configurations suggests that
optimizations specific to YOLOv5 may be better aligned with the
characteristics of the dataset used. These could include better feature
representation at different scales or more effective background class
suppression. Additionally, the refined pre-training on ImageNet might have
provided YOLOv5 models with a more robust feature set for the given detection
tasks, highlighting the importance of the pre-training dataset and procedure.

The convergence of YOLOv5’s architectural maturity with a well-prepared
training dataset challenges the presupposition that newer models naturally
yield better performance. It suggests that the interplay between dataset
intricacies and model architectures plays a more pivotal role than previously
acknowledged. Our study suggests a paradigm shift in model selection criteria
for robotic applications, with an emphasis on real-world testing and
performance analysis over hypothetical model superiority.

##  V Conclusion

This study has critically evaluated the performance of YOLOv5 and YOLOv8
models in robotic applications, with a particular focus on the nuanced
dynamics of advanced object detection frameworks. Through an extensive
ablation study and a rigorous analysis of optimization techniques such as
class balancing, hard negative mining, domain adaptation, post-processing
techniques, and model size adjustment, it has been demonstrated that while
newer models like YOLOv8 offer significant technological advancements, they do
not universally outperform older versions across all metrics. In fact, YOLOv5
models frequently exhibited superior precision and reliability in detection
tasks, underscoring the importance of model selection based on specific
operational contexts rather than solely on model novelty.

The findings of this research not only highlight the critical role of tailored
model training and optimization for robotic applications but also suggest a
paradigm shift towards more application-focused approaches in the deployment
of object detection models. For future work, we aim to explore the integration
of these findings into real-world robotic systems, further testing the models
in dynamic environments to refine their efficiency and accuracy. This would
not only help in advancing the state-of-the-art in robotic navigation and
interaction but also in ensuring that the benefits of technological
advancements are fully leveraged to enhance the operational capabilities of
autonomous systems in complex and varied settings.

## References

  * [1] Z. E. Babi, N. A. Khomami, M. T. Masouleh, and A. Kalhor, “Autonomous robotic assembly and sequence planning based on yolov8,” in _2023 11th RSI International Conference on Robotics and Mechatronics (ICRoM)_. IEEE, 2023, pp. 840–846.
  * [2] Z. Diao, P. Guo, B. Zhang, D. Zhang, J. Yan, Z. He, S. Zhao, C. Zhao, and J. Zhang, “Navigation line extraction algorithm for corn spraying robot based on improved yolov8s network,” _Computers and Electronics in Agriculture_ , vol. 212, p. 108049, 2023.
  * [3] T. Wu and Y. Dong, “Yolo-se: Improved yolov8 for remote sensing object detection and recognition,” _Applied Sciences_ , vol. 13, no. 24, p. 12977, 2023.
  * [4] M. Mahasin and I. A. Dewi, “Comparison of cspdarknet53, cspresnext-50, and efficientnet-b0 backbones on yolo v4 as object detector,” _International Journal of Engineering, Science and Information Technology_ , vol. 2, no. 3, pp. 64–72, 2022.
  * [5] J. Redmon and A. Farhadi, “Yolo9000: better, faster, stronger,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 7263–7271.
  * [6] A. B. et al, “Yolov4: Optimal speed and accuracy of object detection,” _arXiv preprint arXiv:2004.10934_ , 2020.
  * [7] P. Agrawal, G. Jain, S. Shukla, S. Gupta, D. Kothari, R. Jain, and N. Malviya, “Yolo algorithm implementation for real time object detection and tracking,” in _IEEE Students Conference on Engineering and Systems (SCES)_ , 2022.
  * [8] Y. Ge, X. Yue, and L. Meng, “Yolo-gg: A slight object detection model for empty-dish recycling robot,” in _International Conference on Advanced Mechatronic Systems (ICAMechS)_ , 2022.
  * [9] G. Jocher, A. Stoken, J. Borovec, L. Changyu, A. Hogan, L. Diaconu, F. Ingham, J. Poznanski, J. Fang, L. Yu _et al._ , “ultralytics/yolov5: v3. 1-bug fixes and performance improvements,” _Zenodo_ , 2020.
  * [10] M. Tan, R. Pang, and Q. V. Le, “Efficientdet: Scalable and efficient object detection,” in _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , 2020, pp. 10 781–10 790.
  * [11] Y. I. Putra, A. H. Alasiry, A. Darmawan, H. Oktavianto, and Z. M. E. Darmawan, “Camera-based object detection and identification using yolo method for indonesian search and rescue robot competition,” in _5th International Seminar on Research of Information Technology and Intelligent Systems (ISRITI)_ , 2022.
  * [12] Y. Zuo, J. Wang, and J. Song, “Application of yolo object detection network in weld surface defect detection,” in _IEEE 11th Annual International Conference on CYBER Technology in Automation, Control, and Intelligent Systems (CYBER)_ , 2021.
  * [13] I. Khokhlov, E. Davydenko, I. Osokin, I. Ryakin, A. Babaev, V. Litvinenko, and R. Gorbachev, “Tiny-yolo object detection supplemented with geometrical data,” in _IEEE 91st Vehicular Technology Conference (VTC2020-Spring)_ , 2020.
  * [14] A. Howard, M. Sandler, G. Chu, L.-C. Chen, B. Chen, M. Tan, W. Wang, Y. Zhu, R. Pang, V. Vasudevan _et al._ , “Searching for mobilenetv3,” in _Proceedings of the IEEE/CVF international conference on computer vision_ , 2019, pp. 1314–1324.
  * [15] J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei, “Imagenet: A large-scale hierarchical image database,” in _2009 IEEE conference on computer vision and pattern recognition_. Ieee, 2009, pp. 248–255.


11institutetext: J. H. Chaudhry 22institutetext: T. Giannini 33institutetext:
Z. Stevens 44institutetext: The University of New Mexico, Albuquerque, NM
87131, USA. 44email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>55institutetext: S. J. Tavener 66institutetext: Colorado State University,
Fort Collins, CO 80523, USA. 66email<EMAIL_ADDRESS>77institutetext: D. Estep 88institutetext: Simon Fraser University, Burnaby,
BC V5A 1S6, Canada. 88email<EMAIL_ADDRESS>
# Error estimation for the time to a threshold value in evolutionary partial
differential equations

Jehanzeb H. Chaudhry  Donald Estep  Trevor Giannini  Zachary Stevens  Simon J.
Tavener

###### Abstract

We develop an a posteriori error analysis for a numerical estimate of the time
at which a functional of the solution to a partial differential equation (PDE)
first achieves a threshold value on a given time interval. This quantity of
interest (QoI) differs from classical QoIs which are modeled as bounded linear
(or nonlinear) functionals of the solution. Taylor’s theorem and an adjoint-
based a posteriori analysis is used to derive computable and accurate error
estimates in the case of semi-linear parabolic and hyperbolic PDEs. The
accuracy of the error estimates is demonstrated through numerical solutions of
the one-dimensional heat equation and linearized shallow water equations
(SWE), representing parabolic and hyperbolic cases, respectively.

##  1 Introduction

Numerical solutions of parabolic and hyperbolic differential equations are
essential to the study of physical phenomena that evolve over space and time.
To meet the requirements of uncertainty quantification, adjoint-based methods
for estimating the error in a quantity of interest that can be expressed as a
linear or nonlinear functional of the solution are well developed. However,
the _time_ at which a functional of the solution achieves a threshold value,
which we refer to as the time to an “event”, is often the primary problem for
a computational study of a physical system. Examples include: the time at
which the concentration of a chemical species reaches a critical value at a
particular location, the time at which the temperature at a particular
location drops below a critical value, and the time at which a traveling wave
reaches a certain location given specific initial values and topography.
Unfortunately, the error in an estimated time to a particular event cannot be
quantified using standard approaches.

Consider a semi-linear evolutionary PDE of the form

$u_{t}(x,t)+Lu(x,t)=f(u,x,t),\quad x\in\Omega,\;t\in(0,T],$ (1)

with appropriate initial and boundary conditions. Here $u_{t}=\frac{\partial
u}{\partial t}$, $L$ is a differential operator in $x$ that is linear in $u$,
and $f(u,x,t)$ is a differentiable function. While $u,L$ and $f$ may be scalar
or vector valued, these cases are not distinguished in our analyses since this
will be obvious from the context.

Classical a posteriori error analysis considers quantities of interest (QoIs)
that can be expressed as bounded functionals of the solution and utilizes
variational analysis, generalized Green’s functions and computable residuals
Est95 ; AinOde00 ; EstLarWil00 ; BecRan01 ; GilSul02 ; BanRan03 ; Bar04 ;
EstHolMik02 ; CaoPet04 ; houston2017adjoint . Recent extensions to multiscale
and multiphysics systems include CarEstTav08 ; Est09 ; EstGinTav12 ;
JohChaCarEstGinLarTav15 ; Cha18 . Finite difference and volume methods, and
explicit time integration schemes may be analyzed by reformulating these
schemes as equivalent finite element methods, see e.g., ColEstTav14 ;
ChaEstGinShaTav15 ; ChaEstGinTav15 ; ColEstTav15 ; ChaEstTavCarSan16 ;
ChaColSha17 ; ChaShaWil19 . Nonlinear QoIs are typically handled by
linearization around a computed solution, e.g., BecRan01 . The key point
motivating this paper is that the _time_ a functional of the solution of a PDE
first achieves a threshold value cannot be expressed as a bounded linear
functional of $u$, nor can it be trivially linearized.

Let $G(u;t)$ be a linear functional of $u(x,t)$ which is implicitly dependent
on $t$ through $u$, and let $R$ be a threshold value of $G(u;t)$. Assume there
are one or more times $t^{\star}$ during the interval $(0,T]$ for which
$G(u;t^{\star})=R$ is satisfied and define the time $H(u,\tau)$ as

$H(u,\tau)=\min_{t\in(\tau,T]}\underset{t}{\arg}(G(u;t)=R).$ (2)

Here $\tau$ is specified in order to obtain different occurrences of the event
$G(u;t)=R$. For example, $H(u,0)$ is the time on the interval (0,T] of the
_first_ occurrence of the event of the $G(u;t)=R$. Similarly, the time the
_second_ occurrence of the event $G(u;t)=R$ is $H(u,\tau)$, where
$H(u,0)\leq\tau<H(u,H(u,0))$, i.e., for any $\tau$ between the times of first
and second occurrence of the event. The quantity of interest $Q(u)$ for a
fixed $\tau$ is defined as

$Q(u)=H(u,\tau).$ (3)

In recent work ChaEstSteTav21 , an a posteriori error analysis was performed
in order to estimate the error in the time to an event (the time of the first
occurrence of a functional achieving a threshold value), in the context of
ordinary differential equations (ODEs). In the current work, the analysis of
ChaEstSteTav21  is extended to evolutionary PDEs and to the QoI specified by
(3). An error estimate is derived using Taylor series linearized around the
approximate (computed) value of the QoI. The error estimate, which contains
terms involving the unknown error of the solution to the PDE, is made
computable through the use of solutions to certain adjoint problems. In
comparison with the earlier work of ChaEstSteTav21 , the current analysis
requires the solution of an additional adjoint problem arising due to the
presence of the spatial operator $L$. We apply the analysis to the one-
dimensional heat equation and to the one-dimensional linearized shallow water
equations (SWE).

Numerical schemes to solve the semi-linear evolutionary PDE (1) are developed
in §2. A variational formulation is introduced in §2.1 and a space-time
Galerkin finite element discretization in §2.2. The error in QoI (3) is
defined in §3 and the _a posteriori_ error analysis for this QoI is presented
in §3.1. A computable estimate is developed in §3.3 by taking advantage of
classical error estimates which are reviewed in §3.2. Error estimates for two
model problem are developed in §4, the one-dimensional heat equation in §4.1
and the one-dimensional linearized shallow water equations in §4.2,
respectively. Numerical results for both the linear and nonlinear heat
equations and for the linearized shallow water equations with a range of
different bottom topographies are presented in §5.

##  2 Numerical methods

Adjoint-based, a posteriori error estimates utilize the variational (or weak)
form of PDEs. A variational formulation for the semi-linear evolutionary PDE
(1) is introduced in §2.1 and a space-time Galerkin finite element
discretization is presented in §2.2.

###  2.1 Variational formulation

Let $W$ be a generic Hilbert space over the domain $\Omega$ with inner product
$\left(\cdot,\cdot\right)_{W}$ and the induced W-norm $\|\cdot\|_{W}$. Let
$L^{2}(0,T;W)$ denote the space of functions $w$ whose $W$-norm is square
integrable $\forall t\in[0,T]$. The variational form of (1) is: Find $u\in
L^{2}(0,T;W)$ such that

$\left(u_{t}(\cdot,t),v(\cdot,t)\right)+\left(L_{1}u(\cdot,t),L_{2}v(\cdot,t)\right)=\left(f(u,\cdot,t),v(\cdot,t)\right),\quad\forall
v\in W,\quad\forall t\in(0,T],$ (4)

where $f(u,\cdot,t)\in W$ for any $t\in(0,T]$ and
$\left(a,b\right)=\int_{\Omega}a\cdot b\;{\rm d}\Omega$ is the spatial $L^{2}$
inner product over $\Omega$. The two operators $L_{1},L_{2}$ are linear
differential operators such that

$L_{2}^{*}L_{1}u=Lu.$ (5)

The operator $L_{2}^{*}$ is the formal adjoint of $L_{2}$ which satisfies the
property that

$\left(u,L_{2}v\right)=\left(L_{2}^{*}u,v\right),\qquad\forall u\in W,v\in W.$
(6)

Under suitable conditions on the operators $L_{1}$ and $L_{2}$ (and hence
$L$), (4) is well-posed; the reader is referred to evans10  for more details.
Although the differential operator $L_{1}$ may be of lower order than $L$, if
a solution $u$ to (4) is sufficiently differentiable, it also satisfies (1).

###  2.2 Space-time Galerkin finite element method

Let $\mathcal{T}_{h}$ be a shape-regular triangulation of the spatial domain
$\Omega$ where $h$ denotes the diameter of the largest element. For the
numerical examples in §5, $\Omega$ is one-dimensional and $\mathcal{T}_{h}$ is
a uniform set of $N_{x}$ subintervals of the domain where $h$ denotes the
width of the subintervals. Let $V^{q_{s}}$ be the standard space of continuous
piecewise polynomials of degree $q_{s}$ on $\mathcal{T}_{h}$. This space may
be scalar or vector valued depending upon the context.

The time domain $[0,T]$ is discretized into $N_{t}$ uniform subintervals of
width $p$ denoted $\mathcal{T}_{p}$, i.e.,
$\mathcal{T}_{p}=\\{I_{0},\dots,I_{n},\dots,I_{N_{t}-1}\\}$ where
$I_{n}=[t_{n},t_{n+1}]$. Let $l_{j}(t)$ be the $j$th degree Lagrange
polynomial in time on the interval $I_{n}$. For the space-time slab
$S_{n}=\Omega\times I_{n}$ the solution space is defined as

$W_{n}^{q_{t},q_{s}}=\left\\{w(x,t)\,|\,w(x,t)={\sum_{j=0}^{q_{t}}l_{j}(t)v_{j}(x),\;v_{j}(x)}\in
V^{q_{s}},\,(x,t)\in S_{n}\right\\}.$ (7)

We define the continuous Galerkin method, cG($q_{t},q_{s}$) as: Find a
continuous in time $U\in W_{n}^{q_{t},q_{s}}$ for $n=0,\dots,N_{t}-1$, such
that

$\int_{t_{n}}^{t_{n+1}}\Big{[}\Big{(}U_{t}(\cdot,t),v(\cdot,t)\Big{)}+\Big{(}L_{1}U(\cdot,t),L_{2}v(\cdot,t)\Big{)}\Big{]}\;{\rm
d}t=\int_{t_{n}}^{t_{n+1}}\Big{(}f(\cdot,t),v(\cdot,t)\Big{)}\;{\rm d}t,\hskip
14.22636pt\forall v\in W_{n}^{q_{t}-1,q_{s}}.$ (8)

##  3 _A posteriori_ error estimation

Let $u(x,t)$ be the true solution of (1) and $G(u;t)$ be a time-dependent
linear functional of $u(x,t)$ that may be expressed in terms of the inner
product,

$G(u;t)=\int_{\Omega}w(x)\cdot u(x,t)\;{\rm
d}x=\left(w(\cdot),u(\cdot,t)\right),$ (9)

where $w(x)$ is some differentiable weight function. Assume that
$\left.w(x)\right|_{\partial\Omega}=0$. For a chosen threshold $R$ and fixed
value of $\tau$ in (2), denote the true value of the QoI

$t_{t}:=Q(u).$ (10)

Let $U(x,t)$ be an approximation of the solution satisfying (8) and

$t_{c}:=Q(U),$ (11)

the corresponding computed value of the QoI. Figure 1 illustrates these
definitions. We denote the error in the approximate _solution_ as
$e(x,t)=u(x,t)-U(x,t)$. As noted in §1, the a posteriori error analysis for
functionals of $e(x,t)$ is well established. However, we seek to estimate the
error

$e_{Q}=t_{t}-t_{c},$ (12)

and specifically, to calculate a _computable_ estimate $\nu\approx e_{Q}$. We
define the _effectivity ratio_ of this estimate to be

$\rho_{\rm eff}=\frac{\nu}{e_{Q}},$ (13)

where clearly the goal is to construct an error estimate for which $\rho_{\rm
eff}\approx 1$. We first derive an error representation in §3.1. Standard _a
posteriori techniques_ (see §3.2) are then employed to approximate the error
in three linear functionals of $e(x,t)$. The errors in these three functionals
are calculated as the inner product of solutions to three distinct adjoint
problems and computable residuals as described in §3.3.

(a)  $G(u;t)$ for $t\in[0,200]$

(b)  $G(u;t)$ and $G(U;t)$ for $t\in[89,91]$

Figure 1: (a) Reference values of the functional $G(u;t)$ for §5.4. The QoI is
the time for the third occurrence of the event $G(u;t)=R$. (b) Detail near the
time of third occurrence of the event $G(u;t)=R$ showing the true and
approximate values of the functional, $G(u;t)$ and $G(U;t)$ respectively, and
the corresponding true and approximate values of the QoI, $Q(u)$ and $Q(U)$.
The error is $e_{Q}=Q(u)-Q(U)$.

###  3.1 Error estimation for the time to an event $G(u;t)=R$

An adjoint-based _a posteriori_ error estimate for $e_{Q}$ as defined in (12)
is presented in Theorem 3.1 below.

###### Theorem 3.1

Consider the QoI defined by (3), where the functional $G$ is defined by (9),
and $t_{t}$, $t_{c}$, and the error in the QoI, $e_{Q}$, are defined by (10),
(11) and (12) respectively. Let

$\mathcal{D}(U,t_{c})=\left(L_{2}w(\cdot),L_{1}U(\cdot,t_{c})\right)-\left(w(\cdot),f(U,\cdot,t_{c})\right)-\left(\nabla_{u}f^{\top}(U,\cdot,t_{c})w(\cdot),e(\cdot,t_{c})\right)+\left(L_{2}w(\cdot),L_{1}e(\cdot,t_{c})\right),$
(14)

where $e(x,t)=u(x,t)-U(x,t)$, and $\nabla_{u}f$ is a square matrix with
entries $[\nabla_{u}f]_{ij}=\partial f_{i}/\partial u_{j}$. For semi-linear
time evolution PDEs (1), the error in the QoI is given by

$\displaystyle
e_{Q}=t_{t}-t_{c}=\frac{\left(w(\cdot),e(\cdot,t_{c})\right)+\left(w(\cdot),\mathcal{R}_{1}(\cdot,t_{c},t_{t})\right)}{\mathcal{D}(U,t_{c})-\left(w(\cdot),\mathcal{R}_{2}(u,U,\cdot,t_{c})\right)},$
(15)

where the two remainder terms are
$\mathcal{R}_{1}(x,t_{c},t_{t})=\mathcal{O}\left((t_{t}-t_{c})^{2}\right)$ and
$\mathcal{R}_{2}(u,U,x,t_{c})=\mathcal{O}\left((u-U)^{2}\right)$.

###### Proof

Using the definitions of $t_{t}$ and $t_{c}$,

$G(u(t_{t}))-G(U(t_{c}))=R-R=0.$ (16)

From (16) and (9)

$\left(w(\cdot),u(\cdot,t_{t})-U(\cdot,t_{c})\right)=0.$

Expanding the functional $G(u(t_{t}))$ about $t_{c}$ by expanding $u(x,t_{t})$
around $t_{c}$ using Taylor Series,

$\displaystyle 0$
$\displaystyle=\left(w(\cdot),u(\cdot,t_{c})+(t_{t}-t_{c})u_{t}(\cdot,t_{c})+\mathcal{R}_{1}(\cdot,t_{c},t_{t})-U(\cdot,t_{c})\right)$
(17)
$\displaystyle=\left(w(\cdot),e(\cdot,t_{c})\right)+(t_{t}-t_{c})\left(w(\cdot),u_{t}(\cdot,t_{c})\right)+\left(w(\cdot),\mathcal{R}_{1}(\cdot,t_{c},t_{t})\right),$

where $e(x,t)=u(x,t)-U(x,t)$ and the remainder
$\mathcal{R}_{1}(x,t_{c},t_{t})$ is
$\mathcal{O}\left((t_{t}-t_{c})^{2}\right)$. Rearranging (17),

$t_{t}-t_{c}=-\frac{\left(w(\cdot),e(\cdot,t_{c})\right)+\left(w(\cdot),\mathcal{R}_{1}(\cdot,t_{c},t_{t})\right)}{\left(w(\cdot),u_{t}(\cdot,t_{c})\right)}.$
(18)

From the weak formulation of the PDE (4),

$\displaystyle\left(w(\cdot),u_{t}(\cdot,t_{c})\right)$
$\displaystyle=-\left(L_{2}w(\cdot),L_{1}u(\cdot,t_{c})\right)+\left(w(\cdot),f(u,\cdot,t_{c})\right)$
(19)
$\displaystyle=-\left(L_{2}w(\cdot),L_{1}U(\cdot,t_{c})\right)-\left(L_{2}w(\cdot),L_{1}e(\cdot,t_{c})\right)+\left(w(\cdot),f(u,\cdot,t_{c})\right).$

Using Taylor’s Theorem once more to expand $f(u,x,t_{c})$ around $U(x,t_{c})$,

$\displaystyle\left(w(\cdot),f(u,\cdot,t_{c})\right)$
$\displaystyle=\left(w(\cdot),f(U,\cdot,t_{c})\right)+\left(\nabla_{u}f^{\top}(U,\cdot,t_{c})w(\cdot),e(\cdot,t_{c})\right)+\left(w(\cdot),\mathcal{R}_{2}(u,U,\cdot,t_{c})\right),$
(20)

where the remainder $\mathcal{R}_{2}(u,U,x,t_{c})$ is of order
$\mathcal{O}\left((u-U)^{2}\right)$. Substituting (20) into (19) and then into
(18) yields the final result.

###### Corollary 1

If the differential equation (1) is linear, so that $f=f(x,t)$ does not depend
on $u$, then the error representation is

$t_{t}-t_{c}=\frac{\left(w(\cdot),e(\cdot,t_{c})\right)+\left(w(\cdot),\mathcal{R}_{1}(\cdot,t_{c},t_{t})\right)}{\left(L_{2}w(\cdot),L_{1}U(\cdot,t_{c})\right)-\left(w(\cdot),f(\cdot,t_{c})\right)+\left(L_{2}w(\cdot),L_{1}e(\cdot,t_{c})\right)}.$
(21)

###  3.2 Error estimation for a linear functional

In order to define a computable error estimate, standard adjoint-based _a
posteriori_ error analysis is employed to approximate the three terms in (15)
that are linear functionals of $e(x,t)$. The standard analysis is presented
for completeness in Theorem 3.2 below.

###### Theorem 3.2

Given a numerical solution $U(x,t)$ to (1) and data $\psi(x)$, for any
$\hat{t}\in(0,T]$ the error $\left(\psi(\cdot),e(\cdot,\hat{t})\right)$ is
given by

$\left(\psi(\cdot),e(\cdot,\hat{t})\right)=\left(\phi(\cdot,0),e(\cdot,0)\right)+\int_{0}^{\hat{t}}\left(\phi(\cdot,t),f(U,\cdot,{t})-U_{t}(\cdot,{t})\right)\;{\rm
d}t-\int_{0}^{\hat{t}}\left(L_{2}\phi(\cdot,{t}),L_{1}U(\cdot,{t})\right)\;{\rm
d}t$ (22)

where $\phi(x,t)$ is the solution of the adjoint problem

$\displaystyle-\phi_{t}(x,t){+L^{*}\phi(x,t)}$
$\displaystyle=A_{u,U}^{*}\phi(x,t),\quad$ $\displaystyle x\in\Omega,\quad$
$\displaystyle t\in[0,\hat{t}),$ (23) $\displaystyle\phi(x,t)$
$\displaystyle=0,\qquad\qquad\;$ $\displaystyle x\in\partial\Omega,\quad$
$\displaystyle t\in[0,\hat{t}),$ $\displaystyle\phi(x,\hat{t})$
$\displaystyle=\psi(x),\qquad\quad$ $\displaystyle x\in\Omega.$

The operator $A_{u,U}^{*}$ is the adjoint of the linear operator

$A_{u,U}=\int_{0}^{1}\frac{df}{dz}(z,x,t)\;{\rm d}s,$ (24)

where $z=su+(1-s)U$.

###### Proof

Multiplying the adjoint equation (23) by the error $e(x,t)=u(x,t)-U(x,t)$, and
integrating over the space-time domain $\Omega\times[0,\hat{t}]$,

$0=\int_{0}^{\hat{t}}\left(\phi_{t}(\cdot,t),e(\cdot,t)\right)\;{\rm
d}t-\int_{0}^{\hat{t}}\left(L^{*}\phi(\cdot,t),e(\cdot,t)\right)\;{\rm
d}t+\int_{0}^{\hat{t}}\left(A_{u,U}^{*}\phi(\cdot,t),e(\cdot,t)\right)\;{\rm
d}t.$ (25)

Integrating the first term of (25) by parts in time and enforcing the initial
condition $\phi(x,\hat{t})=\psi(x)$,

$\int_{0}^{\hat{t}}\left(\phi_{t}(\cdot,t),e(\cdot,t)\right)\;{\rm
d}t={\left(\psi(\cdot),e(\cdot,\hat{t})\right)-\left(\phi(\cdot,0),e(\cdot,0)\right)-\int_{0}^{\hat{t}}\left(\phi(\cdot,t),e_{t}(\cdot,t)\right)\;{\rm
d}t.}$

From (5), $L^{*}=L_{1}^{*}L_{2}$, and the second term of (25) becomes

$\int_{0}^{\hat{t}}\left(L^{*}\phi(\cdot,t),e(\cdot,t)\right)\;{\rm
d}t=\int_{0}^{\hat{t}}\left(L_{1}^{*}L_{2}\phi(\cdot,t),e(\cdot,t)\right)\;{\rm
d}t=\int_{0}^{\hat{t}}\left(L_{2}\phi(\cdot,t),L_{1}e(\cdot,t)\right)\;{\rm
d}t.$ (26)

The operator $A_{u,U}$ has the property

$A_{u,U}e(x,t)=f(u,x,t)-f(U,x,t),$ (27)

hence, using the property of the adjoint (6) and (27), the third term of (25)
becomes

$\int_{0}^{\hat{t}}\left(A_{u,U}^{*}\phi(\cdot,t),e(\cdot,t)\right)\;{\rm
d}t=\int_{0}^{\hat{t}}\left(\phi(\cdot,t),A_{u,U}e(\cdot,t)\right)\;{\rm
d}t=\int_{0}^{\hat{t}}\left(\phi(\cdot,t),f(u,\cdot,t)-f(U,\cdot,t)\right)\;{\rm
d}t.$ (28)

Combining (25), (26), and (28) yields,

$\displaystyle 0$
$\displaystyle={\left(\psi(\cdot),e(\cdot,\hat{t})\right)-\left(\phi(\cdot,0),e(\cdot,0)\right)-\int_{0}^{\hat{t}}\left(\phi(\cdot,t),e_{t}(\cdot,t)\right)\;{\rm
d}t},$ $\displaystyle\hskip
28.45274pt-\int_{0}^{\hat{t}}\left(L_{2}\phi(\cdot,t),L_{1}e(\cdot,t)\right)\;{\rm
d}t+\int_{0}^{\hat{t}}\left(\phi(\cdot,t),f(u,\cdot,t)-f(U,\cdot,t)\right)\;{\rm
d}t.$

Since $e(x,t)=u(x,t)-U(x,t)$

$\displaystyle\left(\psi(\cdot),e(\cdot,\hat{t})\right)$
$\displaystyle=\left(\phi(\cdot,0),e(\cdot,0)\right)+\int_{0}^{\hat{t}}\left(\phi(\cdot,t),u_{t}(\cdot,t)-U_{t}(\cdot,t)\right)\;{\rm
d}t$ $\displaystyle\hskip
28.45274pt{+}\int_{0}^{\hat{t}}\left(L_{2}\phi(\cdot,t),L_{1}u(\cdot,t)\,{-}\,L_{1}U(\cdot,t)\right)\;{\rm
d}t\;{-}\int_{0}^{\hat{t}}\left(\phi(\cdot,t),f(u,\cdot,t)-f(U,\cdot,t)\right)\;{\rm
d}t.$

Noting that the exact solution $u(x,t)$ satisfies (4) produces the desired
result.

In practice, since operator $A_{u,U}$ requires the true solution to (1), it is
approximated by

$A_{u,U}\approx\nabla_{u}f(U,x,t).$

With this, the right side of the adjoint equation (23) is approximated by

$A_{u,U}^{*}\phi(x,t)\approx\left(\nabla_{u}f(U,x,t)\right)^{\top}\phi(x,t).$

###### Remark 1

Theorem 3.2 is formulated in terms of the strong form of the adjoint problem
(23), but clearly it is adequate for $\phi$ to satisfy the weak form of the
problem

$-\left(\phi(\cdot,t),v(\cdot)\right)+\left(L_{2}\phi(\cdot,t),L_{1}v(\cdot)\right)=\left(\phi(\cdot,t),A_{u,U}v(\cdot)\right)\qquad\forall
v\in H_{0}^{1}(\Omega),\forall t\in(0,\hat{t}].$ (29)

###  3.3 A computable error estimate for the time to an event $G(u;t)=R$

The magnitudes of the remainders $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ in
the error representation (15) decrease more rapidly than other terms as the
accuracy of the numerical solution increases (e.g., by increasing
$N_{x},N_{t},q_{t}$ or $q_{s}$ in § 2.2). Setting both
$\left(w(\cdot),\mathcal{R}_{1}(\cdot,t_{c},t_{t})\right)\approx 0$ and
$\left(w(\cdot),\mathcal{R}_{2}(u,U,\cdot,t_{c})\right)\approx 0$ gives a
first approximation

$\displaystyle t_{t}-t_{c}$
$\displaystyle\approx\frac{\left(w(\cdot),e(\cdot,t_{c})\right)}{\mathcal{D}(U,t_{c})},$
(30)

where $\mathcal{D}(U,t_{c})$ is defined in (14). The terms in the numerator
and denominator containing the unknown error $e(x,t)$ may be estimated using
Theorem 3.2. However, each of these expressions have different stability
properties, their estimation requires the solution to distinct adjoint
problems as given below.

First adjoint problem To estimate $\left(w(\cdot),e(\cdot,t_{c})\right)$, (23)
is solved backwards from $t_{c}$ with initial condition
$\psi(x)=\psi^{(1)}(x)=w(x)$. Substituting the solution $\phi^{(1)}$ in (22)
provides a computable estimate for

$\mathcal{E}_{1}\approx\left(\psi^{(1)}(\cdot),e(\cdot,t_{c})\right)=\left(w(\cdot),e(\cdot,t_{c})\right).$

Second adjoint problem To estimate
$\left(L_{2}w(\cdot),L_{1}e(\cdot,t_{c})\right)$, we assume $w(x)$ is smooth
enough for $L^{*}w$ to be well-defined and solve (23) backwards from $t_{c}$
with initial condition $\psi(x)=\psi^{(2)}(x)=L^{*}w(x)$. Substituting the
solution $\phi^{(2)}$ in (22) provides a computable estimate for

$\mathcal{E}_{2}\approx\left(\psi^{(2)}(\cdot),e(\cdot,t_{c})\right)=\left(L^{*}w(\cdot),e(\cdot,t_{c})\right)=\left(L_{2}w(\cdot),L_{1}e(\cdot,t)\right).$

Third adjoint problem To estimate
$\left(\nabla_{u}f^{\top}(U,\cdot,t_{c})w(\cdot),e(\cdot,t_{c})\right)$, (23)
is solved backwards with initial condition
$\psi(x)=\psi^{(3)}(x)=\nabla_{u}f^{\top}(U,x,t_{c})w(x)$. Substituting the
solution $\phi^{(3)}$ in (22) provides a computable estimate for

$\mathcal{E}_{3}\approx\left(\psi^{(3)}(\cdot),e(\cdot,t_{c})\right)=\left(\nabla_{u}f^{\top}(U,\cdot,t_{c})w(\cdot),e(\cdot,t_{c})\right).$

The third adjoint problem is only required for nonlinear problems. If the
right side in (1) is $f=f(x,t)$, the gradient is $\nabla_{u}f=0$ and
$\mathcal{E}_{3}\equiv 0$.

###### Remark 2

If the numerical solution $U$ is not sufficiently accurate, the $n$-th event
for the approximate solution $U$ may be closest to the $m$-th event for the
true solution $u$ where $n\neq m$. If this is the case, approximation by
Taylor series around $t_{c}$ and $U(x,t_{c})$ is inappropriate and the error
estimator may not be accurate. See Remark 2 in ChaEstSteTav21 .

##  4 Model problems

We apply the abstract theory developed in §3 to the heat equation and the 1D
linearized shallow water equations.

###  4.1 Heat equation

The heat equation models the diffusion of heat or chemical species. We
consider the heat equation in a domain $\Omega\subset\mathbb{R}^{n}$,

$\displaystyle u_{t}(x,t)-\nabla^{2}u(x,t)$ $\displaystyle=f(u,x,t),\qquad$
$\displaystyle x\in\Omega,\qquad t\in(0,T],$ (31) $\displaystyle u(x,t)$
$\displaystyle=0,\qquad$ $\displaystyle x\in\partial\Omega,\
\quad{t\in(0,T],}$ $\displaystyle u(x,0)$ $\displaystyle=u_{0}(x),\qquad$
$\displaystyle x\in\Omega,$

where $u(x,t)$ is the temperature of the medium or concentration of a chemical
species.

Weak form

The weak form is: Find $u\in L^{2}(0,T;H_{0}^{1}(\Omega))$ such that

$\left(u_{t},v\right)+\left(\nabla u,\nabla
v\right)=\left(f,v\right)\quad\forall v\in H_{0}^{1}(\Omega),\qquad
t\in(0,T].$

When written in standard form as $u_{t}+Lu=f$, the operator $L$ is
$-\nabla^{2}$. For any function $v\in H_{0}^{1}(\Omega)$,

$\left(Lu,v\right)=\left(-\nabla^{2}u,v\right)=\left(\nabla u,\nabla
v\right),$

implying

$L_{1}=\nabla=L_{2}.$

Error estimate

The error estimate (30) applied to the heat equation (31) becomes

$t_{t}-t_{c}\approx\frac{\left(w(\cdot),e(\cdot,t_{c})\right)}{\left(\nabla
w(\cdot),\nabla
U(\cdot,t_{c})\right)-\left(w(\cdot),f(U,\cdot,t_{c})\right)-\left(\nabla_{u}f^{\top}(U,\cdot,t_{c})w(\cdot),e(\cdot,t_{c})\right)+\left(\nabla
w(\cdot),\nabla e(\cdot,t_{c})\right)}.$

The three adjoint problems presented in §3.3 all have the same form,

$\left.\begin{aligned}
-\phi^{(i)}_{t}(x,t)-\nabla^{2}\phi^{(i)}(x,t)&=\nabla_{u}f^{\top}\phi,&&x\in\Omega,&&&t\in(t_{c},0]\\\
\phi^{(i)}(x,t)&=0,&&x\in\partial\Omega,&&&t\in(t_{c},0]\\\
\phi^{(i)}(x,t_{c})&=\psi^{(i)}(x),&&x\in\Omega\end{aligned}\qquad\right\\}\quad
i=1,2,3,$

but different initial conditions:

$\displaystyle\psi^{(1)}(x)$ $\displaystyle=w(x),$
$\displaystyle\psi^{(2)}(x)$ $\displaystyle={-}\nabla^{2}w{(x)},$
$\displaystyle\psi^{(3)}(x)$
$\displaystyle=\nabla_{u}f^{\top}(U,x,t_{c})w(x).$

###  4.2 1D linearized shallow water equation (SWE)

The shallow water equations model wave propagation in a fluid domain in which
the horizontal scale is large compared to the vertical scale. Applications of
the SWE arise frequently in the study of the ocean and atmosphere. In
particular, the linearized SWE model have been used to model tsunami wave
propagation and inundation in coastal areas, flooding from a dam break, and
flows and waves in the atmosphere, see e.g., Buhler98 ; CarYeh05 . A
posteriori error analysis for QoIs that are functionals of the solution to the
linearized SWE has been conducted previously in DavLeV16 .

The nonlinear shallow water equations in one dimension on a domain
$\Omega=[0,D]$ are (see DavLeV16 )

$\left.\begin{aligned} \eta_{t}(x,t)+\mu_{x}(x,t)&=f_{1}(x,t),\\\
\mu_{t}(x,t)+\left(\frac{\mu(x,t)^{2}}{\eta(x,t)}+\frac{1}{2}g\eta(x,t)^{2}\right)_{x}+g\eta
B_{x}(x)&=f_{2}(x,t),\end{aligned}\qquad\right\\}\qquad x\in\Omega,\qquad
t\in(0,T].$ (32)

where $\eta(x,t)$ corresponds to the water surface elevation and $\mu(x,t)$
corresponds to the momentum of the fluid. The constant $g$ is the
gravitational acceleration constant, $f(x,t)=(f_{1}(x,t),f_{2}(x,t))^{\top}$
is the forcing (typically $0$) and $B(x)$ is the bathymetry (bottom surface
profile).

Linearizing (32) around a flat fluid surface with elevation $\bar{\eta}$ and
momentum $\bar{\mu}=0$,

$\displaystyle\left.\begin{aligned} \eta_{t}(x,t)+\mu_{x}(x,t)&=f_{1}(x,t),\\\
\mu_{t}(x,t)+g\bar{h}(x)\eta_{x}(x,t)&=f_{2}(x,t),\\\
\end{aligned}\qquad\right\\}\qquad$ $\displaystyle x\in\Omega,\hskip
5.69054ptt\in(0,T],$ $\displaystyle\left.\begin{aligned} \mu(x,t)&=0,\hskip
15.93355pt\end{aligned}\qquad\right\\}\hskip 22.05087pt$ $\displaystyle
x\in\partial\Omega,t\in(0,T],$ (33) $\displaystyle\left.\begin{aligned}
\eta(x,0)&=\eta_{0}(x),\hskip 7.8245pt\\\ \mu(x,0)&=\mu_{0}(x),\\\
\end{aligned}\qquad\right\\}\qquad$ $\displaystyle x\in\Omega,$

where $\bar{h}(x)=\bar{\eta}-B(x)$, $\eta_{0}(x)$ is the initial surface
profile, $\mu_{0}(x)$ is the initial momentum profile, and $\eta$ satisfies
reflective boundary conditions (see DavLeV16 ).

Weak form

The weak form of (4.2) is: Find $\eta\in L^{2}(0,T;H^{1}(\Omega))$ and $\mu\in
L^{2}(0,T;H_{0}^{1}(\Omega))$

$\left(\eta_{t},v_{1}\right)+\left(\mu_{t},v_{2}\right)+\left(\mu_{x},v_{1}\right)+g\bar{h}\left(\eta_{x},v_{2}\right)=\left(f_{1},v_{1}\right)+\left(f_{2},v_{2}\right),\quad\forall{v_{1}\in
H^{1}(\Omega),}\,v_{2}\in H^{1}_{0}(\Omega),\quad t\in(0,T).$

When written in standard form as $u_{t}+Lu=f$, where
$u(x,t)=(\eta(x,t),\mu(x,t))^{\top}$ and the operator $L$ is

$L=\begin{pmatrix}0&\frac{\partial}{\partial x}\\\
g\bar{h}(x)\frac{\partial}{\partial x}&0\end{pmatrix},$

and since $L$ contains only first order derivatives, the operators $L_{1}$ and
$L_{2}$ are

$L_{1}=L=\begin{pmatrix}0&\frac{\partial}{\partial x}\\\
g\bar{h}(x)\frac{\partial}{\partial x}&0\end{pmatrix},\qquad
L_{2}=\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}.$

Error estimate

For the linearized shallow water equations (4.2), Corollary 1 is applicable
and the error estimate (30) becomes

$t_{t}-t_{c}\approx\frac{\left(w(\cdot),e(\cdot,t_{c})\right)}{\left(w(\cdot),A(\cdot)U_{x}(\cdot,t_{c})\right)-\left(w(\cdot),f(\cdot,t_{c})\right)+\left(w(\cdot),A(\cdot)e_{x}(\cdot,t_{c})\right)},$

where

$u=\begin{pmatrix}\eta\\\ \mu\end{pmatrix},\quad f=\begin{pmatrix}f_{1}\\\
f_{2}\end{pmatrix},\quad A(x)=\begin{pmatrix}0&1\\\
g\bar{h}(x)&0\end{pmatrix}.$

The two non-trivial adjoint problems take the form

$\left.\begin{aligned}
{-}\phi^{(i)}_{t}(x,t){-(}A^{\top}(x)\phi^{(i)}{)}_{x}(x,t)&=0,&&x\in\Omega,&&&t\in(t_{c},0],\\\
\phi^{(i)}(x,t)&=0,\quad&&x\in\partial\Omega,&&&t\in(t_{c},0],\\\
\phi^{(i)}(x,t_{c})&=\psi^{(i)}(x),&&x\in\Omega\end{aligned}\qquad\right\\}\quad
i=1,2,$

where

$\phi^{(i)}=\begin{pmatrix}\phi^{(i)}_{1}\\\ \phi^{(i)}_{2}\end{pmatrix}$

with initial conditions

$\displaystyle\psi^{(1)}(x)$ $\displaystyle=w(x),$
$\displaystyle\psi^{(2)}(x)$ $\displaystyle={-}(A^{\top}(x)w(x))_{x}.$

##  5 Numerical results

All numerical examples were computed using a uniform one-dimensional spatial
mesh of length $D$ with $N_{x}$ finite elements and a uniform temporal mesh of
length $T$ with $N_{t}$ elements. The forward problems were solved using a
space-time Galerkin method as described §2.2. Unless otherwise noted, the same
degree of approximation was used for the basis functions in both space and
time. In all examples, the adjoint problems were solved using the same space-
time mesh as the forward problem, but with a Galerkin finite element method
_two_ degrees higher in both space and time. We choose to solve the adjoint
problem with this precision in order to demonstrate the accuracy of the
estimates. In practice, when estimates are used to drive an adaptive strategy,
lower accuracy may be acceptable.

###  5.1 1D linear heat equation

Consider the heat equation in one dimension,

$\displaystyle u_{t}-u_{xx}$ $\displaystyle=\sin(\pi
x)(\pi^{2}\cos(t)-\sin(t)),\quad$ $\displaystyle x\in[0,1],$ $\displaystyle
t\in(0,0.5],$ (34) $\displaystyle u(0,t)$ $\displaystyle=u(1,t)=0,$
$\displaystyle t\in(0,0.5],$ $\displaystyle u(x,0)$ $\displaystyle=\sin(\pi
x),$ $\displaystyle x\in[0,1].$

The right-hand-side forcing function has been chosen so that (34) has the
analytical solution $u(x,t)=\cos(t)\sin(\pi x)$. Let the functional $G(u;t)$
be

$G(u;t)=\left(w(\cdot),u(\cdot,t)\right)=\left(\sin(\pi(\cdot)),u(\cdot,t)\right),$

and the threshold value be $R=0.47$. The true values of the functional
$G(u;t)$ are provided in Figure 2a and the threshold value $R$ appears as a
dashed horizontal line. The QoI is the first occurrence of the event
$G(u;t)=R$, that is, $t_{t}=H(u,0)$.

(a) §5.1 and §5.2.

(b) §5.3.

Figure 2: Reference values of the functional $G(u;t)$ for §5.1, §5.2, and
§5.3.

The “initial” conditions for the adjoint problems in §4.1 are

$\displaystyle\psi^{(1)}(x)$ $\displaystyle=w(x)=\sin(\pi x),$
$\displaystyle\psi^{(2)}(x)$ $\displaystyle={-}w_{xx}{(x)}=\pi^{2}\sin(\pi
x),$ $\displaystyle\psi^{(3)}(x)$
$\displaystyle=\nabla_{u}f^{\top}(U,x,t_{c})w(x)=0.$

The problem (34) was solved using the cG(1,1) method on a sequence of
increasingly refined uniform meshes. The computed time $t_{c}$ to the first
occurrence of the event $G(u;t)=R$, the error in the QoI $e_{Q}$, and the
effectivity ratio $\rho_{\rm eff}$, of the _a posteriori_ error estimator are
presented in Table 1. The error decreases as the mesh is refined as expected,
and the effectivity ratios are close to 1 in all cases.

$N$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$
---|---|---|---
50 | 0.346346 | 1.820 $\times 10^{-3}$ | 1.003
100 | 0.347711 | 4.546 $\times 10^{-4}$ | 1.001
200 | 0.348053 | 1.129 $\times 10^{-4}$ | 1.000
400 | 0.348138 | 2.839 $\times 10^{-5}$ | 1.000
Table 1: Computed time $t_{c}$ to the first occurrence of the event
$G(u;t)=R$, error in the QoI $e_{Q}$, and effectivity ratio $\rho_{\rm eff}$,
for §5.1, using cG(1,1) to solve the forward problem and cG(3,3) to solve the
adjoint problems. The exact time to the first occurrence of the event is
$t_{t}=0.34816603$.

###  5.2 1D nonlinear heat equation

Consider the nonlinear heat equation in one dimension

$\displaystyle u_{t}-u_{xx}$ $\displaystyle=-u^{2}+\sin(\pi
x)(-\sin(t)+\pi^{2}\cos(t)+\cos^{2}(t)\sin(\pi x)),\quad$ $\displaystyle
x\in(0,1],$ $\displaystyle t\in(0,0.5],$ (35) $\displaystyle u(0,t)$
$\displaystyle=u(1,t)=0$ $\displaystyle t\in(0,0.5],$ $\displaystyle u(x,0)$
$\displaystyle=\sin(\pi x)$ $\displaystyle x\in[0,1].$

The right-hand-side “forcing” function has been chosen so that the nonlinear
heat equation (35) has the same analytical solution as the linear heat
equation (34) in §5.1, namely $u(x,t)=\cos(t)\sin(\pi x)$. The functional
$G(u;t)$ is again chosen to be

$G(u;t)=\left(w(\cdot),u(\cdot,t)\right)=\left(\sin(\pi(\cdot)),u(\cdot,t)\right),$

and the threshold value to be $R=0.47$. The QoI is the first occurrence of the
event $G(u;t)=R$, that is, $t_{t}=H(u,0)$. The true values of the functional
are graphed in Figure 2a.

The “initial” conditions for the adjoint problems in §4.1 are

$\displaystyle\psi^{(1)}(x)$ $\displaystyle=w(x)=\sin(\pi x),$
$\displaystyle\psi^{(2)}(x)$ $\displaystyle={-}w_{xx}{(x)}=\pi^{2}\sin(\pi
x),$ $\displaystyle\psi^{(3)}(x)$
$\displaystyle=\nabla_{u}f^{\top}(U,x,t_{c})w(x)=-2\sin(\pi x)U(x,t_{c}).$

Note in particular, that unlike for the previous, linear case in §5.1, the
third adjoint problem is non-trivial. The problem (35) was solved using the
cG(1,1) methods using a sequence of increasingly refined uniform meshes. The
computed time $t_{c}$ to the first occurrence of the event $G(u;t)=R$, the
error in the QoI $e_{Q}$, and the effectivity ratio $\rho_{\rm eff}$, of the
_a posteriori_ error estimator are presented in Table 2 which clearly
indicates that the error decreases as the mesh is refined and the effectivity
ratios remains close to 1 in all cases.

$N$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$
---|---|---|---
50 | 0.346531 | 1.635 $\times 10^{-3}$ | 1.002
100 | 0.347757 | 4.087 $\times 10^{-4}$ | 1.001
200 | 0.348065 | 1.015 $\times 10^{-4}$ | 1.000
400 | 0.348140 | 2.553 $\times 10^{-5}$ | 1.000
Table 2: Computed time $t_{c}$ to the first occurrence of the event
$G(u;t)=R$, error in the QoI $e_{Q}$, and effectivity ratio $\rho_{\rm eff}$,
for §5.2, using cG(1,1) to solve the forward problem and cG(3,3) to solve the
adjoint problems. The exact time to the first occurrence of the event is
$t_{t}=0.34816603$.

###  5.3 1D linearized SWE: Manufactured solution

Consider the linearized SWEs (4.2) over the space-time domain
$(x,t)\in[0,10]\times[0,1]$, with right hand side

$f=\begin{pmatrix}f_{1}\\\ f_{2}\end{pmatrix}=\begin{pmatrix}-\sin(t)\sin(\pi
x)+\pi\cos(t)\cos(\pi x)\\\ -\sin(t)\sin(\pi x)+\pi g\bar{h}(x)\cos(t)\cos(\pi
x)\end{pmatrix}.$ (36)

In order to construct an analytical exact solution, the bathymetry was chosen
to be a constant, $B(x)=-10$ as shown in Figure 3a and Dirichlet boundary
conditions were imposed at the boundaries of the domain, specifically
$u(0,t)=u(10,t)=(2,0)^{\top}$. The initial surface elevation and momentum were
chosen as $u(x,0)=(\eta_{0}(x),\mu_{0}(x))^{\top}=(2+\sin(\pi x),\sin(\pi
x))^{\top}$. The initial surface elevation $\eta_{0}(x)$ is shown in Figure
3b. The rest-state elevation was $\bar{\eta}=2$.

(a) $B(x)$

(b) $\eta_{0}(x)$

(c) $w_{2}(x)$

Figure 3: The bathymetry $B(x)$, initial wave height $\eta_{0}(x)$, and
nonzero component of the weight function $w_{2}(x)$ for §5.3.

Finally, the forcing term (36) was chosen so that (4.2) has the analytical
solution

$u(x,t)=\begin{pmatrix}\eta(x,t)\\\
\mu(x,t)\end{pmatrix}=\begin{pmatrix}2+\cos(t)\sin(\pi x)\\\ \cos(t)\sin(\pi
x)\end{pmatrix}.$

In order to define the functional $G(u;t)$, let the weight function be
$w(x)=(0,w_{2}(x))$ where

$w_{2}(x)=\qquad\left\\{\qquad\begin{aligned} 0\qquad\qquad\qquad&x\leq 5,\\\
10(x-5)^{2}(x-6)^{2}\quad 5<&x<6,\\\ 0\qquad\qquad\qquad&x\geq
6.\end{aligned}\qquad\right\\}$

The non-zero component of the weight function $w_{2}(x)$ is depicted in Figure
3c. The true values of the functional
$G(u;t)=\left(w(\cdot),u(\cdot,t)\right)$ are provided in Figure 2b, where the
threshold value $R=-0.19$ appears as a dashed horizontal line. The QoI is the
first occurrence of the event $G(u;t)=R$, that is, $t_{t}=H(u,0)$.

The initial conditions for the adjoint problems in §4.2 are

$\displaystyle\psi^{(1)}(x)$ $\displaystyle=(0,w_{2}(x))^{\top},$
$\displaystyle\psi^{(2)}(x)$
$\displaystyle=({-}g(\bar{h}(x)w_{2}(x))_{x},0)^{\top}.$

The linear SWE were solved using the cG(2,2) method (§2.2) on a sequence of
increasingly refined uniform meshes. The computed time $t_{c}$ to the first
occurrence of the event $G(u;t)=R$, the error in the QoI $e_{Q}$, and the
effectivity ratio $\rho_{\rm eff}$, of the _a posteriori_ error estimator are
presented in Table 3.

$N$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$
---|---|---|---
50 | 0.820996 | -1.120 $\times 10^{-3}$ | 0.997
100 | 0.819917 | -4.035 $\times 10^{-5}$ | 1.000
200 | 0.819883 | -6.940 $\times 10^{-6}$ | 1.000
400 | 0.819878 | -1.851 $\times 10^{-6}$ | 1.000
Table 3: Computed time $t_{c}$ to the first occurrence of the event
$G(u;t)=R$, error in the QoI $e_{Q}$, and effectivity ratio $\rho_{\rm eff}$,
for §5.3, using cG(2,2) to solve the forward problem and cG(4,4) to solve the
adjoint problems. The exact time to the first occurrence of the event is
$t_{t}=0.81987644$.

###  5.4 1D linearized SWE: Constant bathymetry

Consider the linearized SWEs (4.2) over the space-time domain
$(x,t)\in[0,400]\times[0,200]$, with the more realistic right hand side
$f=(0,0)^{\top}$. The bathymetry was chosen to be the constant $B(x)=-0.1$ as
depicted in Figure 4a and the rest-state height was $\bar{\eta}=1$. Dirichlet
boundary conditions were imposed on the momentum only, specifically
$\mu(0,t)=\mu(D,t)=0$. The initial conditions were
$u(x,0)=(\eta_{0}(x),0)^{\top}$, where

$\eta_{0}(x)=\qquad\left\\{\qquad\begin{aligned}
0,\qquad\qquad\qquad\qquad&x\leq 100,\\\
\frac{0.4}{390625}(x-100)^{2}(x-150)^{2},\quad 100<&x<150,\\\
0,\qquad\qquad\qquad\qquad&x\geq 150.\end{aligned}\qquad\right\\}$

The initial surface elevation $\eta_{0}(x)$ is shown in Figure 4b.

In order to define the functional $G(u;t)$ we set the weight function
$w=(w_{1}(x),0)^{\top}$ where

$w_{1}(x)=\qquad\left\\{\qquad\begin{aligned} 0,\qquad\qquad\qquad\qquad&x\leq
160,\\\ \frac{1}{200000}(x-160)^{2}(x-200)^{2},\quad 160<&x<200,\\\
0,\qquad\qquad\qquad\qquad&x\geq 200.\end{aligned}\qquad\right\\}$

The non-zero component of the weight function $w_{1}(x)$ is shown in Figure
4c.

(a) $B(x)$

(b) $\eta_{0}(x)$

(c) $w_{1}(x)$

Figure 4: The bathymetry $B(x)$, initial wave height $\eta_{0}(x)$, and
nonzero component of the weight function $w_{1}(x)$ for §5.4.

The linear SWE were solved using the cG(2,2) method (§2.2), and the surface
elevation $\eta(x,t)$ at several time is shown in Figure 5.

(a) $t=0$

(b) $t=20$s

(c) $t=40$s

(d) $t=65$s

(e) $t=90$s

(f) $t=107.5$s

(g) $t=125$s

(h) $t=165$s

(i) $t=200$s

Figure 5: Evolution of the surface elevation $\eta(x,t)$ for §5.4. The initial
profile generates left and right traveling waves which eventually reflect off
the left-hand and right-hand boundaries. The arrows indicate the direction of
propagation.

Since an analytical solution $u$ is unavailable, a reference solution $u$ was
computed using a mesh-size of $N_{\rm ref}=800$ with the cG(3,3) method.
Reference values of the functional $G(u;t)$ are graphed in Figure 6a, where
the threshold value of $R=2$ is indicated by the dashed horizontal line

(a) §5.4.

(b) §5.5

(c) §5.6.

Figure 6: Reference values of the $G(u;t)$ for §5.4, §5.5 and §5.6.

The QoI defined by the functional $G(u;t)$ is essentially the time at which
the surface elevation at a narrow region centered at $x=180$ achieves a
specified value. Recall that by an event we mean $G(\cdot;t)$ achieving the
threshold value $R$. The first event occurs when the incident wave arrives
near $x=180$ and $G(u;t)$ increases above the threshold value. The second
event occurs when the incident wave travels beyond the narrow region centered
at $x=180$ and $G(u;t)$ drops below the threshold value. The third and fourth
events occur when the reflected wave passes through the narrow region centered
at $x=180$ and $G(u;t)$ first increases above and then drops below the
threshold value.

The initial conditions for the adjoint problems in §4.2 are

$\displaystyle\psi^{(1)}(x)$ $\displaystyle=(w_{1}(x),0)^{\top},$
$\displaystyle\psi^{(2)}(x)$ $\displaystyle=(0,{-}w_{1,x}(x))^{\top}.$

The computed times $t_{c}$ to the first, second and third occurrences of the
event $G(u;t)=R$, errors in the QoI $e_{Q}$, and effectivity ratios $\rho_{\rm
eff}$, are presented in Table 4, where the effectivity ratios are all seen to
be close to 1 even for the later events.

| First Event | Second Event | Third Event
---|---|---|---
$N$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$
50 | 13.4083 | -5.849 $\times 10^{-2}$ | 1.001 | 20.2736 | -1.204 $\times 10^{-1}$ | 1.008 | 90.0299 | -5.370 $\times 10^{-1}$ | 0.991
100 | 13.3634 | -1.363 $\times 10^{-2}$ | 1.000 | 20.1599 | -6.778 $\times 10^{-3}$ | 1.003 | 89.5458 | -5.291 $\times 10^{-2}$ | 1.000
200 | 13.3495 | 2.779 $\times 10^{-4}$ | 1.000 | 20.1529 | 2.538 $\times 10^{-4}$ | 0.990 | 89.4961 | -3.164 $\times 10^{-3}$ | 0.988
400 | 13.3499 | -1.459 $\times 10^{-4}$ | 1.000 | 20.1531 | 6.998 $\times 10^{-5}$ | 1.002 | 89.4931 | -2.166 $\times 10^{-4}$ | 0.967
Table 4: Computed times $t_{c}$ to the first, second and third occurrences of
the event $G(u;t)=R$, error in the QoI $e_{Q}$, and effectivity ratio
$\rho_{\rm eff}$, for §5.4, using cG(2,2) to solve the forward problem and
cG(4,4) to solve the adjoint problems. The first, second, and third events
occur at $t_{t}=13.349799$, $t_{t}=20.153171$, and $t_{t}=89.492912$
respectively.

###  5.5 1D linearized SWE: Continental shelf

Consider the linearized SWEs (4.2) over the space-time domain
$(x,t)\in[0,400000]\times[0,4200]$. In this example, the bathymetry $B(x)$ was
chosen to be

$B(x)=\qquad\left\\{\qquad\begin{aligned} -200,\qquad\qquad\qquad&x\leq
25000,\\\ -0.152x+3600,\quad 25000<&x<50000,\\\ -4000,\qquad\qquad\qquad&x\geq
50000,\end{aligned}\qquad\right\\}$

as depicted in Figure 7a and represents a continental shelf. The rest-state
height $\bar{\eta}=1$. Dirichlet boundary conditions were applied to the
momentum only, specifically $\mu(0,t)=\mu(D,t)=0$. The initial condition was
$u(x,0)=(\eta_{0}(x),0)^{\top}$, where

$\eta_{0}(x)=\qquad\left\\{\qquad\begin{aligned}
0,\qquad\qquad\qquad\qquad\qquad&x\leq 100000,\\\
\frac{0.4}{25000^{4}}(x-100000)^{2}(x-150000)^{2},\quad 100000<&x<150000,\\\
0,\qquad\qquad\qquad\qquad\qquad&x\geq 150000.\end{aligned}\qquad\right\\}$

The initial surface elevation $\eta_{0}(x)$ is depicted in Figure 7b.

In order to define the functional $G(u;t)$ (9), the weight function
$w=(w_{1}(x),0)^{\top}$ was

$w_{1}(x)=\qquad\left\\{\qquad\begin{aligned}
0,\qquad\qquad\qquad\qquad\qquad&x\leq 10000,\\\
\frac{1}{7500^{4}}(x-10000)^{2}(x-25000)^{2},\quad 10000<&x<25000,\\\
0,\qquad\qquad\qquad\qquad\qquad&x\geq 25000.\end{aligned}\qquad\right\\}$

The non-zero component of the weight function $w_{1}(x)$ is depicted in Figure
7c.

(a) $B(x)$

(b) $\eta_{0}(x)$

(c) $w_{1}(x)$

Figure 7: The bathymetry $B(x)$, initial wave height $\eta_{0}(x)$, and
nonzero component of the weight function $w_{1}(x)$ for §5.5.

A reference solution $u$ was computed using a mesh with $N_{\rm ref}=1280$
sub-intervals with the cG(3,3) method. The corresponding values of the
functional $G(u;t)$ are plotted in Figure 6b and the critical value of
$R=1000$ is presented as a dashed horizontal line. The QoIs we consider are
the times for the first, second and third occurrences of the event $G(u;t)=R$.

The initial conditions for the adjoint problems in §4.2 are

$\displaystyle\psi^{(1)}(x)$ $\displaystyle=(w_{1}(x),0)^{\top},$
$\displaystyle\psi^{(2)}(x)$ $\displaystyle=(0,{-}w_{1,x}(x))^{\top}.$

The forward problem was solved using the cG(2,2) method (§2.2) on a sequence
of increasingly refined uniform meshes. The computed times $t_{c}$ to the
first, second and third occurrences of the event $G(u;t)=R$, errors in the QoI
$e_{Q}$, and effectivity ratios $\rho_{\rm eff}$, are presented in Table 5.

| First Event | Second Event | Third Event
---|---|---|---
$N$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$
80 | 651.267 | -3.647 $\times 10^{-1}$ | 1.002 | 838.675 | -2.456 $\times 10^{-1}$ | 0.970 | 1438.541 | 9.255 $\times 10^{-1}$ | 0.925
160 | 650.975 | -7.329 $\times 10^{-2}$ | 1.000 | 838.421 | 8.597 $\times 10^{-3}$ | 0.929 | 1439.384 | 1.008 $\times 10^{-1}$ | 0.979
320 | 650.907 | -4.970 $\times 10^{-3}$ | 1.000 | 838.428 | 8.944 $\times 10^{-4}$ | 1.006 | 1439.479 | 6.161 $\times 10^{-3}$ | 0.986
640 | 650.902 | 3.672 $\times 10^{-4}$ | 1.000 | 838.429 | -2.890 $\times 10^{-4}$ | 1.012 | 1439.484 | 1.046 $\times 10^{-3}$ | 1.042
Table 5: Computed times $t_{c}$ to the first, second and third occurrences of
the event $G(u;t)=R$, error in the QoI $e_{Q}$, and effectivity ratio
$\rho_{\rm eff}$, for §5.5, using cG(2,2) to solve the forward problem and
cG(4,4) to solve the adjoint problems. The first, second, and third events
occur at $t_{t}=650.79647$, $t_{t}=838.42910$, and $t_{t}=1439.4849$
respectively.

The effect of the spatial discretization was investigated by repeating the
above calculations using cG(2,1). Note that the degree of the spatial basis
functions used to solve the adjoint problems was also reduced by one. As
anticipated, as is shown in Table 6, the errors $e_{Q}$ for the times to the
first three events $G(u;t)=R$ are much larger than when solving using cG(2,2).
However, the effectivity ratios remain close to 1 as desired.

| First Event | Second Event | Third event
---|---|---|---
$N$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$
80 | 643.874 | 6.922 $\times 10^{0}$ | 1.023 | 846.194 | -7.681 $\times 10^{0}$ | 0.999 | 1465.863 | -2.648 $\times 10^{1}$ | 1.050
160 | 650.852 | -5.536 $\times 10^{-2}$ | 1.014 | 840.674 | -2.161 $\times 10^{0}$ | 0.999 | 1442.571 | -3.192 $\times 10^{0}$ | 0.994
320 | 650.802 | -5.360 $\times 10^{-3}$ | 1.002 | 838.683 | -1.702 $\times 10^{-1}$ | 0.998 | 1439.675 | -2.964 $\times 10^{-1}$ | 1.003
640 | 650.795 | 1.889 $\times 10^{-3}$ | 0.997 | 838.521 | -8.420 $\times 10^{-3}$ | 0.995 | 1439.392 | -1.001 $\times 10^{-2}$ | 1.002
Table 6: Computed times $t_{c}$ to the first, second and third occurrences of
the event $G(u;t)=R$, error in the QoI $e_{Q}$, and effectivity ratio
$\rho_{\rm eff}$, for §5.5, using cG(2,1) to solve the forward problem and
cG(4,3) to solve the adjoint problems. The first, second, and third events
occur at $t_{t}=650.79647$, $t_{t}=838.51300$, and $t_{t}=1439.3789$
respectively.

###  5.6 1D linearized SWE: Reef

Consider the homogeneous, linearized SWEs (4.2) over the space-time domain
$(x,t)\in[0,400000]\times[0,4200]$. The bathymetry $B(x)$ was a simplified
representation of an obstruction on the sea floor, e.g., a reef, and is given
as

$B(x)=\qquad\left\\{\qquad\begin{aligned} -4000,\qquad\qquad\qquad\qquad&x\leq
25000,\\\ \frac{-50}{25000^{2}}(x-200000)(x-250000,)\quad 200000<&x<250000,\\\
-4000,\qquad\qquad\qquad\qquad&x\geq 50000,\end{aligned}\qquad\right\\}$

which is depicted in Figure 8a. The rest-state $\bar{\eta}=1$. The Dirichlet
boundary conditions were again $\mu(0,t)=\mu(D,t)=0$, and the initial
condition was $u(x,0)=(\eta_{0}(x),0)^{\top}$, where

$\eta_{0}(x)=\qquad\left\\{\qquad\begin{aligned}
0,\qquad\qquad\qquad\qquad\qquad&x\leq 100000,\\\
\frac{0.4}{25000^{4}}(x-100000)^{2}(x-150000)^{2},\quad 100000<&x<150000,\\\
0,\qquad\qquad\qquad\qquad\qquad&x\geq 150000.\end{aligned}\qquad\right\\}$

The initial wave height $\eta_{0}(x)$ is presented in Figure 8b.

In order to define the functional $G(u;t)$, the weight function
$w=(w_{1}(x),0)^{\top}$ was given by

$w_{1}(x)=\qquad\left\\{\qquad\begin{aligned}
0,\qquad\qquad\qquad\qquad\qquad&x\leq 10000,\\\
\frac{1}{7500^{4}}(x-10000)^{2}(x-25000)^{2},\quad 10000<&x<15000,\\\
0,\qquad\qquad\qquad\qquad\qquad&x\geq 25000.\end{aligned}\qquad\right\\}$

and the non-zero component of the weight function $w_{1}(x)$ is depicted in
Figure 8c.

(a) $B(x)$

(b) $\eta_{0}(x)$

(c) $w_{1}(x)$

Figure 8: The bathymetry $B(x)$, initial wave height $\eta_{0}(x)$, and
nonzero component of the weight function $w_{1}(x)$ for §5.6.

The reference solution was computed on a mesh with $N_{\rm ref}=1280$ sub-
intervals with the cG(3,3) method. The reference values of the functional
$G(u;t)$ are plotted in Figure 6c where the critical value of $R=1000$ is
presented as a dashed horizontal line. While this appears more complicated
than Figures 6a and 6b this is merely because the weight function is located
close to a domain boundary and so $G(u;t)$ does not have time to “recover” and
drop back to a small value before the reflected wave impinges on the critical
region. The QoIs we consider are the times for the first, second and third
occurrences of the event $G(u;t)=R$.

The initial conditions for the adjoint problems in §4.2 are

$\displaystyle\psi^{(1)}(x)$ $\displaystyle=(w_{1}(x),0)^{\top},$
$\displaystyle\psi^{(2)}(x)$ $\displaystyle=(0,{-}w_{1,x}(x))^{\top}.$

The forward problem was solved using the cG(2,2) method (§2.2) on a series of
increasingly fine uniform meshes. The computed times $t_{c}$ to the first,
second and third occurrences of the event $G(u;t)=R$, errors in the QoI
$e_{Q}$, and effectivity ratios $\rho_{\rm eff}$, are presented in Table 7,
where the effectivity ratios are all close to one.

| First event | Second event | Third event
---|---|---|---
$N$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$ | $t_{c}$ | $e_{Q}$ | $\rho_{\rm eff}$
80 | 489.353 | -3.398 $\times 10^{0}$ | 1.016 | 615.540 | -9.503 $\times 10^{0}$ | 1.260 | 663.188 | -6.692 $\times 10^{0}$ | 0.934
160 | 485.870 | 8.481 $\times 10^{-2}$ | 1.000 | 607.105 | -1.069 $\times 10^{0}$ | 1.010 | 656.738 | -2.412 $\times 10^{-1}$ | 0.981
320 | 485.951 | 4.482 $\times 10^{-3}$ | 1.000 | 606.055 | -1.869 $\times 10^{-2}$ | 1.043 | 656.478 | 1.899 $\times 10^{-2}$ | 0.942
640 | 485.954 | 7.146 $\times 10^{-4}$ | 1.001 | 606.038 | -1.194 $\times 10^{-3}$ | 1.016 | 656.501 | -4.511 $\times 10^{-3}$ | 0.960
Table 7: Computed times $t_{c}$ to the first, second and third occurrences of
the event $G(u;t)=R$, error in the QoI $e_{Q}$, and effectivity ratio
$\rho_{\rm eff}$, for §5.6, using cG(2,2) to solve the forward problem and
cG(4,4) to solve the adjoint problems. The first, second, and third events
occur at $t_{t}=485.95515$, $t_{t}=606.03632$, and $t_{t}=656.49655$
respectively.

##  6 Conclusions

We have developed and implemented an accurate adjoint-based a posteriori error
estimates for the time to an event, by which we mean the time for a functional
of the solution to a time-dependent PDEs to achieve a threshold value. These
estimates are based on the use of Taylor’s Theorem on the functional defining
the QoI, and multiple applications of standard adjoint-based error analysis,
to achieve a computable estimate. Implementations for the 1D heat equation and
the linearized shallow water equation provide examples of both parabolic and
hyperbolic PDEs. The error estimates provides a basis for developing
refinement strategies to reduce the error in the time to the event. A more
detailed analysis is required to partition the error in to contributions
arising from the independent spatial and temporal discretizations (see e.g.,
chaudhry2021posteriori ) and to determine the overall effect on the error from
terms arising in both the numerator and denominator of (30).

## Acknowledgments

J. Chaudhry’s and Z. Stevens’s work is supported by the NSF-DMS 1720402. D.
Estep’s and S. Tavener’s work is supported by NSF-DMS 1720473. D. Estep’s work
is supported by NSERC.

## References

  * (1) M. Ainsworth and T. Oden.  _A posteriori_ error estimation in finite element analysis.  John Wiley-Teubner, 2000.
  * (2) W. Bangerth and R. Rannacher.  Adaptive Finite Element Methods for Differential Equations.  Birkhauser Verlag, 2003.
  * (3) T. J. Barth.  _A posteriori_ Error Estimation and Mesh Adaptivity for Finite Volume and Finite Element Methods, volume 41 of Lecture Notes in Computational Science and Engineering.  Springer, New York, 2004.
  * (4) R. Becker and R. Rannacher.  An optimal control approach to _a posteriori_ error estimation in finite element methods.  Acta Numerica, pages 1–102, 2001.
  * (5) Oliver Bühler.  A shallow-water model that prevents nonlinear steepening of gravity waves.  Journal of the Atmospheric Sciences, 1998.
  * (6) Yang Cao and Linda Petzold.  _A posteriori_ error estimation and global error control for ordinary differential equations by the adjoint method.  SIAM Journal on Scientific Computing, 26(2):359–374, 2004.
  * (7) V. Carey, D. Estep, and S. Tavener.  _A posteriori_ analysis and adaptive error control for multiscale operator decomposition solution of elliptic systems I: Triangular systems.  SIAM Journal on Numerical Analysis, 47(1):740–761, 2008.
  * (8) George F. Carrier and Harry Yeh.  Tsunami propagation from a finite source.  Computer Modeling in Engineering & Sciences, 2005.
  * (9) J. H. Chaudhry, D. Estep, V. Ginting, and S. Tavener.  _A posteriori_ analysis for iterative solvers for non-autonomous evolution problems.  SIAM Journal on Uncertainty Quantification, 3, 2015.
  * (10) Jehanzeb Chaudhry, Donald Estep, and Simon Tavener.  _A posteriori_ error analysis for a space-time parallel discretization of parabolic partial differential equations.  arXiv preprint arXiv:2111.00606, 2021.
  * (11) J.H. Chaudhry.  _A posteriori_ analysis and efficient refinement strategies for the Poisson–Boltzmann equation.  SIAM Journal on Scientific Computing, 40(4):A2519—-A2542, 2018\.
  * (12) J.H. Chaudhry, J.B. Collins, and J.N. Shadid.  _A posteriori_ error estimation for multi-stage Runge-Kutta IMEX schemes.  Applied Numerical Mathematics, 117:36–49, Jul 2017.
  * (13) J.H. Chaudhry, D. Estep, V. Ginting, J.N. Shadid, and S. Tavener.  _A posteriori_ error analysis of IMEX multi-step time integration methods for advection–diffusion–reaction equations.  Computer Methods in Applied Mechanics and Engineering, 285:730–751, 2015.
  * (14) J.H. Chaudhry, D. Estep, S. Tavener, V. Carey, and J. Sandelin.  _A posteriori_ error analysis of two-stage computation methods with application to efficient discretization and the Parareal algorithm.  SIAM Journal on Numerical Analysis, 54(5):2974–3002, 2016.
  * (15) J.H. Chaudhry, J.N. Shadid, and T. Wildey.  _A posteriori_ analysis of an IMEX entropy-viscosity formulation for hyperbolic conservation laws with dissipation.  Applied Numerical Mathematics, 135, 2019.
  * (16) J.H. Chaudry, D. Estep, Z. Stevens, and S. Tavener.  Error estimation and uncertainty quantification for first time to a threshold value.  BIT Numerical Mathematics, 61:275–307, 2021.
  * (17) James B. Collins, Don Estep, and Simon Tavener.  _A posteriori_ error estimation for the Lax–Wendroff finite difference scheme.  Journal of Computational and Applied Mathematics, 263:299–311, 2014\.
  * (18) J.B. Collins, D. Estep, and S. Tavener.  _A posteriori_ error analysis for finite element methods with projection operators as applied to explicit time integration techniques.  BIT Numerical Mathematics, 55(4):1017–1042, 2015.
  * (19) B.N. Davis and R.J. LeVeque.  Adjoint methods for guiding adaptive mesh refinement in tsunami modeling.  Pure and Applied Geophysics, 173:4055–4074, 2016.
  * (20) D. Estep.  _A posteriori_ error bounds and global error control for approximation of ordinary differential equations.  SIAM J. Numer. Anal., 32(1):1–48, 1995.
  * (21) D. Estep.  Error estimates for multiscale operator decomposition for multiphysics models.  In J. Fish, editor, Multiscale methods: bridging the scales in science and engineering. Oxford University Press, USA, 2009.
  * (22) D. Estep, V. Ginting, and S. Tavener.  _A posteriori_ analysis of a multirate numerical method for ordinary differential equations.  Computer Methods in Applied Mechanics and Engineering, 223:10–27, 2012.
  * (23) D. Estep, M. Holst, and D. Mikulencak.  Accounting for stability: _a posteriori_ error estimates based on residuals and variational analysis.  Comm. Numer. Methods Engrg., 18:15–30, 2002.
  * (24) D. Estep, M. Larson, and R. Williams.  Estimating the error of numerical solutions of systems of reaction-diffusion equations.  Memoirs of the American Mathematical Society, 696, 07 2000.
  * (25) Lawrence C. Evans.  Partial differential equations.  American Mathematical Society, 2010.
  * (26) M.B. Giles and E. Süli.  Adjoint methods for pdes: _a posteriori_ error analysis and postprocessing by duality.  Acta Numerica, 11(1):145–236, 2002.
  * (27) Paul Houston.  Adjoint error estimation and adaptivity for hyperbolic problems.  In Handbook of Numerical Analysis, volume 18, pages 233–261. Elsevier, 2017.
  * (28) A. Johansson, J.H. Chaudhry, V. Carey, D. Estep, V. Ginting, M. Larson, and S.J. Tavener.  Adaptive finite element solution of multiscale PDE–ODE systems.  Computer Methods in Applied Mechanics and Engineering, 287:150–171, 2015.