# Analysis of the Impact of Unforced Errors in Tennis

Hashan Peiris, Nirodha Epasinghege Dona and Tim B. Swartz  H. Peiris is a PhD
candidate and T. Swartz is Professor, Department of Statistics and Actuarial
Science, Simon Fraser University, 8888 University Drive, Burnaby BC, Canada
V5A1S6. N. Epasinghege Dona is Senior Data Analyst, BC Health. Swartz has been
partially supported by the Natural Sciences and Engineering Research Council
of Canada. The work has been carried out with support from the CANSSI
(Canadian Statistical Sciences Institute) Collaborative Research Team (CRT) in
Sports Analytics.

###### Abstract

This paper considers the impact of unforced errors in sport. Although the
proposed methods are applicable to various sports, we demonstrate the approach
in the context of professional tennis. The value of the approach is that we
can provide estimates of the points lost, games lost, sets lost and matches
lost due to unforced errors. The methods are based on a bootstrapping
procedure which also yields standard errors for the estimates. The approach is
valuable in terms of player evaluation, and can also be used for training
purposes where it is possible to assess the quantification of improvement
based on fewer unforced errors.

Keywords: ATP and WTA tours, bootstrapping, tennis analytics.

##  1 INTRODUCTION

The term “unforced error” is often mentioned in sport. It is mistake made by a
player or by a team that is viewed as unnecessary and caused entirely by the
player or the team. For example, in American football, a coaching staff has
complete control over the assignment of players on the field, and therefore,
the penalty “too many men on the field” would be regarded as an unforced
error. Since unforced errors are controllable, they are regrettable, and are
often discussed and criticized.

Despite the widespread acknowledgment that unforced errors ought to be
minimized, the analysis of unforced errors is a topic that has received
minimal attention in the sports analytics literature. In badminton, Yadav and
Shukla (2011) considered 30 singles matches in university matches, and
observed that match winners had significantly fewer unforced errors than
losers. Again in badminton, Barreira and Chiminazzo (2020) considered 57 men’s
singles matches from the 2016 Olympic Games, and observed that match winners
had significantly fewer unforced errors than losers. In tennis, Brody (2006)
outlines how players can reduce unforced errors by adhering to various
principles of physics. Also in tennis, Attray and Attray (2021) provide a
theoretical framework for studying unforced errors. They suggest that the
minimum number of additional points needed to win a game given an unforced
error is a measure of the impact of an unforced error. Given varying match
circumstances (i.e., the score), unforced errors have impacts of varying
magnitude.

We suggest that there are two reasons for the paucity of sophisticated studies
involving unforced errors in matches:

  1. 1.

The determination of whether a mistake was an unforced error or a forced error
is often subjective, and typical datasets do not provide the distinction.

  2. 2.

It is not always obvious how to evaluate the impact of an unforced error. Had
the unforced error not occurred, what follows in the match is a
counterfactual. That is, without the unforced error, it is unclear how the
remainder of the match would have unfolded.

In this paper, we consider the impact of unforced errors in tennis.
Fortunately, we have access to a dataset that categorizes tennis errors as
either forced or unforced. The categorization is carried out by trained
volunteers who watch video recordings of matches. The training aspect is
important as it reduces subjectivity and improves the reliability of the
datasets. We overcome the counterfactual problem by simulating the remainder
of a match following an unforced error using a bootstrapping procedure. We
utilize the probabilities of what may have happened had the unforced error not
occurred. By repeatedly simulating the remainder of the match had the unforced
error not occurred, we obtain the expected impact of the unforced error.
Further, we can investigate “what-if” scenarios. For example, we might pose
the question, how would a player have fared in a match had the player reduced
unforced errors by 10%? Answers to these types of questions are valuable in
terms of player evaluation and player development.

Although unforced errors do not appear to have been analyzed in tennis, tennis
has seen an upsurge of work in analytics. One of the first serious statistical
contributions to tennis analytics was the work by Klassen and Magnus (2001)
which concerned an investigation of the iid assumption that points are
independent and identically distributed. They concluded that there is a
positive correlation between successive outcomes and that servers are less
likely to win a point in important situations.

Baker and McHale (2014, 2017) used the classical Bradley-Terry model to
compare the performance of tennis players. McHale and Morton (2011) employed
the same model for the purpose of forecasting tennis matches. Epasinghege
Dona, Gill and Swartz (2024) investigate player characteristics in relation to
to the length of the rally. A comprehensive text on tennis analytics for
Wimbledon is given by Klaasen and Magnus (2014).

Sports analytics has been assisted by the availability of tracking data. With
tracking data, player and ball locations are recorded with high frequency
(i.e., spatio-temporal data), and these detailed datasets have contributed to
explorations of many sporting problems that were previously unimaginable.
Gudmundsson and Horton (2017) provide a review on spatio-temporal analyses
used in invasion sports where player tracking data are available. In tennis,
there are a growing number of papers that provide statistical analyses that
rely on tracking data. Tea and Swartz (2023) investigate serving tendencies
which may allow a player to anticipate the nature of the opponent’s serve. The
approach relies on hierarchical models in a Bayesian framework. Kovalchik and
Albert (2022) also used a Bayesian framework where they investigate serve
returns by introducing a semiparametric mixture model. Other papers that use
tracking data but focus exclusively on the serve include Mecheri et al. (2016)
and Wei et al. (2015). The text by Albert et al. (2017) provides a flavour for
sports statistics across major sports.

In Section 2, we describe the data and associated issues. In Section 3, we
describe the proposed bootstrapping methods that permit the quantification of
the impact of unforced errors. Specifically, we provide estimates of the
expected points lost, games lost, sets lost and matches lost due to unforced
errors. In Section 4, we explore some descriptive statistics and we illustrate
the proposed methods in the context of a match between tennis professionals.
We conclude with a short discussion in Section 5.

##  2 TENNIS DATA

The data utilized in this paper are taken from the nearly 1.5 million serves
in the men’s (ATP) and women’s (WTA) professional tennis matches. The data are
available from the Match Charting Project
(https://github.com/JeffSackmann/tennis_MatchChartingProject) maintained by
Jeff Sackman. The data involve matches from 1970 through 2022 and contain
shot-by-shot outcomes involving 710 distinct men and 489 distinct women. This
public dataset provides information on shot type, shot direction, depth of
returns, types of errors, and more. The data were collected by trained
volunteers after watching the video recordings of matches. To the best of our
knowledge, there is no other source of publicly-available data of this type.
At present, these data appear to be underutilized with respect to statistical
analyses.

The data cover matches from all the major Grand Slam events, the Davis Cup and
many minor tournaments. The best players of this time period are included in
the dataset. Since the charting of matches was at the discretion of
volunteers, data are skewed in favour of later years. Also, better players
have a higher representation since they frequently reached the latter rounds
of tournaments.

###  2.1 Data Management Issues

When analyzing data, it is obviously important that the data are accurate.
Therefore, we carried out various procedures to check data accuracy.

In the MatchChartingProject dataset, rows correspond to points awarded.
Therefore, we augmented the dataset to include all serves. For example,
whenever a second serve occurred, this implied that there was a first serve
that resulted in a fault, for which no point was awarded.

We observed minimal numbers of obvious data entry errors. One example concerns
the rallyCount variable where there were five non-numeric characters in the
ATP data. We simply removed the corresponding rows when this occurred.

##  3 METHODS

Consider a match or a series of matches between player A and player B. We wish
to consider the impact of the unforced errors committed by player A in these
matches. Note that we only define unforced errors as those that lead directly
to a point obtained by player B; therefore first serve faults by player A are
not considered unforced errors.

We categorize the serves into four datasets according to:

  1. 1.

player A serving to player B on first serve

  2. 2.

player A serving to player B on second serve

  3. 3.

player B serving to player A on first serve

  4. 4.

player B serving to player A on second serve

We first consider the scenario where player A commits unforced errors at the
same rate that occurred in the datasets. We do not provide the details of
scoring in tennis but note that players accumulate points which lead to games
won which lead to sets won which lead to matches won. The bootstrapping
algorithm for this scenario is easily described and begins with one of the
players initiating service in Game 1 of Set 1. As the simulation of a match
progresses, the score of the match is updated and the server is determined
according to the rules of tennis. For each serve, the simulation of a serve
result proceeds as follows:

  * •

determine the relevant dataset based on the score of the match

  * •

randomly sample one serve from the dataset

    * –

record the outcome based on the sampled serve

The algorithm is simple and allows us to simulate a hypothetical match where
player A has made unforced errors at the same rate as in the datasets. We can
then repeat the algorithm over many matches to obtain the expected points won
by player A, the expected games won by player A, the expected sets won by
player A and the expected matches won by player A.

We then consider a second and unrealistic scenario where player A does not
make unforced errors. This permits the calculation of the same expected values
as before. Therefore, the differences in the expected values under the two
scenarios provides a quantification of the impact of unforced errors. The
second scenario of no unforced errors provides us with a baseline for the
maximum amount of improvement that a player can make.

The simulation algorithm in the second scenario is slightly more complicated
since it is unknown which player would have won the point had an unforced
error not occurred. For each serve in this scenario, the simulation of a serve
result proceeds as follows:

  * •

determine the relevant dataset based on the score of the match

  * •

randomly sample one serve from the dataset

    * –

if player A does not commit an unforced error, record the outcome based on the
sampled serve

    * –

if player A does commit an unforced error, record the outcome as generated
from the probabilities in Table 1

Table 1 provides player point probabilities derived from the work of
Epasinghege Dona, Gill and Swartz (2024). The probabilities correspond to A
winning the point had A not made an unforced error on touch $t$, and are
estimated from all serves in the ATP dataset. We define the serve as the first
touch $t=1$, and therefore unforced errors only occur in rallies corresponding
to touches $t=2,3,...$. From Table 1, we observe a monotonic increasing
pattern in the probabilities of the even number touches $t=2,4,6,8,10$. This
is logical since even number touches correspond to situations where A was the
receiver, and the service advantage to A dissipates as the rally lengthens.
Similarly, we observe a monotonic decreasing pattern in the probabilities of
the odd number touches $t=3,5,7,9$. This is also logical since odd number
touches correspond to situations where A was the server, and the service
disadvantage to A dissipates as the rally lengthens. We also note that the
probabilities converge to the common value 0.573 as the rally lengthens. This
also implies that the service advantage eventually disappears. It is
reasonable that the probabilities in Table 1 exceed 0.5 since A has avoided
the unforced error, and this puts B in a position where B must play the ball
and consequently may make a mistake.

Therefore, in the second simulation procedure where an unforced error is made,
we note the touch number where the unforced error was made, then assume that
the unforced error did not occur, and complete the outcome of the rally using
the probabilities in Table 1.

t | Probability
---|---
2 | 0.535
3 | 0.599
4 | 0.558
5 | 0.586
6 | 0.569
7 | 0.575
8 | 0.571
9 | 0.573
10 | 0.573
Table 1: Estimated probabilities of player A winning the point had player A
not made an unforced error on touch $t$. For $t>10$, we use the probabilities
corresponding to $t=10$.

The above algorithm allows us to quantify the impact of unforced errors in
terms of recognizable tennis statistics. However, we are also able to carry
out a what-if analysis for training purposes. Specifically, we can assess the
impact of player A reducing their unforced errors by a given probability, say
$x$. Here, we adjust the second bootstrapping algorithm where in a given serve
simulation, with probability $1-x$, we randomly select from the instances
where player A made an unforced error. In these cases, we replace the outcome
of the serve according to the probabilities in Table 1.

##  4 ANALYSIS

###  4.1 Descriptive Statistics

As the topic of unforced errors in tennis has not received much attention in
the academic literature, it is instructive to first explore our dataset using
some simple descriptive techniques. In what follows, we observe some features
of unforced errors that have not previously been investigated.

To begin, we consider all serves in the ATP and WTA datasets except for first
serve faults. Recall that we do not regard first serve faults as unforced
errors since they result in a second serve opportunity.

Accordingly, consider a rally that involves $t^{*}$ touches. The server would
have played the ball $t_{s}=(t^{*}+1)/2$ times if $t^{*}$ is odd, and the
server would have played the ball $t_{s}=t^{*}/2$ times if $t^{*}$ is even.
Similarly, the receiver would have played the ball $t_{r}=(t^{*}-1)/2$ times
if $t^{*}$ is odd, and the receiver would have played the ball $t_{r}=t^{*}/2$
times if $t^{*}$ is even. Therefore, we calculate a player’s unforced error
rate as the number of unforced errors made by the player divided by the number
of times that they played the ball. This is calculated over all opportunities
during a player’s professional career.

In Figure 1, we provide histograms for the unforced error rate for players in
the ATP tour (Figure 1a) and the WTP tour (Figure 1b). These histograms
consider only players who have played at least 10 matches. We observe that for
most players, the unforced error rates in tennis lie in the interval (5%,12%).
It seems that the unforced error rate is similar for men and women.

Figure 1: Histograms of unforced error rates for players on the ATP tour and
the WTA tour.

In Table 2, we provide the five players with the lowest unforced error rates
and the five players with the highest unforced error rates. This is done for
both the ATP and WTA tours. On the ATP tour, we observe that the five players
with the highest unforced error rates have not experienced great success. On
the other hand, we observe that five players with the lowest unforced error
rates include some famous players including Mats Wilander and Bjorn Borg. We
also note that Rafael Nadal has the sixth lowest unforced error rate of 5.4%.
This provides evidence concerning the importance of minimizing unforced
errors. Furthermore, amongst the 153 players in the ATP dataset, the 16
players with the most matches all had unforced error rates less than 8.2%. If
a player has played a lot of matches, then they have experienced sustained
success. This again emphasizes that low unforced error rates are related to
success. The same general patterns are observed on the WTA tour where
prominent players Chris Evert and Caroline Wozniacki appear on the low rate
list.

Tour | Player | Low Rates | Player | High Rates
---|---|---|---|---
ATP | Mats Wilander | 3.6 | Reilly Opelka | 12.3
ATP | Bjorn Borg | 4.9 | Aslan Karatsev | 12.1
ATP | Jaume Munar | 5.2 | Benoit Paire | 11.9
ATP | Gilles Simon | 5.2 | Peter Gojowczyk | 11.8
ATP | Daniel Evans | 5.3 | Nikoloz Basilashvili | 11.3
WTA | Sara Sorribes Tormo | 4.1 | Jelena Ostapenko | 13.1
WTA | Sara Errani | 4.8 | Madison Keys | 13.0
WTA | Monica Niculescu | 5.0 | Petra Kvitova | 12.6
WTA | Chris Evert | 5.3 | Julia Goerges | 12.3
WTA | Caroline Wozniacki | 5.3 | Dayana Yastremska | 11.9
Table 2: Unforced error rate probabilities for players on the ATP and WTA
tours.

From Figure 1, we see that on average, the unforced error rate across players
is roughly 8.0% (ATP) and 8.6% (WTA). As significant as this statistic may
appear, unforced errors may be viewed as even more detrimental to success in
tennis. Recall, every point must terminate with either a forced or an unforced
error. When we look at the number of serves in tennis that end by an unforced
error, the rates are 30.3% and 32.9% on the ATP and WTA tours, respectively.

We know that there are different issues in tennis with respect to the touch
number. For example, with big servers, touch $t=2$ is a great challenge for
the receiver as the main focus is to simply return the serve. In Figure 2, we
plot the unforced error rates for both the ATP and WTA tours according to
touch number. Given that longer rallies are less common, we have truncated the
plots. We have also separated the plots according to touches $t=3,5,\ldots,13$
for the server and touches $t=2,4,\ldots,12$ for the receiver. We generally
observe that the server’s unforced error rate is steady across the touch
number $t$. On the other hand, the receiver’s unforced error rate is lower for
the initial return $t=2$ and then appears to stabilize. This may be because
the receiver’s main objective with respect to the serve is to simply return
the serve instead of hitting a winning shot. We observe similar patterns in
both the men’s and women’s tours.

(a) ATP tour

(b) WTA tour

Figure 2: Unforced error rates on the ATP and WTA tours with respect to touch
number.

It is also interesting to investigate how the game has changed with respect to
unforced errors. In Figure 3, we plot the unforced error rates in the ATP and
WTA tours versus year. We observe that the unforced error rate has gradually
increased over time in both the men’s and women’s games. We do not have a
definitive explanation for this. One conjecture relates to technology changes
in tennis, especially relating to racquets. With such advances, there has been
an increase in the speed of the game which in turn may make it more difficult
to keep the ball in play, and consequently increase unforced errors.

Figure 3: Unforced error rates on the ATP and WTA tours over time.

###  4.2 Bootstrapping Investigation

We now illustrate the proposed methods by considering two prominent players,
Roger Federer and Rafael Nadal. They have enjoyed an intense rivalry,
beginning in 2004, and have played against one another in 9 Grand Slam finals.
In Table 3, we provide some descriptive comparisons of the two players. In
comparison with Figure 1, we note that Federer’s career unforced error rate of
8.2% is about average for an ATP player whereas Nadal’s rate of 5.4% is
extremely low. We note that Nadal has made tentative plans to retire in 2025.

Player | ATP Start (Age) | ATP End | ATP Matches | UFE
---|---|---|---|---
Roger Federer | 1998 (17) | 2022 | 1526 | 0.082
Rafael Nadal | 2001 (14) | 2024-active | 1300+ | 0.054
Table 3: Information concerning the singles careers of Federer and Nadal with
career unforced error rates (UFE).

We now carry out an analysis involving Federer and Nadal with the results
shown in Table 4. We show actual results between the two players, simulated
results using Federer’s historic unforced error rate, simulated results with
Federer’s unforced error rate reduced by 10% and simulated results with
Federer’s unforced errors eliminated (i.e. reduced by 100%). These results are
based on the 35 ATP matches recorded between Federer and Nadal in our dataset.
We note that official records indicate that these two players actually faced
each other 40 times; hence the dataset is missing 5 matches. In the
simulations, we do not use the career unforced error rates. Instead, we use
the unforced error rates (UFE) exhibited between the two players in the 35
matches. In these matches, Federer’s UFE was 8.9% and Nadal’s UFE was 5.7%.
The inter-personal UFE rates between the two players are very close to their
career rates. Also, for simulated match results, a win is defined as the first
player to have won three sets. The results in Table 4 are based on simulations
from $N=3000$ matches that took roughly two hours on a laptop computer.

From Table 4, in actual competitions, we observe that Federer won 49.3% of the
points against Nadal, demonstrating the competitiveness of the matches.
However, at the match level, Federer won only 16 of the 40 matches, i.e. only
40% of the time. This demonstrates something about Nadal’s ability to come
through in clutch situations. This is further confirmed by Federer’s
simulation match results which are better than his actual match results. We
have confidence in the simulation results since the differences between
Federer’s actual points won and simulated points won is minor. Most
importantly, we quantify the extent to which Federer’s unforced errors
impacted his results against Nadal. Had Federer been able to reduce his
unforced error rate by even a small margin (10%), his match results would have
increased substantially to 53.7%, placing him above Nadal. And as an upper
bound for improvement, we consider the simulation case where Federer
eliminated all unforced errors. In this unrealistic scenario, Federer would
have been completely dominant. This exercise illustrates the critical impact
of unforced errors in tennis.

Cases | Points Won | Games Won | Matches Won
---|---|---|---
Federer - Actual | 49.3 | 47.7 | 40.0
Federer - Simulated Historic UFE Rate | 49.5 (0.07) | 48.9 (0.18) | 44.5 (0.90)
Federer - Simulated 10% UFE Reduction | 50.3 (0.07) | 50.7 (0.18) | 53.7 (0.91)
Federer - Simulated 100% UFE Reduction | 58.6 (0.08) | 67.3 (0.20) | 98.7 (0.21)
Table 4: Federer percentage results against Nadal under various scenarios.
Bootstrap standard errors are given in parentheses.

##  5 DISCUSSION

Whereas unforced errors in tennis are a much discussed topic, the
quantification of the impact of unforced errors is something that appears to
be missing from the literature. This paper uses bootstrapping methods to
address this novel problem. The consequent understanding of a player’s
unforced error rate can then be used for player evaluation and training.

In Section 4.2, we investigated the the impact of unforced errors between two
giants of men’s tennis, Roger Federer and Rafael Nadal. This was a convenient
investigation since these two players had an extensive history of playing
against one another, and hence, their unforced error rates could be well-
assessed. In future investigations, one may wish to study pairs of players
with less of a competitive history. In these cases, one could use their
performances against the field (i.e. the typical player) to assess unforced
error rates. Doing so would allow us to compare any pair of players, even
those from different eras.

Another avenue for future investigation may involve digging deeper in the
circumstances in which unforced errors are committed. In live matches,
unforced error rates are routinely reported, but knowing when these errors
occur is relevant to performance. For example, they may occur on high or low
leverage points, early or late in a rally, etc.

##  6 REFERENCES

Albert, J.A., Glickman, M.E., Swartz, T.B. and Koning, R.H., Editors (2017).
Handbook of Statistical Methods and Analyses in Sports, Chapman & Hall/CRC
Handbooks of Modern Statistical Methods, Boca Raton.

Attray, H. and Attray, S. (2021). Cost of an unforced error in tennis - A
statistical approach. International Tennis Federation, Issue 84, 12-14.

Baker, R.D. and McHale, I. (2014). A dynamic paired comparisons model: Who is
the greatest tennis player? European Journal of Operational Research, 236 (2),
677-684.

Baker, R.D. and McHale, I. (2017). An empirical Bayes model for time-varying
paired comparisons rankings: Who is the greatist women’s tennis player?
European Journal of Operational Research, 258 (1), 328-333.

Barreira, J. and Chiminazzo, J.G.C. (2020). Who, how and when to perform
winner points and unforced errors in badminton matches? An analysis of men’s
single matches in the 2016 Olympic Games. International Journal of Performance
Analysis in Sport, 20(4), 610-619.

Brody, H. (2006). Unforced errors and error reduction in tennis. British
Journal of Sports Medicine, 40(5), 397-400.

Epasinghege Dona, N., Gill, P.S. and Swartz, T.B. (2024). What does rally
length tell us about player characteristics in tennis? Journal of the Royal
Statistical Society, Series A, https://doi.org/10.1093/jrssa/qnae027

Gudmundsson, J. and Horton, M. (2017). Spatio-temporal analysis of team
sports. ACM Computing Surveys, 50(2), Article 22.

Klassen, F.J.G.M. and Magnus, J.R. (2001). Are points in tennis independent
and identically distributed? Evidence from a dynamic binary panel model.
Journal of the American Statistical Association, 96(454), 500-509.

Klassen, F. and Magnus, J.R. (2014). Analyzing Wimbledon: The Power of
Statistics, University Press Scholarship Online,
https://doi.org/10.1093/acprof:oso/9780199355952.001.0001

Kovalchik, S.A. and Albert, J. (2022). A statistical model of serve return
impact patterns in professional tennis. Accessed October 13, 2022 at
https://arxiv.org/abs/2202.00583

McHale, I. and Morton, A. (2011). A Bradley-Terry type model for forecasting
tennis match results. International Journal of Forecasting, 27 (2), 619-630.

Mecheri, S., Rioult, F., Mantel, B., Kauffmann, F. and Benguigui, N. (2016).
The serve impact in tennis: First large-scale study of big hawk-eye data.
Statistical Analysis and Data Mining, 9(5), 310-325.

Tea, P. and Swartz, T.B. (2022). The analysis of serve decisions in tennis
using Bayesian hierarchical models. Annals of Operations Research, 325(1),
633-648.

Wei, X., Lucey, P., Morgan, S., Carr, P., Reid, M. and Sridharan, S. (2015).
Predicting serves in tennis using style priors. Proceedings of the 21th ACM
SIGKDD International Conference on Discovery and Data Mining, 2207-2215.

Yadav, S.K. and Shukla, Y.M. (2011). Analysis of unforced errors in relation
to performance in singles in badminton. International Journal of Physical
Education, 4(2), 117-119.


# Explicit high-order generalized-$\alpha$ methods for isogeometric analysis
of structural dynamics

Pouria Behnoudfar<EMAIL_ADDRESS>Gabriele Loli
Alessandro Reali  Giancarlo Sangalli  Victor M. Calo  Curtin Institute for
Computation & School of Earth and Planetary Sciences, Curtin University, Kent
Street, Bentley, Perth, WA 6102, Australia  Dipartimento di Matematica “F.
Casorati”, Università di Pavia, 27100 Pavia, Italy  Department of Civil
Engineering and Architecture, Università di Pavia, 27100 Pavia, Italy  School
of Electrical Engineering, Computing and Mathematical Sciences, Curtin
University, P.O. Box U1987, Perth, WA 6845, Australia

###### Abstract

We propose a new family of high-order explicit generalized-$\alpha$ methods
for hyperbolic problems with the feature of dissipation control. Our approach
delivers $2k,\,\left(k\in\mathbb{N}\right)$ accuracy order in time by solving
$k$ matrix systems explicitly and updating the other $2k$ variables at each
time-step. The user can control the numerical dissipation in the discrete
spectrum’s high-frequency regions by adjusting the method’s coefficients. We
study the method’s spectrum behaviour and show that the CFL condition is
independent of the accuracy order. The stability region remains invariant
while we increase the accuracy order. Next, we exploit efficient
preconditioners for the isogeometric matrix to minimize the computational
cost. These preconditioners use a diagonal-scaled Kronecker product of
univariate parametric mass matrices; they have a robust performance with
respect to the spline degree and the mesh size, and their decomposition
structure implies that their application is faster than a matrix-vector
product involving the fully-assembled mass matrix. Our high-order schemes
require simple modifications of the available implementations of the
generalized-$\alpha$ method. Finally, we present numerical examples
demonstrating the methodology’s performance regarding single- and multi-patch
IGA discretizations.

###### keywords:

Generalized-$\alpha$ method, high-order time integrator, Explicit method,
dissipation control, CFL condition, hyperbolic, preconditioner, Isogeometric
analysis

††journal: Arxiv

##  1 Introduction

Chung and Hulbert [1] introduced the generalized-$\alpha$ method for
hyperbolic systems arising in structural dynamics; later, Hulbert and Chung
[2] presented an explicit version of this time-marching method. The method is
second-order accurate in time and provides user-control on the high-frequency
numerical dissipations. Likewise, other well-known methods like, e.g., the
Newmark-$\beta$ [3] and the HHT-$\alpha$ [4] methods have both implicit and
explicit formulations and are limited to the second-order accuracy.
Nevertheless, the Newmark method does not control the dissipation while, using
the HHT-$\alpha$, one obtains dissipative solutions in low-frequency regions.
The generalized-$\alpha$ method produces an algorithm that combines high-
frequency and low-frequency dissipations optimally; see [1, 2].

These methods, including the generalized-$\alpha$ method, are limited to the
second-order accuracy in time. In contrast, high-order methods like Lax-
Wendroff, Runge-Kutta, Adams-Moulton, and backward differentiation schemes
(see [5]) or high-order IGA collocation methods for dynamics (see, e.g., [6])
lack explicit control over the numerical dissipation of high frequencies.
Thus, we propose an explicit $k$-step generalized-$\alpha$ method that
delivers $2k$ accuracy in time for second derivatives problems in time. Our
method obtains high-order accuracy in the time with optimal control of the
resulting system’s spectral behaviour. Thus, we build an algorithm that
consists of $3k$ equations; for each three-equation set, we solve an explicit
system for a variable and update the other two. We then study the resulting
amplification matrix’s spectral properties to determine the CFL condition and
introduce user-defined parameters that control the numerical dissipation.
Later, we prove that the CFL condition of our method is independent of
temporal accuracy.

The explicit generalized-$\alpha$ method, at each time step, factors the
isogeometric mass matrix; thus, we use a preconditioned conjugate gradients
(PCG) as an iterative solver using $\boldsymbol{\mathcal{M}}$ [7], which is
easy to implement, extremely efficient and robust. We build the preconditioner
using a diagonal scaling of the parametric mass matrix (i.e., matrix
associated with the $d$-dimensional cube pre-image of the parametric space)
for single-patch geometries. We combine the preconditioners defined above for
each patch for multi-patch geometries using an additive Schwarz domain-
decomposition method. In [7], the authors proved that these preconditioners
are efficient to apply and are robust with respect to the mesh size; in the
single-patch case, they show that the behaviour improves as the problem size
grows (mesh refinement).

We present numerical simulations to demonstrate the performance of the time-
marching scheme and its preconditioner. We provide numerical evidence on the
time integrator’s high-order accuracy and its dispersion properties.
Additionally, we show the optimal convergence in the spatial domain for
several problems. The outline for the remainder of the paper is as follows.
Section 2 describes the hyperbolic problem we consider and introduces the
spatial discretizations to obtain the matrix formulation of the problem.
Section 3 presents our fourth-order explicit generalized-$\alpha$ method;
therein, we also analyze the method’s stability, its CFL condition, and its
temporal accuracy. Section 4 generalizes the method to $2k^{th}$-order of
accuracy. Then, we introduce our solver in Section 5 with details on single-
and multi-patch isogeometric analysis. We verify the solver’s convergence and
its computational performance in Section 6 numerically. Section 7 describes
our contributions and further applications.

##  2 Problem Statement

We start with an initial boundary-value hyperbolic problem, a model problem
for structural dynamics:

$\displaystyle\begin{cases}\ddot{u}(\boldsymbol{x},t)-\nabla\cdot(\omega^{2}\,\nabla
u(\boldsymbol{x},t))+c\left(\dot{u}(\boldsymbol{x},t)\right)&=f(\boldsymbol{x},t),\qquad\,\,\,(\boldsymbol{x},t)\in\Omega\times[0,T],\\\
u(\boldsymbol{x},t)&=u_{D},\qquad\qquad\boldsymbol{x}\in\partial\Omega,\\\
u(\boldsymbol{x},0)&=u_{0},\qquad\,\qquad\boldsymbol{x}\in\Omega,\\\
\dot{u}(\boldsymbol{x},0)&=v_{0},\qquad\,\,\qquad\boldsymbol{x}\in\Omega,\end{cases}$
(1)

Let $\Omega\subset\mathbb{R}^{d},d=1,2,3,$ be an open bounded domain. The
operator $\nabla$ is the spatial gradient and a superscript dot denotes a time
derivative such that
$\dot{u}(\boldsymbol{x},t)=\frac{\partial{u}(\boldsymbol{x},t)}{\partial t}$
and
$\ddot{u}(\boldsymbol{x},t)=\frac{\partial^{2}{u}(\boldsymbol{x},t)}{\partial
t^{2}}$. $c\left(\dot{u}(\boldsymbol{x},t)\right)$ models linear damping. The
source $f$, propagation speed $\omega$, initial data $u_{0},\,v_{0}$, and
Dirichlet boundary conditions $u_{D}$ are given and assumed regular enough for
the problem to admit a weak solution. In order to derive our numerical method
for (1), we first obtain a semi-discretized problem by discretizing in space,
then, we deploy our explicit generalized-$\alpha$ method to have a fully
discretized system.

###  2.1 Spatial discretization

Adopting a Galerkin method (in particular, isogeometric analysis), the matrix
problem resulting from the semi-discretization of (1) reads:

$M\ddot{U}+C\dot{U}+KU=F,$ (2)

where $M$, $C$, and $K$ are the mass, damping, and stiffness matrices,
respectively. $U$ denotes the vector of the unknowns, and $F$ is the source
vector. The initial conditions also read:

$U(0)=U_{0},\qquad V(0)=V_{0},$ (3)

where $U_{0}$ and $V_{0}$ represent the given vectors initial conditions
corresponding to $u_{0,h}$ and $\dot{u}_{0,h}$, respectively. In the next
section, we propose our numerical technique to deal with the time derivative
$\ddot{U}$ and $\dot{U}$ in (2) with the accuracy of order $2k$ in the
temporal domain with $k\in\mathbb{N}$.

######  Remark 1.

Herein, we propose an explicit generalized-$\alpha$ scheme by considering a
general spatial discretization, leading to the matrix problem (2). Therefore,
the use of isogeometric analysis does not limit the method’s applicability;
effectively, our method applies to the spatial discretization of a time-
dependent semi-discretized problem.

######  Remark 2.

In problem (1) and accordingly (2), for simplicity, we only consider constant
$\omega$ and assume that the solution $u(\cdot\,,t)$ satisfies homogeneous
boundary condition. One requires slight modifications of the discrete bilinear
and linear functions for the cases of heterogeneous propagation speed and non-
homogeneous boundary conditions [8, 9].

###  2.2 Time-discretization

To obtain a fully discrete problem (2), we adopt an appropriate time marching
scheme to deal with $\ddot{u}_{h}$ and $\dot{u}_{h}$; in the next section, we
propose a new high-order explicit generalized-$\alpha$ method.

##  3 Explicit generalized-$\alpha$ method

Consider a partitioning of the time interval $[0,T]$ as
$0=t_{0}<t_{1}<\cdots<t_{N}=T$ with a grid size $\tau_{n}=t_{n+1}-t_{n}$. We
approximate $U(t_{n}),\,\dot{U}(t_{n}),\,\ddot{U}(t_{n})$ using
$U_{n},\,V_{n},\,A_{n}$, respectively. Tthe explicit generalized-$\alpha$
method with second-order accuracy in time solves (2) at time-step $t_{n+1}$,

$\displaystyle MA_{n+\alpha_{m}}+CV_{n}+KU_{n}$
$\displaystyle=F_{n+\alpha_{f}},$ (4a) $\displaystyle V_{n+1}$
$\displaystyle=V_{n}+\tau A_{n}+\tau\gamma\llbracket A_{n}\rrbracket,$ (4b)
$\displaystyle U_{n+1}$ $\displaystyle=U_{n}+\tau
V_{n}+\frac{\tau^{2}}{2}A_{n}+\tau^{2}\beta\llbracket A_{n}\rrbracket,$ (4c)

where

$\displaystyle F_{n+\alpha_{f}}$ $\displaystyle=F(t_{n+\alpha_{f}}),$ (5)
$\displaystyle A_{n+\alpha_{m}}$ $\displaystyle=A_{n}+\alpha_{m}\llbracket
A_{n}\rrbracket$ $\displaystyle\llbracket A_{n}\rrbracket$
$\displaystyle=A_{n+1}-A_{n}.$

To retrieve the unknown $A$ at the initial state, we solve

$A_{0}=M^{-1}(F_{0}-KU_{0}-CV_{0}).$ (6)

At each time step, we first compute $A_{n+1}$ using (4a), and then evaluate
$V_{n+1}$ and $U_{n+1}$ from (4b) and (4c), respectively. The method (4) has a
truncation error in time of $\mathcal{O}(\tau^{3})$ (for further details, see
[10, 11]). We extend the accuracy, assuming sufficient temporal regularity of
the solution; thus, we use a Taylor expansion with higher-order terms. We
introduce new variables $\mathcal{L}^{a}(A_{n})$ to approximate the $a$-th
order derivative of $A_{n}$ in time. Therefore, for example, to derive a
fourth-order explicit generalized-$\alpha$ method, we solve

$\displaystyle MA_{n}^{\alpha_{1}}+CV_{n}+KU_{n}$
$\displaystyle=F_{n+\alpha_{f1}},$ (7) $\displaystyle
M\mathcal{L}^{3}(A_{n})^{\alpha_{2}}+C\mathcal{L}^{2}(A_{n})+K\mathcal{L}^{1}(A_{n})$
$\displaystyle=F^{(3)}_{n+\alpha_{f2}},$

with updating conditions

$\displaystyle U_{n+1}$ $\displaystyle=U_{n}+\tau
V_{n}+\frac{\tau^{2}}{2}A_{n}+\frac{\tau^{3}}{6}\mathcal{L}^{1}(A_{n})+\frac{\tau^{4}}{24}\mathcal{L}^{2}(A_{n})+\frac{\tau^{5}}{120}\mathcal{L}^{3}(A_{n})+\beta_{1}\tau^{2}P_{n},$
(8) $\displaystyle V_{n+1}$ $\displaystyle=V_{n}+\tau
A_{n}+\frac{\tau^{2}}{2}\mathcal{L}^{1}(A_{n})+\frac{\tau^{3}}{6}\mathcal{L}^{2}(A_{n})+\frac{\tau^{4}}{24}\mathcal{L}^{3}(A_{n})+\gamma_{1}\tau
P_{n},$ $\displaystyle\mathcal{L}^{1}(A_{n+1})$
$\displaystyle=\mathcal{L}^{1}(A_{n})+\tau\mathcal{L}^{2}(A_{n})+\frac{\tau^{2}}{2}\mathcal{L}^{3}(A_{n})+\tau^{2}\beta_{2}\llbracket\mathcal{L}^{3}(A_{n})\rrbracket,$
$\displaystyle\mathcal{L}^{2}(A_{n+1})$
$\displaystyle=\mathcal{L}^{2}(A_{n})+\tau\mathcal{L}^{3}(A_{n})+\tau\gamma_{2}\llbracket\mathcal{L}^{3}(A_{n})\rrbracket,$

where

$\displaystyle P_{n}$
$\displaystyle=A_{n+1}-A_{n}-\tau\mathcal{L}^{1}(A_{n})-\frac{\tau^{2}}{2}\mathcal{L}^{2}(A_{n})-\frac{\tau^{3}}{6}\mathcal{L}^{3}(A_{n}),$
(9) $\displaystyle A_{n}^{\alpha_{1}}$
$\displaystyle=A_{n}+\tau\mathcal{L}^{1}(A_{n})+\frac{\tau^{2}}{2}\mathcal{L}^{2}(A_{n})+\frac{\tau^{3}}{6}\mathcal{L}^{3}(A_{n})+\alpha_{1}P_{n},$
$\displaystyle\mathcal{L}^{3}(A_{n})^{\alpha_{2}}$
$\displaystyle=\mathcal{L}^{3}(A_{n})+\alpha_{2}\llbracket\mathcal{L}^{3}(A_{n})\rrbracket.$

We approximate $\frac{\partial^{(a)}}{\partial
t^{(a)}}f(\cdot\,,\,n+\alpha_{f2})$ as $F^{(a)}_{n+\alpha_{f2}}$. Following
(6), one can readily obtain the initial data of the unknowns using the given
information on $U_{0}$ and $V_{0}$ as

$\displaystyle A_{0}$ $\displaystyle=M^{-1}\left(F_{0}-KU_{0}-CV_{0}\right),$
(10) $\displaystyle\mathcal{L}^{1}(A_{0})$
$\displaystyle=M^{-1}\left(F_{0}^{(1)}-KV_{0}-CA_{0}\right),$
$\displaystyle\mathcal{L}^{2}(A_{0})$
$\displaystyle=M^{-1}\left(F_{0}^{(2)}-KA_{0}-C\mathcal{L}^{1}(A_{0})\right),$
$\displaystyle\mathcal{L}^{3}(A_{0})$
$\displaystyle=M^{-1}\left(F_{0}^{(3)}-K\mathcal{L}^{1}(A_{0})-C\mathcal{L}^{2}(A_{0})\right).$

Next, we derive the corresponding coefficients that deliver the desired
accuracy and discuss the method’s stability.

###  3.1  Order of accuracy in time

We now determine the parameters $\gamma_{1}$ and $\gamma_{2}$ that guarantee
the fourth-order accuracy of (7); thus, we use the following result.

######  Theorem 1.

Assuming that the solution is sufficiently smooth with respect to time, the
method (7) is fourth-order accurate in time given

$\gamma_{1}=\frac{1}{2}-\alpha_{f1}+\alpha_{1},\qquad\qquad\gamma_{2}=\frac{1}{2}-\alpha_{f2}+\alpha_{2}.$
(11)

###### Proof.

Following [12], we determine the amplification matrix by associating with
$\lambda=\theta^{2}$ the eigenvalues of the matrix $M^{-1}K$. The damping
term’s eigenvalues $M^{-1}C$ are $\zeta\theta$ where $\zeta$ is the damping
coefficient. Without loss of the proof’s generality, we set $\zeta=0$ and
ignore damping. Then, substituting (8) into (7), we obtain a system of
equations at each time step as

$A\mathbb{U}_{n+1}=B\mathbb{U}_{n}+\mathbb{F}_{n+\alpha_{f}},$ (12)

letting $\mathbb{U}_{n}^{T}=\begin{bmatrix}U_{n},\ \tau V_{n},\
\tau^{2}A_{n},\ \tau^{3}\mathcal{L}^{1}(A_{n}),\
\tau^{4}\mathcal{L}^{2}(A_{n}),\
\tau^{5}\mathcal{L}^{3}(A_{n})\end{bmatrix}^{T}$,

$\displaystyle A$ $\displaystyle=\begin{bmatrix}1&0&-\beta_{1}&0&0&0\\\
0&1&-\gamma_{1}&0&0&0\\\ 0&0&\alpha_{1}&0&0&0\\\ 0&0&0&1&0&-\beta_{2}\\\
0&0&0&0&1&-\gamma_{2}\\\ 0&0&0&0&0&\alpha_{2}\end{bmatrix},$ (13)
$\displaystyle B$
$\displaystyle=\begin{bmatrix}1&1&\frac{1}{2}-\beta_{1}&\frac{1}{6}-\beta_{1}&\frac{1}{24}-\frac{\beta_{1}}{2}&\frac{1}{120}-\frac{\beta_{1}}{6}\\\
0&1&1-\gamma_{1}&\frac{1}{2}-\beta_{1}&\frac{1}{6}-\frac{\beta_{1}}{2}&\frac{1}{24}-\frac{\beta_{1}}{6}\\\
-\tau^{2}\lambda&-\tau^{2}\lambda&\alpha_{1}-1-\frac{1}{2}\tau^{2}\lambda&\alpha_{1}-1-\frac{1}{6}\tau^{2}\lambda&\frac{1}{2}(\alpha_{1}-1)&\frac{1}{6}(\alpha_{1}-1)\\\
0&0&0&1&1&\frac{1}{2}-\beta_{2}\\\ 0&0&0&0&1&1-\gamma_{2}\\\
0&0&0&-\tau^{2}\lambda&0&\alpha_{2}-1\end{bmatrix}$

and $\mathbb{F_{n+\alpha_{f}}}$ consists of the forcing terms
$F_{n+\alpha_{f1}}$ and $(F^{(3)}_{n+\alpha_{f2}})$. Thus, the amplification
matrix $G$ becomes

$G=A^{-1}B,$ (14)

with $G$ being an upper-block triangular matrix written as:

$G=\begin{bmatrix}\Lambda_{1}&\Xi\\\ \boldmath{0}&\Lambda_{2}\end{bmatrix},$
(15)

where

$\displaystyle\Lambda_{1}$
$\displaystyle=\dfrac{1}{\alpha_{1}}\begin{bmatrix}\alpha_{1}-{\beta_{1}\tau^{2}\lambda}&\alpha_{1}-{\beta_{1}\tau^{2}\lambda}&\frac{1}{2}\left({\alpha_{1}-\beta_{1}\left(\tau^{2}\lambda+2\right)}\right)\\\
-{\gamma_{1}\tau^{2}\lambda}&\alpha_{1}-{\gamma_{1}\tau^{2}\lambda}&\alpha_{1}-\frac{\gamma_{1}}{2}\left(\tau^{2}\lambda+2\right)\\\
-{\tau^{2}\lambda}&-{\tau^{2}\lambda}&\alpha_{1}-1-\frac{\tau^{2}}{2}\lambda\\\
\end{bmatrix},$ (16) $\displaystyle\Lambda_{2}$
$\displaystyle=\dfrac{1}{\alpha_{2}}\begin{bmatrix}\alpha_{2}-{\beta_{2}\tau^{2}\lambda}&\alpha_{2}&\frac{\alpha_{2}}{2}-{\beta_{2}}\\\
-{\gamma_{2}\tau^{2}\lambda}&\alpha_{2}&\alpha_{2}-{\gamma_{2}}\\\
-{\tau^{2}\lambda}&0&{\alpha_{2}-1}\end{bmatrix}.$

Following [11, 13], firstly, we obtain the high-order unknowns
$\mathcal{L}^{1}(A_{n})$ and $\mathcal{L}^{2}(A_{n})$, being second-order
accurate. Secondly, we set the parameters such that the upper block on the
diagonal also delivers second-order accuracy. Then, we obtain the fourth-order
accuracy by adding the high-order terms through the upper off-diagonal block
to our solution. Therefore, we study each diagonal block to derive the related
parameters. For this aim, for any arbitrary amplification matrix, we can state
the following,

$G_{0}\mathcal{L}^{1}(A_{n+1})-G_{1}\mathcal{L}^{1}(A_{n})+G_{2}\mathcal{L}^{1}(A_{n-1})-G_{3}\mathcal{L}^{1}(A_{n-2})=0,$
(17)

where the coefficients are invariants of the amplification matrix as
$G_{0}=1$, $G_{1}$ is the trace of $G$, $G_{2}$ is the sum of principal minors
of $G$, and $G_{3}$ is the determinant of $G$. Using a Taylor series
expansion, we have:

$\displaystyle\mathcal{L}^{1}(A_{n+1})$
$\displaystyle=\mathcal{L}^{1}(A_{n})+\tau\mathcal{L}^{2}(A_{n})+\frac{\tau^{2}}{2}\mathcal{L}^{3}(A_{n})+\mathcal{O}(\tau^{3}),$
(18) $\displaystyle\mathcal{L}^{1}(A_{n-1})$
$\displaystyle=\mathcal{L}^{1}(A_{n})-\tau\mathcal{L}^{2}(A_{n})+\frac{\tau^{2}}{2}\mathcal{L}^{3}(A_{n})+\mathcal{O}(\tau^{3}),$
$\displaystyle\mathcal{L}^{1}(A_{n-2})$
$\displaystyle=\mathcal{L}^{1}(A_{n})-2\tau\mathcal{L}^{2}(A_{n})+{2\tau^{2}}\mathcal{L}^{3}(A_{n})+\mathcal{O}(\tau^{3}).$

Setting $\gamma_{2}=\frac{1}{2}-\alpha_{f2}+\alpha_{2}$, we obtain that
$\mathcal{L}^{1}(A_{n+1}),\mathcal{L}^{2}(A_{n+1})$ are second-order accurate
in time. Next, we rewrite (18) for the unknown $U_{n+1}$ as

$\displaystyle U_{n+1}$ $\displaystyle=U_{n}+\tau
V_{n}+\frac{\tau^{2}}{2}A_{n}+\mathcal{R},$ (19) $\displaystyle U_{n-1}$
$\displaystyle=U_{n}-\tau V_{n}+\frac{\tau^{2}}{2}A_{n}+\mathcal{R},$
$\displaystyle U_{n-2}$ $\displaystyle=U_{n}-2\tau
V_{n}+4\frac{\tau^{2}}{2}A_{n}+\mathcal{R},$

where $\mathcal{R}$ is a function of
$\mathcal{L}^{1}(A_{n}),\mathcal{L}^{2}(A_{n})$ and $\mathcal{L}^{3}(A_{n})$.
We neglect this term in the analysis as it is a residual. Thus, we can prove
that the remaining terms are second-order accurate in time. Then, we add the
residuals to the second-order accurate solution, to produce a truncation error
of $\mathcal{O}(\tau^{5})$ and consequently, we obtain a fourth-order accurate
scheme in time, which completes the proof. ∎

###  3.2 Stability analysis and CFL condition

We bound the spectral radius of the amplification matrix by one to deliver a
stable time-marching scheme. For this, firstly, we calculate the eigenvalues
of the matrix $G$ in (14) for the case $\Theta:=\tau^{2}\lambda\to 0$.
Therefore, in this case, the diagonal blocks in (16) are:

$\displaystyle\Lambda_{1}$
$\displaystyle=\begin{bmatrix}1&1&\frac{1}{2}-\frac{\beta_{1}}{\alpha_{1}}\\\
0&1&1-\frac{\gamma_{1}}{\alpha_{1}}\\\
0&0&1-\frac{1}{\alpha_{1}}\end{bmatrix},\ \Lambda_{2}$
$\displaystyle=\begin{bmatrix}1&1&\frac{1}{2}-\frac{\beta_{2}}{\alpha_{2}}\\\
0&1&1-\frac{\gamma_{2}}{\alpha_{2}}\\\
0&0&1-\frac{1}{\alpha_{2}}\end{bmatrix}.$ (20)

Then, $G$’s eigenvalues are:

$\lambda_{1}=\lambda_{2}=\lambda_{3}=\lambda_{4}=1,\qquad\lambda_{5}=1-\frac{1}{\alpha_{1}},\qquad\lambda_{6}=1-\frac{1}{\alpha_{2}}.$
(21)

The boundedness of $\lambda_{5}$ and $\lambda_{6}$ in (21) implies that
${\alpha_{1}}\geq\frac{1}{2}$ and ${\alpha_{2}}\geq\frac{1}{2}$.

Now, we derive the CFL condition for our method. In the analysis of the
implicit generalized-$\alpha$ methods of second- and higher-order accuracy [1,
11], we analyze the discrete system’s eigenvalue distribution in the limit
$\Theta\to\infty$ and set the method’s free parameters such that all
eigenvalues are equal to a real constant $\rho_{\infty}\in[0,1]$; this
parameter controls the numerical dissipation. In our explicit method, we need
to find the method’s conditional stability region. Thus, we exploit the
amplification matrix’s upper-triangular structure (15) to explicitly compute
the eigenvalues of $G$, which is equivalent to finding the eigenvalues of each
diagonal block in (16). Therefore, the characteristic polynomials for the
diagonal blocks are:

$\displaystyle\Lambda_{1}:\qquad 2\alpha_{1}\lambda_{1}^{3}$
$\displaystyle+{\lambda_{1}^{2}\left(\Theta(2\beta_{1}+2\gamma_{1}+1)-6\alpha_{1}+2\right)}$
(22)
$\displaystyle+{\lambda_{1}\left(\Theta(-4\beta_{1}-2\gamma_{1}+1)+6\alpha_{1}-4\right)}+{2\beta_{1}\Theta-2\alpha_{1}+2}=0,$
$\displaystyle\Lambda_{2}:\qquad 2\alpha_{2}\lambda_{2}^{3}$
$\displaystyle+{\lambda_{2}^{2}\left(2\Theta\beta_{2}-6\alpha_{2}+2\right)}+{\lambda_{2}\left(\Theta(-4\beta_{2}+2\gamma_{2}+1)+6\alpha_{2}-4\right)}$
$\displaystyle+{\Theta(2\beta_{2}+1-2\gamma_{2})-2\alpha_{2}+2}=0,$

where two roots of each characteristic polynomial in (22) are the principal
roots while the third is spurious (unphysical). We require the two principal
roots to be complex conjugates except in the high-frequency regions. This
requirement maximizes high-frequency dissipation while setting two eigenvalues
to one in the low-frequency range to improve its approximation accuracy. Thus,
we set $\beta_{1}$ and $\beta_{2}$ such that the complex parts of the
principal roots of the blocks $\Lambda_{1}$ and $\Lambda_{2}$, respectively,
vanish in high-frequency regions; this setting changes the largest eigenvalue
from one to a user-defined value, which results in an approach similar to the
implicit generalized-$\alpha$ methods. Additionally, we define the limit in
which this bifurcation happens at block $i,\,i=1,\,2$, by $\Omega_{bi}$. We
define the critical stability limit for our explicit method as $\Omega_{si}$
and show the stability region of block $i$.

We find the parameter values using the characteristic equation corresponding
to the diagonal block $i$ of the amplification matrix as:

$\sum^{3}_{j=0}\left(\tilde{a}_{i}^{j}+\tilde{y}_{i}^{j}\Theta^{2}\right)\lambda_{i}^{(3-j)}=0,$
(23)

where $\tilde{a}_{i}$ and $\tilde{y}_{i}$ are functions of the parameters
$\gamma_{i}$ and $\beta_{i}$. Each block $i$ has three eigenvalues; we set two
of them to $\rho_{bi}$ and one becomes $\rho_{si}$. Therefore, we rewrite the
characteristic polynomial (23) as

$\left(\lambda_{i}^{2}+\rho_{bi}^{2}+2\lambda_{i}\rho_{bi}\right)\left(\lambda_{i}+\rho_{si}\right)=0.$
(24)

Next, to have all three roots with real values at the bifurcation limit, we
equate (24) to (22) to obtain:

$\displaystyle\Omega_{b1}$
$\displaystyle=2-2\rho_{b1}-\rho_{s1}+\rho_{s1}\rho_{b1}^{2},$ (25)
$\displaystyle\Omega_{b2}$
$\displaystyle=2+2\rho_{b2}+\rho_{s2}-\rho_{s2}\rho_{b2}^{2},$
$\displaystyle\alpha_{1}$
$\displaystyle=\frac{2-(1+\rho_{b1})\rho_{s1}}{(-1+\rho_{b1})(-1+\rho_{s1})}$
$\displaystyle\alpha_{2}$
$\displaystyle=\frac{2+(1-\rho_{b2})\rho_{s2}}{(1+\rho_{b2})(1+\rho_{s2})}.$

######  Remark 3.

Herein, we constrain $\rho_{si}\leq\rho_{bi}$ which results in
$\Omega_{bi}\leq\Omega_{si}$. We can maximize the bifurcation region
$\Omega_{bi}$ by setting $\rho_{si}=\rho_{bi}$ and consequently, obtain a one-
parameter family of algorithms, Figure 1 shows this numerically.

Figure 1 shows how the bifurcation and stability regions change as the user-
defined parameter changes.

(a) $\rho_{s1}=\rho_{s2}=\rho_{b2}=\,C$.

(b) $\rho_{s2}=\rho_{s1}=\rho_{b1}=\,C$.

Figure 1: Effect of user-defined parameters on the bifurcation regions
$\Omega_{b1}$ and $\Omega_{b2}$.

To set the complex part of the eigenvalues equal to zero in the high-frequency
regions, we define the parameters $\beta_{i}$ as

$\displaystyle\beta_{1}$
$\displaystyle=\frac{(1+\rho_{b1})(-1+\rho_{b1}\rho_{s1})^{2}}{(-1+\rho_{b1})^{2}(-1+\rho_{s1})(-2+\rho_{s1}+\rho_{b1}\rho_{s1})},$
(26) $\displaystyle\beta_{2}$
$\displaystyle=\frac{-5-3\rho_{b2}-4\rho_{s2}+2\rho_{b2}\rho_{s2}+2\rho_{b2}^{2}\rho_{s2}-\rho_{s2}^{2}+\rho_{b2}\rho_{s2}^{2}}{(1+\rho_{b2})^{2}(-2-3\rho_{s2}+\rho_{b2}\rho_{s2}-\rho_{s2}^{2}+\rho_{b2}\rho_{s2}^{2})}.$

Finally, we find the critical values at which each block’s spectral radius
becomes larger than one. For this, we set $\Theta=\Omega_{si}$ and
$\lambda_{i}=1$ in (23) and solve the resulting equation. Thus, we introduce
the critical values $\Omega_{si}$ as

$\displaystyle\Omega_{s1}$
$\displaystyle=\frac{4(1-\rho_{b1})(2-\rho_{b1}\rho_{s1}-\rho_{s1})(3+\rho_{b1}-\rho_{s1}-3\rho_{b1}\rho_{s1})}{2(5-\rho_{b1}^{2})+(5-13\rho_{b1}-\rho_{b1}^{2}-\rho_{b1}^{3})\rho_{s1}-(1+\rho_{b1})^{3}\rho_{s1}^{2}},$
(27) $\displaystyle\Omega_{s2}$
$\displaystyle=\frac{4(1+\rho_{b2})(2-\rho_{b2}\rho_{s2}+\rho_{s2})(3-\rho_{b2}+\rho_{s2}-3\rho_{b2}\rho_{s2})}{2(5-\rho_{b2}^{2})+(5-13\rho_{b2}-\rho_{b2}^{2}+\rho_{b2}^{3})\rho_{s2}-(1-\rho_{b2})^{3}\rho_{s2}^{2}}.$

######  Remark 4.

Following closely [14, 2], $\alpha_{f1}$ and $\alpha_{f2}$ are free parameters
with respect to spectral radius. Therefore, we set $\alpha_{f1}=1$ to deliver
fourth-order accuracy. Accordingly, to optimally combine low- and high-
frequency dissipation, we set $\alpha_{f2}=0$.

Therefore, by setting $0\leq\rho_{b1},\rho_{b2}\leq 1$, we control the
system’s spectral radius $\rho$ and, consequently, the high-frequency
numerical damping. Figure 2 shows how the user controls $\rho$ and the
stability region $\Omega_{s}$; setting $\rho_{b}=\rho_{s}=0.99$ leads to the
largest stability region $\Omega_{s}=4(\theta\tau)^{2}=4\Theta$, equivalent to
the stability region of the second-order central difference method [12].

Figure 2: Spectral radius $\rho$ behaviour for different $\rho_{b}$ values,
where $\rho_{b1}=\rho_{b2}=\rho_{s1}=\rho_{s2}$.

##  4 General $2k^{th}$-order accuracy in time

This section extends our approach to deliver $2k^{th}$-order accurate methods
with $k\geq 2$. That is, we introduce our $2k^{th}$-order explicit
generalized-$\alpha$ method as:

$\displaystyle MA_{n}^{\alpha_{1}}+CV_{n}+KU_{n}$
$\displaystyle=F_{n+\alpha_{f1}},$ (28) $\displaystyle
M\mathcal{L}^{3j-3}(A_{n})^{\alpha_{j}}+C\mathcal{L}^{3j-4}(A_{n})+K\mathcal{L}^{3j-5}(A_{n})$
$\displaystyle=F^{(3j-3)}_{n+\alpha_{fj}},,\qquad j=2,\cdots,k-1,$
$\displaystyle
M\mathcal{L}^{3k-3}(A_{n})^{\alpha_{k}}+C\mathcal{L}^{3k-4}(A_{n})+K\mathcal{L}^{3k-5}(A_{n})$
$\displaystyle=F^{(3k-3)}_{n+\alpha_{fk}},$

and updating the system using the following

$\displaystyle U_{n+1}$ $\displaystyle=U_{n}+\tau
V_{n}+\frac{\tau^{2}}{2}A_{n}+\frac{\tau^{3}}{6}\mathcal{L}^{1}(A_{n})+\cdots+\frac{\tau^{3k-1}}{(3k-1)!}\mathcal{L}^{3k-3}(A_{n})+\beta_{1}\tau^{2}P_{n,1},$
(29) $\displaystyle V_{n+1}$ $\displaystyle=V_{n}+\tau
A_{n}+\frac{\tau^{2}}{2}\mathcal{L}^{1}(A_{n})+\cdots+\frac{\tau^{3k-2}}{(3k-2)!}\mathcal{L}^{3k-3}(A_{n})+\tau\gamma_{1}P_{n,1},$
$\displaystyle\mathcal{L}^{3j-5}(A_{n+1})$
$\displaystyle=\mathcal{L}^{3j-5}(A_{n})+\tau\mathcal{L}^{3j-4}(A_{n})+\cdots+\frac{\tau^{3(k-j)+2}}{(3(k-j)+2)!}\mathcal{L}^{3k-3}(A_{n})+\tau^{2}\beta_{j}P_{n,j},$
$\displaystyle\mathcal{L}^{3j-4}(A_{n+1})$
$\displaystyle=\mathcal{L}^{3j-4}(A_{n})+\cdots+\frac{\tau^{3(k-j)+1}}{(3(k-j)+1)!}\mathcal{L}^{3k-3}(A_{n})+\tau\gamma_{j}P_{n,j},\qquad
j=2,\cdots,k-1,$ $\displaystyle\mathcal{L}^{3k-5}(A_{n+1})$
$\displaystyle=\mathcal{L}^{3k-5}(A_{n})+\tau{L}^{3k-4}(A_{n})+\frac{\tau^{2}}{2}\mathcal{L}^{3k-3}(A_{n})+\tau^{2}\beta_{k}\llbracket\mathcal{L}^{3k-3}(A_{n}\rrbracket,$
$\displaystyle\mathcal{L}^{3k-4}(A_{n+1})$
$\displaystyle=\mathcal{L}^{3k-4}(A_{n})+\tau\mathcal{L}^{3k-3}(A_{n})+\tau\gamma_{k}\llbracket\mathcal{L}^{3k-3}(A_{n}\rrbracket,$

with

$\displaystyle P_{n,1}$
$\displaystyle=A_{n+1}-A_{n}-\tau\mathcal{L}^{1}(A_{n})-\cdots-\frac{\tau^{3k-3}}{(3k-3)!}\mathcal{L}^{3k-3}(A_{n}),$
(30) $\displaystyle A_{n}^{\alpha_{1}}$
$\displaystyle=A_{n}+\tau\mathcal{L}^{1}(A_{n})+\cdots+\frac{\tau^{3k-3}}{(3k-3)!}\mathcal{L}^{3k-3}(A_{n})+\alpha_{1}P_{n,1},$
$\displaystyle P_{n,j}$
$\displaystyle=\mathcal{L}^{3j-3}(A_{n+1})-\mathcal{L}^{3j-3}(A_{n})-\cdots-\frac{\tau^{3(k-j)}}{3(k-j)!}\mathcal{L}^{3k-3}(A_{n}),$
$\displaystyle A_{n}^{\alpha_{j}}$
$\displaystyle=\mathcal{L}^{3j-3}(A_{n})+\cdots+\frac{\tau^{3(k-j)}}{3(k-j)!}\mathcal{L}^{3k-3}(A_{n})+\alpha_{j}P_{n,j},\qquad
j=2,\cdots,k-1,$ $\displaystyle\mathcal{L}^{3k-3}(A_{n})^{\alpha_{k}}$
$\displaystyle=\mathcal{L}^{3k-3}(A_{n})+\alpha_{k}\llbracket\mathcal{L}^{3k-3}(A_{n})\rrbracket.$

Letting $k=2$ recovers the fourth-order explicit generalized-$\alpha$ method
of Section 3. Using similar arguments to the proof of Theorem 1, we can
establish higher-order schemes in the form of (28) and (29) by using high-
order Taylor series. Then, to seek $2k^{th}$-order of accuracy, we substitute
(29) into (28) and find a matrix system as

$L{\boldsymbol{U}}_{n+1}=R{\boldsymbol{U}}_{n}+\mathbb{F}_{n+\alpha_{f}}.$
(31)

Therefore, the amplification matrix of the $2k^{th}$-order accurate scheme
becomes:

$G_{k}=L^{-1}R=\begin{bmatrix}\mathbb{\Lambda}_{1}&\mathbb{\Xi}_{12}&\cdots&\cdots&\mathbb{\Xi}_{1k}\\\
\boldsymbol{0}&\mathbb{\Lambda}_{2}&\mathbb{\Xi}_{23}&\cdots&\mathbb{\Xi}_{2k}\\\
\vdots&&\ddots&\\\
\boldsymbol{0}&\boldsymbol{0}&\cdots&\mathbb{\Lambda}_{k-1}&\mathbb{\Xi}_{k-1}\\\
\boldsymbol{0}&\boldsymbol{0}&\cdots&\boldsymbol{0}&\mathbb{\Lambda}_{k}\end{bmatrix},$
(32)

with

$\displaystyle\mathbb{\Lambda}_{j}$
$\displaystyle=\dfrac{1}{\alpha_{j}}\begin{bmatrix}\alpha_{j}-{\beta_{j}\Theta}&\alpha_{j}-{\beta_{j}\Theta}&\frac{1}{2}\left({\alpha_{j}-\beta_{j}\left(\Theta+2\right)}\right)\\\
-{\gamma_{j}\Theta}&\alpha_{j}-{\gamma_{j}\Theta}&\alpha_{j}-\frac{\gamma_{j}}{2}\left(\Theta+2\right)\\\
-{\Theta}&-{\Theta}&\alpha_{j}-1-\frac{\Theta}{2}\\\ \end{bmatrix},\qquad
j=1,\cdots,k-1$ (33) $\displaystyle\mathbb{\Lambda}_{k}$
$\displaystyle=\dfrac{1}{\alpha_{k}}\begin{bmatrix}\alpha_{k}-{\beta_{k}\Theta}&\alpha_{k}&\frac{\alpha_{k}}{2}-{\beta_{k}}\\\
-{\gamma_{k}\Theta}&\alpha_{k}&\alpha_{k}-{\gamma_{k}}\\\
-{\Theta}&0&{\alpha_{k}-1}\end{bmatrix}.$

######  Remark 5.

The amplification matrix (32) is an upper-triangular block matrix; thus, we
neglect the non-diagonal block contributions in the eigenvalue analysis.

######  Theorem 2.

Assuming $u_{h}(t)$ and $\dot{u}_{h}(t)$ have sufficient regularity in time,
our semi-discrete method (28)-(30) for advancing (2) is $2k^{th}$-order
accurate in time when

$\displaystyle\gamma_{j}$
$\displaystyle=\frac{1}{2}-\alpha_{fj}+\alpha_{j},\qquad
j=1,\cdots,k-1,\quad\gamma_{k}$
$\displaystyle=\frac{1}{2}-\alpha_{fk}+\alpha_{k}.$ (34)

###### Proof.

The amplification matrix (32) is an upper-diagonal block matrix; each block is
a $3\times 3$ matrix. The last two diagonal blocks $\Lambda_{k-1}$ and
$\Lambda_{k}$ have similar entries to those of the amplification matrix of the
fourth-order method. The other diagonal blocks are analogous to
$\Lambda_{k-1}$. Hence, separately for each block, we determine the relevant
terms for (17) and consider the higher-order terms to obtain second-order
accuracy. Consequently, after solving the whole system and adding the higher-
order terms to the unknowns $u_{h}$ and $v_{h}$, we have a truncation error of
$\mathcal{O}(\tau^{2k+1})$. ∎

###  4.1 CFL condition and dissipation control

We define the system’s spectral behaviour, the bifurcation regions, and its
CFL conditions; thus, we study the amplification matrix’s eigenvalues (32).
Accordingly, we calculate each diagonal block’s eigenvalues. The first $k-1$
blocks have identical structures; hence, following Section 3, we propose these
algorithmic parameters as

$\boxed{\begin{aligned}
\Omega_{bj}&=2-2\rho_{bj}-\rho_{sj}+\rho_{sj}\rho_{bj}^{2},\qquad
j=1,\cdots,k-1,\\\
\alpha_{j}&=\frac{2-(1+\rho_{bj})\rho_{sj}}{(-1+\rho_{bj})(-1+\rho_{sj})},\\\
\beta_{j}&=\frac{(1+\rho_{bj})(-1+\rho_{bj}\rho_{sj})^{2}}{(-1+\rho_{bj})^{2}(-1+\rho_{sj})(-2+\rho_{sj}+\rho_{bj}\rho_{sj})},\\\
\Omega_{bk}&=2+2\rho_{bk}+\rho_{sk}-\rho_{sk}\rho_{bk}^{2},\\\
\alpha_{k}&=\frac{2+(1-\rho_{bk})\rho_{sk}}{(1+\rho_{bk})(1+\rho_{sk})},\\\
\beta_{k}&=\frac{-5-3\rho_{bk}-4\rho_{sk}+2\rho_{bk}\rho_{sk}+2\rho_{bk}^{2}\rho_{sk}-\rho_{sk}^{2}+\rho_{bk}\rho_{sk}^{2}}{(1+\rho_{bk})^{2}(-2-3\rho_{sk}+\rho_{bk}\rho_{sk}-\rho_{sk}^{2}+\rho_{bk}\rho_{sk}^{2})},\end{aligned}}$
(35)

and the critical values $\Omega_{sj}$ are

$\displaystyle\Omega_{sj}$
$\displaystyle=\frac{4(1-\rho_{bj})(2-\rho_{bj}\rho_{sj}-\rho_{sj})(3+\rho_{bj}-\rho_{sj}-3\rho_{bj}\rho_{sj})}{2(5-\rho_{bj}^{2})+(5-13\rho_{bj}-\rho_{bj}^{2}-\rho_{bj}^{3})\rho_{sj}-(1+\rho_{bj})^{3}\rho_{sj}^{2}},\quad
j=1,\cdots,k-1,$ (36) $\displaystyle\Omega_{sk}$
$\displaystyle=\frac{4(1+\rho_{bk})(2-\rho_{bk}\rho_{sk}+\rho_{sk})(3-\rho_{bk}+\rho_{sk}-3\rho_{bk}\rho_{sk})}{2(5-\rho_{bk}^{2})+(5-13\rho_{bk}-\rho_{bk}^{2}+\rho_{bk}^{3})\rho_{sk}-(1-\rho_{bk})^{3}\rho_{sk}^{2}}.$

######  Remark 6.

We constrain $\rho_{sm}\leq\rho_{bm},\,m=1,\cdots,k$ and maximize the
bifurcation regions $\Omega_{bm}$ by setting $\rho_{sm}=\rho_{bm}=\rho$ with
$0\leq\rho<1$ as a user-defined parameter; our method is a one-parameter
family of time-marching algorithms. Additionally, the stability region is
independent of the accuracy order; we obtain higher-order accuracy without
affecting any features of the second-order algorithm (i.e., preserve stability
regions, bifurcation limit, and dissipation control). Similarly, $\alpha_{fj}$
and $\alpha_{fk}$ are free parameters; thus, we set
$\alpha_{fj}=1,\,\,j=1,\,\cdots,\,k-1$, and $\alpha_{fk}=0$.

######  Remark 7.

The most expensive computational cost of our explicit time marching is the
factorization of the mass matrices. Next, we discuss a state-of-the-art
approach to precondition the mass matrices resulting from isogeometric
analysis that minimize this cost for complex geometries as well as for single-
and multi-patch discretizations.

##  5 Solver

The method we propose has many advantages: for instance, the CFL condition and
the dissipation control are independent of the order of accuracy; however, it
requires solving $k$ mass-matrix systems. Herein, we discuss a method to
accelerate these systems’ solution when using isogeometric analysis. Firstly,
we describe efficient and robust preconditioners for maximum-continuity
isogeometric mass matrices for single and multi-patch geometries. Then, we
adopt an iterative solver (i.e., preconditioned conjugate gradients, PCG) to
calculate the solutions.

###  5.1 B-splines

Given two positive integers $p$ and $m$, consider an open knot vector

$\Xi:=\\{\xi_{1},\ldots,\xi_{m+p+1}\\}$

such that

$\xi_{1}=\ldots=\xi_{p+1}<\xi_{p+2}\leq\ldots\leq\xi_{m}<\xi_{m+1}=\ldots=\xi_{m+p+1},$

where interior repeated knots are allowed with maximum multiplicity $p$. We
assume $\xi_{1}=0$ and $\xi_{m+p+1}=1$. From the knot vector $\Xi$, we define
degree-$p$ B-spline functions using the Cox-De Boor recursive formula: we
start with piecewise constants ($p=0$):

$\widehat{b}_{i,0}(\zeta)=\left\\{\begin{array}[]{ll}1&\text{if
}\xi_{i}\leq\zeta<\xi_{i+1},\\\ 0&\text{otherwise},\end{array}\right.$

for $p\geq 1$, the following recursion defines the B-spline functions

$\widehat{b}_{i,p}(\zeta)=\frac{\zeta-\xi_{i}}{\xi_{i+p}-\xi_{i}}\widehat{b}_{i,p-1}(\zeta)+\frac{\xi_{i+p+1}-\zeta}{\xi_{i+p+1}-\xi_{i+1}}\widehat{b}_{i+1,p-1}(\zeta),$

where $0/0=0$. Each B-spline $\widehat{b}_{i,p}$ depends only on $p+2$ knots,
which we collect in a local knot vector

$\Xi_{i,p}:=\\{\xi_{i},\ldots,\,\xi_{i+p+1}\\},$

is non-negative and its support is the interval $[\xi_{i},\xi_{i+p+1}]$.
Moreover, these B-spline functions constitute a partition of unity, that is

$\displaystyle\sum_{i=1}^{m}\widehat{b}_{i,p}(x)=1,$ $\displaystyle\forall
x\in(0,1).$ (37)

The univariate spline space is

$\widehat{\mathcal{S}}_{h}=\widehat{\mathcal{S}}_{h}([0,1]):=\mathrm{span}\\{\widehat{b}_{i,p}\\}_{i=1}^{m},$

where $h$ denotes the maximal mesh-size. We may drop the degree $p$ from the
notation when it will not lead to confusion (see, e.g., [15, 16]).

We define multivariate B-splines from univariate ones by tensorization, as is
common practice. Let $d$ be the space dimension and consider open knot vectors
${\Xi_{k}=\\{\xi_{k,1},\ldots,\xi_{k,m+p+1}\\}}$ and a set of multi-indices
${\mathbf{I}:=\\{\mathbf{i}=(i_{1},\ldots,i_{d}):\,1\leq i_{l}\leq m\\}}$. For
each multi-index $\mathbf{i}=(i_{1},\ldots,i_{d})$, we introduce the
$d$-variate B-spline,

$\widehat{B}_{\mathbf{i}}(\mathbf{\zeta}):=\widehat{b}[\Xi_{i_{1},p}](\zeta_{1})\ldots\widehat{b}[\Xi_{i_{d},p}](\zeta_{d}).$

The corresponding spline space is

$\widehat{\mathcal{S}}_{h}=\widehat{\mathcal{S}}_{h}([0,1]^{d}):=\mathrm{span}\left\\{B_{\mathbf{i}}:\,\mathbf{i}\in\mathbf{I}\right\\},$

where $h$ is the maximal mesh-size in all dimensions, that is,

$h:=\max_{\begin{subarray}{c}1\leq k\leq d\\\ 1\leq i\leq
m+p+1\end{subarray}}\\{|\xi_{k,i+1}-\xi_{k,i}|\\}.$

######  Assumption 1.

Knot vectors are quasi-uniform; there exists $\alpha>0$, $h$ independent, such
that each non-empty knot span $(\xi_{k,i},\xi_{k,i+1})$ fulfills $\alpha
h\leq\xi_{k,i+1}-\xi_{k,i}$, for $1\leq k\leq d$.

###  5.2 Single-patch geometric space

We consider a single-patch domain $\Omega\subset\mathbb{R}^{d}$, a
$d$-dimensional parametrization $\boldsymbol{F}$,

$\displaystyle\Omega=\boldsymbol{F}(\widehat{\Omega}),$
$\displaystyle\text{with
}\boldsymbol{F}(\boldsymbol{\xi})=\sum_{\mathbf{i}}\boldsymbol{C}_{\mathbf{i}}\widehat{B}_{\mathbf{i}}(\boldsymbol{\xi}),$

where $\boldsymbol{C}_{\mathbf{i}}$ are the control points and
$\widehat{B}_{\mathbf{i}}$ are tensor-product B-spline basis functions defined
on a parametric patch $\widehat{\Omega}:=(0,1)^{d}$, where $\boldsymbol{F}$ is
an invertible map. Following the isoparametric paradigm, isogeometric basis
functions $B_{\mathbf{i}}$ are the push-forward of the parametric basis, that
is, $B_{\mathbf{i}}=\widehat{B}_{\mathbf{i}}\circ\boldsymbol{F}^{-1}$. Thus,
the isogeometric space on $\Omega$ is defined as

$\mathcal{S}_{h}=\mathcal{S}_{h}(\Omega):=\mathrm{span}\left\\{B_{\mathbf{i}}:=\widehat{B}_{\mathbf{i}}\circ\mathbf{F}^{-1}\
:\ \mathbf{i}\in\mathbf{I}\right\\}.$

We introduce a co-lexicographical reordering of the basis functions and write

$\mathcal{S}_{h}=\mathrm{span}\left\\{B_{\mathbf{i}}:\,\mathbf{i}\in\mathbf{I}\right\\}=\mathrm{span}\left\\{B_{i}\right\\}_{i=1}^{N_{\text{dof}}}.$
(38)

###  5.3 Multi-patch B-splines

Following [17], a multi-patch domain $\Omega\subset\mathbb{R}^{d}$ is an open
set, a subdomain union

$\overline{\Omega}=\bigcup_{{r}=1}^{N_{\mathrm{patch}}}\overline{\Omega^{({r})}},$
(39)

where $N_{\mathrm{patch}}$ is the number of subdomains,
$\Omega^{({r})}=\boldsymbol{F}^{({r})}(\widehat{\Omega})$ are the disjoint
patches (subdomain pre-images), each $\boldsymbol{F}^{({r})}$ has a different
spline parametrization, and the super index $({r})$ refers to
$\Omega^{({r})}$. We introduce for each patch $\Omega^{({r})}$, B-spline
spaces

$\displaystyle\widehat{\mathcal{S}}^{({r})}_{h}:=\mathrm{span}\left\\{\widehat{B}^{({r})}_{i}:\,i=1,\ldots,N_{\text{dof}}^{({r})}\right\\}.$

and isogeometric spaces

$\displaystyle\mathcal{S}^{({r})}_{h}:=\mathrm{span}\left\\{B^{({r})}_{i}\ :\
i=1,\ldots,N_{\text{dof}}^{({r})}\right\\}.$

We assume for simplicity that all patches have the same degree $p$. We define
an isogeometric space on $\Omega$ by imposing continuity at the interfaces
between patches, that is

$V_{h}:=\left\\{v\in
C^{0}(\Omega):v|_{\Omega^{({r})}}\in\mathcal{S}^{({r})}_{h}\text{ for
}{r}=1,\ldots,N_{\mathrm{patch}}\right\\}.$ (40)

We assume a suitable conformity to construct a basis for space $V_{h}$; for
all $r,s\in\\{1,\ldots,N_{\mathrm{patch}}\\}$, with $r\neq s$, let
$\Gamma_{{r}{s}}=\partial\Omega^{({r})}\cap\partial\Omega^{({s})}$ be the
interface between the patches $\Omega^{({r})}$ and $\Omega^{({s})}$.

######  Assumption 2.

We assume:

  1. 1.

$\Gamma_{{r}{s}}$ is either a vertex or the image of a full edge or the image
of a full face for both parametric domains.

  2. 2.

For each $B^{({r})}_{\mathbf{i}}\in\mathcal{S}_{h}^{({r})}$ such that
$\mathrm{supp}(B^{({r})}_{\mathbf{i}})\cap\Gamma_{{r}{s}}\neq\emptyset$, there
exists a function $B^{({s})}_{\mathbf{j}}\in\mathcal{S}_{h}^{({s})}$ such that
$B^{({r})}_{\mathbf{i}}|_{\Gamma_{{r}{s}}}=B^{({s})}_{\mathbf{j}}|_{\Gamma_{{r}{s}}}$.

We define, for each patch $\Omega^{({r})}$, an application

$G^{({r})}:\\{1,\ldots,N_{\text{dof}}^{(r)}\\}\rightarrow\mathcal{J}=\\{1,\ldots,\dim(V_{h})\\},$

such that $G^{({r})}(i)=G^{({s})}(j)$ if and only if
$\Gamma_{{r}{s}}\neq\emptyset$ and
$B^{({r})}_{i}|_{\Gamma_{{r}{s}}}=B^{({s})}_{j}|_{\Gamma_{{r}{s}}}$. Moreover,
we define, for each global index $l\in\mathcal{J}$, a set of pairs
$\mathcal{J}_{l}:=\\{({r},i):\,G^{({r})}(i)=l\\}$, which collects local
indices of patch-wise contributions to a global function, and the scalar

$n_{l}:=\\#\mathcal{J}_{l},$ (41)

that expresses the patch multiplicity for the global index $l$. Furthermore,
let

$N_{\mathrm{adj}}:=\max\\{n_{l}:\,l\in\mathcal{J}\\}$ (42)

be the maximum number of adjacent patches (i.e., those with non-empty closure
intersection). We define, for each $l\in\mathcal{J}$, the global basis
function

$B_{l}(\boldsymbol{x}):=\begin{cases}B^{({r})}_{i}(\boldsymbol{x})&\text{ if
}\boldsymbol{x}\in\overline{\Omega^{({r})}}\text{ and
}({r},i)\in\mathcal{J}_{l},\\\ 0&\text{ otherwise},\end{cases}$ (43)

which is continuous due to Assumption 2. Then

$V_{h}=\mathrm{span}\\{B_{l}:\,l\in\mathcal{J}\\}.$ (44)

The set $\\{B_{l}:\,l\in\mathcal{J}\\}$ where $B_{l}$ is defined as in (43),
represents a basis for $V_{h}$. Finally, we define the index set
$\mathcal{J}^{(r)}\subset\mathcal{J}$ such that $l\in\mathcal{J}^{(r)}$ if and
only if $l=G^{(r)}(i)$ for some $i$, where
$\\#\mathcal{J}^{(r)}=N_{\text{dof}}^{(r)}$ and $\mathcal{J}^{(r)}$ are the
index set for $\widehat{\mathcal{S}}_{h}^{(r)}$ and $\mathcal{S}_{h}^{(r)}$,
abusing notation.

###  5.4 Mass preconditioner on a patch

In this section, we briefly revisit the preconditioner described in [7], for
an isogeometric mass matrix associated with a single-patch domain, denoted
$\Omega$, that is

$[\mathbf{M}]_{i,j}=\int_{\widehat{\Omega}}\widehat{B}_{i}\widehat{B}_{j}|\det(D\boldsymbol{F})|.$
(45)

######  Assumption 3.

We assume $\boldsymbol{F}\in C^{1}([0,1]^{d})$ and there exists $\delta>0$
such that for all $\boldsymbol{x}\in[0,1]^{d}$,
$\det(D\boldsymbol{F}^{(r)})\geq\delta$.

Hence, as a preconditioner for the mass matrix $\mathbf{M}$, we consider

$\boldsymbol{\mathcal{M}}:=\mathbf{D}^{\frac{1}{2}}\widehat{\mathbf{D}}^{-\frac{1}{2}}\widehat{\mathbf{M}}\widehat{\mathbf{D}}^{-\frac{1}{2}}\mathbf{D}^{\frac{1}{2}},$
(46)

where

$\displaystyle[\widehat{\mathbf{M}}]_{i,j}:=\int_{\widehat{\Omega}}\widehat{B}_{i}\widehat{B}_{j},$
$\displaystyle\widehat{\mathbf{D}}:=\text{diag}\left(\widehat{\mathbf{M}}\right),$
$\displaystyle\mathbf{D}:=\text{diag}\left(\mathbf{M}\right).$ (47)

We define

$\displaystyle\kappa(\mathbf{A}):=\frac{\lambda_{\mathrm{max}}(\mathbf{A})}{\lambda_{\mathrm{min}}(\mathbf{A})},$
(48)

for a symmetric, positive definite matrix $\mathbf{A}$. This preconditioner
(46) has several good properties (see [7] for details):

  * •

asymptotic exactness, that is

$\lim_{h\rightarrow
0}\kappa(\boldsymbol{\mathcal{M}}^{-\frac{1}{2}}\mathbf{M}\boldsymbol{\mathcal{M}}^{-\frac{1}{2}})=1;$
(49)

  * •

$p$-robustness, in fact numerical tests show a slow growth (almost linear) of
$\kappa(\boldsymbol{\mathcal{M}}^{-\frac{1}{2}}\mathbf{M}\boldsymbol{\mathcal{M}}^{-\frac{1}{2}})$
with respect to the spline degree $p$.

  * •

good behaviour with respect to the spline parametrization $\boldsymbol{F}$;
numerical evidence shows that
$\kappa(\boldsymbol{\mathcal{M}}^{-\frac{1}{2}}\mathbf{M}\boldsymbol{\mathcal{M}}^{-\frac{1}{2}})$
and the number of PCG iteration needed to converge are small even for
distorted geometries.

###  5.5 Mass preconditioner on multi-patch domain

We introduce a mass matrix preconditioner for multi-patch domains, that is

$\displaystyle[M]_{i,j}=\int_{\Omega}B_{i}B_{j},$ $\displaystyle
B_{i},B_{j}\in\mathcal{B}_{\text{mp}},$

where $\Omega$ is the union of $\Omega^{(r)}$, see (39). Following [7], we
combine the single-patch preconditioner of (46) with an additive Schwarz
method. We define a family of local spaces

$\displaystyle
V_{h}^{(r)}:=\mathrm{span}\left\\{\mathcal{B}_{\text{mp}}^{(r)}\right\\},$
$\displaystyle r=1,\ldots,N_{\mathrm{patch}},$ (50)

where we have set

$\mathcal{B}_{\text{mp}}^{(r)}:=\\{B_{l}:\,l\in\mathcal{J}^{(r)}\\},$

with $\mathcal{J}^{(r)}$ defined as in Section 5.3. Therefore, $V_{h}^{(r)}$
is the subspace of $V_{h}$ spanned by the B-splines whose support intersect
$\Omega^{(r)}$. Moreover, following the notation of [18], we consider
restriction operators $R^{(r)}:V_{h}\rightarrow V_{h}^{(r)}$ with
$r=1,\ldots,N_{\mathrm{patch}}$, defined by

$\displaystyle
R^{(r)}\left(\sum_{l\in\mathcal{J}}u_{l}B_{l}\right)=\sum_{l\in\mathcal{J}^{(r)}}u_{l}B_{l}$

and their transposes, in the basis representation,
${{R^{(r)}}^{T}:V_{h}^{(r)}\rightarrow V_{h}}$ correspond, in our case, to the
inclusion of $V_{h}^{(r)}$ into $V_{h}$. We denote with $\mathbf{R}^{(r)}$ and
$\mathbf{R}^{(r)^{T}}$ the rectangular matrices associated to ${R^{(r)}}$ and
${R^{(r)}}^{T}$, respectively.

The additive Schwarz preconditioner (inverse) is

$\boldsymbol{\mathcal{M}}^{-1}_{\text{ad}}:=\sum_{{r}=1}^{N_{\mathrm{patch}}}{\mathbf{R}^{(r)}}^{T}{\boldsymbol{\mathcal{M}}^{(r)}}^{-1}\mathbf{R}^{(r)},$
(51)

where we set

$\boldsymbol{\mathcal{M}}^{(r)}:=\mathbf{D}^{(r)^{\frac{1}{2}}}\widehat{\mathbf{D}}^{(r)^{-\frac{1}{2}}}\widehat{\mathbf{M}}^{(r)}\widehat{\mathbf{D}}^{(r)^{-\frac{1}{2}}}\mathbf{D}^{(r)^{\frac{1}{2}}}$
(52)

and

$\displaystyle[\widehat{\mathbf{M}}^{(r)}]_{i,j}:=\int_{\widehat{\Omega}}\widehat{B}^{(r)}_{i}\widehat{B}^{(r)}_{j},$
$\displaystyle\widehat{\mathbf{D}}^{(r)}:=\text{diag}\left(\widehat{\mathbf{M}}^{(r)}\right),$
$\displaystyle\mathbf{D}^{(r)}:=\text{diag}\left(\mathbf{M}^{(r)}\right),$
(53)

assumption the basis functions ordering of
$\\{\widehat{B}^{(r)}_{i}\\}_{i=1}^{{N_{\text{dof}}^{(r)}}}$ and
$\\{{B}^{(r)}_{i}\\}_{i=1}^{N_{\text{dof}}^{(r)}}$ follows Section 5.3.

######  Assumption 4.

For all $r=1,\ldots,N_{\mathrm{patch}}$, let us assume $\boldsymbol{F}^{(r)}$
fulfils Assumption 3.

An upper bound for the conditioning of the preconditioned system [7] is

$\displaystyle
k\left(\boldsymbol{\mathcal{M}}^{-\frac{1}{2}}_{\text{ad}}\mathbf{M}\boldsymbol{\mathcal{M}}^{-\frac{1}{2}}_{\text{ad}}\right)\leq
CN_{\mathrm{adj}}^{2},$

where the constant $C$ is independent of $h$ and $N_{\mathrm{patch}}$;
equation (42) defines $N_{\mathrm{adj}}$. Numerical experiments show the
method’s $p$-robustness and good behavior with respect to the spline
parametrizations $\boldsymbol{F}^{(r)}$.

###  5.6 Preconditioners application cost

The mass matrices and preconditioners of Section 5.4 and 5.5 are symmetric and
positive definite; thus, we adopt the preconditioned conjugate gradient method
(PCG) to solve the associated linear systems. The application of the single-
patch preconditioner requires the solution of a linear system associated with

$\displaystyle\boldsymbol{\mathcal{M}}=\mathbf{D}^{\frac{1}{2}}\widehat{\mathbf{D}}^{-\frac{1}{2}}\widehat{\mathbf{M}}\widehat{\mathbf{D}}^{-\frac{1}{2}}\mathbf{D}^{\frac{1}{2}}.$

the Kronecker product structure of the matrix implies that the
preconditioner’s application cost is proportional to $N_{\text{dof}}$ [19, 20,
21, 7].

For the multi-patch cases, the application of
$\boldsymbol{\mathcal{M}}^{-1}_{\text{ad}}$, (51), involves, for
$r\in\\{1,\ldots,N_{\mathrm{patch}}\\}$, the application of the operators
$R^{(r)}$ and ${R^{(r)}}^{T}$, whose cost is negligible, and the application
of $\left({\boldsymbol{\mathcal{M}}^{(r)}}\right)^{-1}$, whose cost we analyze
above. In conclusion, also for $\boldsymbol{\mathcal{M}}^{-1}_{\text{ad}}$,
the cost of application is $O(N_{\text{dof}})$ floating-point operations
(FLOPS).

##  6 Numerical results

In this section, firstly, we show the performance of the proposed explicit
generalized-$\alpha$ method. Then, we provide further numerical results on our
solver. All the tests are performed with Matlab R2015a and GeoPDEs toolbox
[22]. To show how our explicit method and solver work, the linear systems are
solved by PCG, with tolerance equal to $10^{-12}$ and with the null vector as
the initial guess. We denote by $n_{\mathrm{sub}}$ the number of subdivisions,
which are the same in each parametric direction and in each patch, by $p$ the
spline degree and by $\tau$ the size of the time grid. Moreover, we underline
that we only consider the splines of maximal regularity. The symbol “*”
denotes the impossibility of the formation of the matrix $\mathbf{M}$, due to
memory requirements.

In all the examples, we consider the following hyperbolic problem

$\displaystyle\begin{cases}\ddot{u}(\boldsymbol{x},t)-\nabla\cdot(\omega^{2}\,\nabla
u(\boldsymbol{x},t))&=f(\boldsymbol{x},t),\qquad\,\,\,(\boldsymbol{x},t)\in\Omega\times(0,T],\\\
u(\boldsymbol{x},t)&=u_{D},\qquad\qquad\boldsymbol{x}\in\partial\Omega,\\\
u(\boldsymbol{x},0)&=u_{0},\qquad\,\qquad\boldsymbol{x}\in\Omega,\\\
\dot{u}(\boldsymbol{x},0)&=v_{0},\qquad\,\,\qquad\boldsymbol{x}\in\Omega,\end{cases}$
(54)

where $\omega$ is assumed to be one except in the problem taken into account
in section 6.2.

###  6.1 Convergence of the generalized-$\alpha$ method

For verifying the accuracy of the $4^{th}$ order explicit generalized-$\alpha$
method, we solve (54) on $\Omega=[0,1]\times[0,1]$, choosing the source term,
Dirichlet boundary condition and initial condition such that the exact
solution reads

$\displaystyle u(x,y,t)=\sin(10\pi x)\sin(10\pi
x)\left[\cos\left(10\sqrt{2}\pi t\right)+\sin\left(10\sqrt{2}\pi
t\right)\right].$

In Figure 3, we show the $4^{th}$ order of convergence in the time domain for
the unknown displacement and velocity at final time $T=0.1$. In this example,
we use $100\times 100$ isogeometric elements with polynomial degree of $p=7$
with continuity $C^{6}$. We also set all the user-defined parameters
$\rho_{bi}=\rho_{si}=\rho$. As the analysis showed before, the obtained
solutions $u^{n}_{h}$ and $v^{n}_{h}$ which respectively approximate
$u(x,t_{n})$ and $\dot{u}(x,t_{n})$ converge with an order of four.

Figure 3: $L^{2}(\Omega)$ norm error at $T=0.1$ for different values of user-
defined parameters.

###  6.2 Dispersion behaviour of the generalized-$\alpha$ method

In this example, we numerically show that choosing higher-order methods leads
to better approximations in low-frequency zones. Again, we consider a
homogeneous Dirichlet boundary condition with an exact solution given by

$u(x,t)=\sin(j\pi x)\cos(\pi t),$ (55)

and setting $\omega=\frac{1}{j}$. We discretize the spatial domain using
$N=400$ elements with polynomials of degree $p=4$ and regularity $C^{3}$. We
refer to Figure 4, where we consider two cases; one is with the time step
$\tau=0.05$ in which the bifurcation region is reached and, therefore, we
obtain distinguishable results by setting different values for $\rho$.
Additionally, it shows the importance of picking $\rho$ wisely. In the other
example, we set the time step $\tau=10^{-3}$. Here, the spectral behaviour of
the system does not change by varying $j$ and results in similar solutions for
different values of $\rho$.

Figure 4: Comparison between the dispersion behaviour of the second and fourth
order methods, $\tau=0.05$ (top) and $\tau=10^{-3}$ (bottom) for final time
$t_{n}=5$.

###  6.3 Preconditioner Performance

In order to analyze the behaviour of the proposed preconditioners and the
proposed fourth-order method, we solve (54) on a regular single-patch domain
(Figure 5(a)), a singular one (Figure 6(a)) and a multi-patch domain (Figure
7(a)), obtained by glueing together seven blade-shaped patches like the one
represented in Figure 5(a). For each of them, source term, Dirichlet boundary
condition and initial condition are chosen such that the exact solution is
always

$\displaystyle
u(x_{1},x_{2},x_{3},t)=\sin(x_{1})\sin(x_{2})\sin(x_{3})\left[\cos(20\pi
t)+\sin(20\pi t)\right].$

(a) single-patch domain (blade geometry)

(b) Optimal convergence of $L_{2}$ error

Figure 5: $L^{2}(\Omega)$ norm relative error at $T=64\cdot\tau$ with
$\tau=10^{-5}$ for the blade geometry.

(a) Singular domain (donut geometry)

(b) Optimal convergence of $L_{2}$ error

Figure 6: $L^{2}(\Omega)$ norm relative error at $T=64\cdot\tau$ with
$\tau=10^{-5}$ for the donut geometry.

(a) multi-patch domain (fan geometry)

(b) Optimal convergence of $L_{2}$ error

Figure 7: $L^{2}(\Omega)$ norm relative error at $T=64\cdot\tau$ with
$\tau=10^{-5}$ for the fan geometry.

We report the mean value, across all the time steps, of the number of
iterations needed by PCG for reaching the given tolerance for the three
different spatial domains. Tables 1 and 3 shows that the number of PCG
iterations is always very low and decreases when the subdivisions are
increased. This is true also for the singularly parametrized Donut domain (see
Table 2), even though this case is beyond the robustness result presented in
Section 5.4.

$n_{\mathrm{sub}}$ | $p=1$ | $p=2$ | $p=3$ | $p=4$
---|---|---|---|---
8 | 7.5 | 10.1 | 10.2 | 10.4
16 | 6.1 | 8.4 | 8.4 | 8.9
32 | 5.5 | 6.5 | 6.5 | 7.2
64 | 5.0 | 5.5 | 5.5 | 6.0
Table 1: Mean values, across all the time steps, of the iterations needed by PCG on the blade, for $\tau=10^{-5}$ and $T=64\cdot\tau$. $n_{\mathrm{sub}}$ | $p=1$ | $p=2$ | $p=3$ | $p=4$
---|---|---|---|---
8 | 6.0 | 8.0 | 9.0 | 10.0
16 | 5.0 | 7.0 | 8.0 | 8.0
32 | 5.0 | 6.0 | 7.0 | 8.0
64 | 5.3 | 6.5 | 7.0 | 7.8
Table 2: Mean values, across all the time steps, of the iterations needed by PCG on the donut, for $\tau=10^{-5}$ and $T=64\cdot\tau$. $n_{\mathrm{sub}}$ | $p=1$ | $p=2$ | $p=3$ | $p=4$
---|---|---|---|---
8 | 22.5 | 32.0 | 37.7 | 41.9
16 | 21.0 | 27.8 | 30.8 | 35.6
32 | 19.0 | 25.5 | 27.9 | 29.9
64 | 17.0 | 22.5 | 26.0 | 27.5
Table 3: Mean values, across all the time steps, of the iterations needed by
PCG on the fan problem, for $\tau=5\cdot 10^{-5}$ and $T=64\cdot\tau$.

Furthermore, for assessing the good behaviour of the isogeometric
discretization, for each spatial domain, we report the relative error in
$L^{2}(\Omega)$ at the final instant $T=6.4\cdot 10^{-4}$ with $\tau=10^{-5}$
for different mesh size and spline degree. In Figure 5(b), we can see that the
rates of convergence are optimal with respect to the mesh size $h\approx
n_{\mathrm{sub}}^{-1}$, i.e., of order $O(h^{p+1})$, for $p=1,2,3,4$, as
expected from standard a priori error estimate. This optimal convergence in
the spatial domain is also true for singular and multi-patch domains. For
this, we refer the reader to Figures 6(b) and 7(b).

##  7 Contributions

We propose a new class of higher-order explicit generalized-$\alpha$ methods
for solving hyperbolic problems that provide dissipation control.
Additionally, the method’s stability region is independent of the accuracy in
time. We obtain $2k^{th}$-order of accuracy by solving $k$ mass-matrix
systems. We adopt a preconditioner designed and built for the isogeometric
mass matrix, which significantly reduces the computational costs. We discuss
several numerical examples that show the stability and performance of our one-
parameter family of explicit time-marching methods.

## Acknowledgement

This publication was also made possible in part by the CSIRO Professorial
Chair in Computational Geoscience at Curtin University and the Deep Earth
Imaging Enterprise Future Science Platforms of the Commonwealth Scientific
Industrial Research Organisation, CSIRO, of Australia. This project has
received funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Sklodowska-Curie grant agreement No
777778 (MATHROCKS). The Curtin Corrosion Centre and the Curtin Institute for
Computation kindly provide ongoing support. The Australian Government Research
Training Program Scholarship supported P. Behnoudfar’s research. Part of this
work was carried over while P. Behnoudfar was invited by Prof. A. Reali in
Pavia, partially supported by MIUR-PRIN project XFAST-SIMS (no. 20173C478N).
G. Loli and G. Sangalli were partially supported by the European Research
Council through the FP7 Ideas Consolidator Grant HIGEOM n.616563, and by the
Italian Ministry of Education, University and Research (MIUR) through the
”Dipartimenti di Eccellenza Program (2018-2022) - Dept. of Mathematics,
University of Pavia”. They are also members of the Gruppo Nazionale Calcolo
Scientifico - Istituto Nazionale di Alta Matematica (GNCS-INDAM). These
supports are gratefully acknowledged.

## References

  * [1] J. Chung and G. M. Hulbert.  A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-$\alpha$ method.  Journal of Applied Mechanics, 60(2):371–375, 1993.
  * [2] G. M. Hulbert and J. Chung.  Explicit time integration algorithms for structural dynamics with optimal numerical dissipation.  Computer Methods in Applied Mechanics and Engineering, 137(2):175 – 188, 1996.
  * [3] N. M. Newmark.  A method of computation for structural dynamics.  Journal of the engineering mechanics division, 85(3):67–94, 1959\.
  * [4] H. M. Hilber, T. J. R. Hughes, and R. L. Taylor.  Improved numerical dissipation for time integration algorithms in structural dynamics.  Earthquake Engineering & Structural Dynamics, 5(3):283–292, 1977\.
  * [5] J. C. Butcher.  Numerical methods for ordinary differential equations.  John Wiley & Sons, 2016.
  * [6] J. A. Evans, R. R. Hiemstra, T. J. R. Hughes, and A. Reali.  Explicit higher-order accurate isogeometric collocation methods for structural dynamics.  Computer Methods in Applied Mechanics and Engineering, 338:208–240, 2018.
  * [7] G. Loli, G. Sangalli, and M. Tani.  Easy and efficient preconditioning of the isogeometric mass matrix.  Computers & Mathematics with Applications, 2021.
  * [8] T. J. R. Hughes, J. A. Evans, and A. Reali.  Finite element and nurbs approximations of eigenvalue, boundary-value, and initial-value problems.  Computer Methods in Applied Mechanics and Engineering, 272:290–320, 2014.
  * [9] T. J. R. Hughes.  The finite element method: linear static and dynamic finite element analysis.  Courier Corporation, 2012.
  * [10] P. Behnoudfar, Q. Deng, and V. M. Calo.  High-order generalized-alpha method.  Applications in Engineering Science, 4:100021, 2020.
  * [11] P. Behnoudfar, Q. Deng, and V. M. Calo.  Higher-order generalized-$\alpha$ methods for hyperbolic problems.  arXiv preprint arXiv:1906.06081, 2019.
  * [12] E. Hairer and G. Wanner.  Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, volume 14.  Springer, 2010.
  * [13] P. Behnoudfar, Q. Deng, and V. M. Calo.  Split generalized-$\alpha$ method: A linear-cost solver for multi-dimensional second-order hyperbolic systems.  Computer Methods in Applied Mechanics and Engineering, 376:113656, 2021.
  * [14] P. Behnoudfar, V. M. Calo, Q. Deng, and P. D. Minev.  A variationally separable splitting for the generalized-$\alpha$ method for parabolic equations.  International Journal for Numerical Methods in Engineering, 121(5):828–841, 2020.
  * [15] J. A. Cottrell, T. J. R. Hughes, and Y. Bazilevs.  Isogeometric analysis: toward integration of CAD and FEA.  John Wiley & Sons, Chichester, 2009.
  * [16] C. De Boor.  A practical guide to splines, Revised Edition, volume 27 of Applied Mathematical Sciences.  Springer-Verlag, New York, 2001.
  * [17] L. Beirão da Veiga, A. Buffa, G. Sangalli, and R. H. Vázquez.  Mathematical analysis of variational isogeometric methods.  Acta Numer., 23:157––287, 2014.
  * [18] A. Toselli and O. Widlund.  Domain decomposition methods—algorithms and theory, volume 34 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, 2005.
  * [19] L. Gao and V. M. Calo.  Fast isogeometric solvers for explicit dynamics.  Computer Methods in Applied Mechanics and Engineering, 274:19–41, 2014.
  * [20] M. Łoś, M. Paszyński, A. Kłusek, and W. Dzwinel.  Application of fast isogeometric l2 projection solver for tumor growth simulations.  Computer Methods in Applied Mechanics and Engineering, 316:1257–1269, 2017.
  * [21] L. Gao and V. M. Calo.  Preconditioners based on the alternating-direction-implicit algorithm for the 2d steady-state diffusion equation with orthotropic heterogeneous coefficients.  Journal of Computational and Applied Mathematics, 273:274–295, 2015\.
  * [22] R. H. Vázquez.  A new design for the implementation of isogeometric analysis in Octave and Matlab: GeoPDEs 3.0.  Comput. Math. Appl., 72(3):523–554, 2016.