with the reference 1D results. Here, we want to investigate and assess the
contribution to the 0D model accuracy obtained by including only a single
component as nonlinear and to compare these different possible scenarios with
the previously considered ones where all components are either linear or
nonlinear.
To this aim, we performed, for both the 37-artery network and ADAN56 model, a
new set of fully 0D simulations, where in each simulation a single nonlinear
term is activated in the 0D models (either the resistance $R$, the inductance
$L$ or the pressure-area relation), while all the other components are kept
linear. The obtained 0D results are then compared to the available fully-
linear and fully-nonlinear 0D results.

Figures 13 and 14 display the error analysis carried out for the 37-artery
network and ADAN56 model, respectively. For each type of 0D model (fully-
linear, with nonlinear resistance, with nonlinear inductance, with nonlinear
pressure-area relation and fully-nonlinear), the errors in the predicted 0D
results are computed with respect to the reference 1D results in the different
relative error metrics (83). In each boxplot, errors are compared over the
entire network and in several specific vascular regions, namely the aorta,
right arm and right leg, to examine whether specific features of the
performance of the different 0D models emerge in any of the vascular districts
considered. In all boxplots, for all error metrics, 0D model types and
vascular regions, we report the median value, the interquartile range (IQR)
and the minimum and maximum values of the error data. The IQR is computed as
the difference between the 75th (upper) and 25th (lower) quartiles of each
error dataset.

From this error analysis, several observations and outcomes are worthy to be
pointed out. Overall, the fully-nonlinear scenario improves the performance of
the corresponding 0D blood flow models and the accuracy of the predicted
results with respect to the fully-linear scenario; in few isolated cases, the
fully-nonlinear and fully-linear results are comparable, with no significant
gain obtained when all the nonlinearities are included in the 0D models, but
also with no relevant worsening. Furthermore, when activating single nonlinear
terms in the 0D models, the nonlinear inductance $L$ results to be the most
determining nonlinear component which brings the major benefits to these
models in terms of reduction of errors with respect to the reference 1D
solution. Indeed, the contribution of this nonlinear parameter alone is
significant in accurately defining both pressure and flow waveforms, with
errors having small median values and small variability (small-sized IQR). The
nonlinear pressure-area relation also plays a considerable role in defining
more accurate pressure results. This fact is evident when looking at the
boxplots for the relative systolic errors, for both arterial networks
considered. Finally, it is worth mentioning that the main motivation for the
introduction of nonlinear components, especially the pressure-area relation,
regards the application of these 0D models to situations in which large
deviations from a baseline state are to be considered, such as hypertension,
an haemorrhage, a collapsed state or postural changes. Indeed, in these non
physiological conditions, the combined contribution of the nonlinear pressure-
area relation and of the nonlinear pressure-dependent parameters is expected
to play an important role in accurately predicting the corresponding
haemodynamic states.

Figure 13: 37-artery network. Boxplots comparing the errors computed in the
different relative error metrics (83) with respect to the reference 1D results
for (in each boxplot, from left to right): fully-linear 0D results (all-LIN),
0D results with nonlinear resistance (NL-R), 0D results with nonlinear
inductance (NL-L), 0D results with nonlinear pressure-area relation (NL-Tl)
and fully-nonlinear 0D results (all-NL). In each boxplot, errors are compared
over the entire network (yellow) and in several specific vascular regions:
aorta (blue), right arm (green) and right leg (red). Relative errors for
pressure $P$ are shown in the top row boxplots (from left to right boxplots:
RMS, systolic and diastolic errors), while relative errors for flow rate $Q$
are reported in the bottom row boxplots. In each box, the box mid-line
represents the median value, the box itself corresponds to the interquartile
range (IQR), while the lines (whiskers) extending from the box indicate the
minimum and maximum values. Figure 14: ADAN56 model. Boxplots comparing the
errors computed in the different relative error metrics (83) with respect to
the reference 1D results for (in each boxplot, from left to right): fully-
linear 0D results (all-LIN), 0D results with nonlinear resistance (NL-R), 0D
results with nonlinear inductance (NL-L), 0D results with nonlinear pressure-
area relation (NL-Tl) and fully-nonlinear 0D results (all-NL). In each
boxplot, errors are compared in the entire network (yellow) and in several
specific vascular regions: aorta (blue), right arm (green) and right leg
(red). Relative errors for pressure $P$ are shown in the top row boxplots
(from left to right boxplots: RMS, systolic and diastolic errors), while
relative errors for flow rate $Q$ are reported in the bottom row boxplots. In
each box, the box mid-line represents the median value, the box itself
corresponds to the interquartile range (IQR), while the lines (whiskers)
extending from the box indicate the minimum and maximum values.

##  6 Concluding remarks and future work

In this paper we have extended the well-known family of lumped-parameter
models for blood flow in a vascular segment, in order to construct 0D models
that preserve important properties of the original 1D model. In particular,
the main features characterizing these models are: (i) the nonlinearity in the
pressure-area relationship and (ii) the nonlinearity in the model parameters
$R$ and $L$, which depend on the average time-varying cross-sectional area.
The different 0D vessel configurations depending on the data prescribed at the
inlet and outlet of the vessel have been described; 0D junctions have been
introduced to couple 0D vessels in a network.
In order to validate the newly derived family of nonlinear 0D models, fully 1D
and fully 0D simulations have been performed and compared for three different
benchmark arterial networks, showing that, even if there are discrepancies
between 0D and 1D results, these nonlinear 0D models are able to capture and
reproduce the essential features of pressure and flow waveforms. Furthermore,
remarkable differences between linear and nonlinear 0D results have been
observed in these physiologically normal cases and we expect that these
differences will be even more relevant when addressing large deviations from
the baseline haemodynamic states, like, for example, when simulating
hypertension or haemorrhagic events.
Our ultimate goal is to construct hybrid 1D-0D networks, where certain vessels
are treated 1D, while others are modelled as 0D vessels. Indeed, our work is
primarily motivated by the attempt of facing the issues of computational
efficiency and execution time, related to the modelling of blood flow in
highly complex networks. Potential applications of the proposed framework will
regard for instance the ADAN model developed by Blanco et al. [4, 5] and the
global, closed-loop, multiscale mathematical model for the human
cardiovascular system of Müller and Toro [25, 26]. By approximating not
necessarily all, but some vessels of such complex networks with these
nonlinear 0D models, defining the different 0D vessel configurations, we
expect that the computational cost and execution time of the resulting hybrid
1D-0D vessel networks will be dramatically reduced with respect to the fully
1D networks, even for long term simulations, and that the results will be
significantly improved with respect to the fully 0D networks, still preserving
the topology of the original 1D network and without losing any essential
information with respect to the full 1D model. For this purpose, future work
will include to find relevant a-priori model selection criteria to define an
adaptive model selection strategy that will allow to determine, given a
network, the most suitable model to be used for each vessel of the network,
either the 1D model or the newly derived nonlinear 0D models.
Moreover, we will consider purely 0D and hybrid 1D-0D couplings in order to
find appropriate coupling conditions between 1D/0D vessels that do not require
any restrictions on the admissible 0D vessel types converging at a junction.
Finally, in the proposed nonlinear lumped-parameter models we will incorporate
additional aspects of the pressure-area relation, such as the viscoelastic
properties of vessel walls, in order to derive more realistic and reliable 0D
models.

## Acknowledgements

Beatrice Ghitti acknowledges the University of Trento for financing her Ph.D.
studentship.
Lucas O. Müller acknowledges funding from the Italian Ministry of Education,
University and Research (MIUR) in the frame of the Departments of Excellence
Initiative 2018–2022 attributed to the Department of Mathematics of the
University of Trento (grant L. 232/2016) and in the frame of the PRIN 2017
project Innovative numerical methods for evolutionary partial differential
equations and applications.
All authors warmly thank Prof. Marco Sabatini (Department of Mathematics,
University of Trento, Italy) for very helpful discussions that contributed to
the completion of the ODE stability analysis presented in this work.

## References

  * [1] J Alastruey, KH Parker, J Peiró, and SJ Sherwin.  Lumped parameter outflow models for 1-D blood flow simulations: Effect on pulse waves and parameter estimation.  Communications in Computational Physics, 4:317–336, 2008.
  * [2] D. S. Berger and J. K. J. Li.  Temporal relationship between left ventricular and arterial system elastances.  IEEE Transactions on Biomedical Engineering, 39(4):404–410, April 1992.
  * [3] J. Blacher and M. E. Safar.  Large-artery stiffness, hypertension and cardiovascular risk in older patients.  Nature Clinical Practice Cardiovascular Medicine, 2:450–5, October 2005.
  * [4] P. J. Blanco, S. M. Watanabe, E. Dari, M. A. R. Passos, and R. A. Feijóo.  Blood flow distribution in an anatomically detailed arterial network model: Criteria and algorithms.  Biomechanics and Modeling in Mechanobiology, 13:1303–1330, March 2014.
  * [5] PJ Blanco, S Watanabe, M Passos, P Lemos, and RA Feijóo.  An anatomically detailed arterial network model for one-dimensional computational hemodynamics.  IEEE transactions on bio-medical engineering, 62:736–53, 2015.
  * [6] Etienne Boileau, P. Nithiarasu, Pablo J. Blanco, Lucas O. Müller, Fredrik Eikeland Fossan, L. R. Hellevik, W. P. Donders, W. Huberts, M. Willemet, and J. Alastruey.  A benchmark study of numerical schemes for one-dimensional arterial blood flow modelling.  International Journal for Numerical Methods in Biomedical Engineering, 31(10), 2015.
  * [7] A. Cappello, G. Gnudi, and C. Lamberti.  Identification of the three-element windkessel model incorporating a pressure-dependent compliance.  Annals of Biomedical Engineering, 23:164–177, 2006.
  * [8] K. Cruickshank, L. Riste, S. G. Anderson, J. S. Wright, G. Dunn, and R. G. Gosling.  Aortic pulse-wave velocity and its relationship to mortality in diabetes and glucose intolerance.  Circulation, 106(16):2085–2090, 2002.
  * [9] S. Epstein, M. Willemet, P. J. Chowienczyk, and J. Alastruey.  Reducing the number of parameters in 1D arterial blood flow modeling: Less is more for patient-specific simulations.  American Journal of Physiology-Heart and Circulatory Physiology, 309(1):H222–H234, 2015.
  * [10] R. Fogliardi, M. Di Donfrancesco, and R. Burattini.  Comparison of linear and nonlinear formulations of the three-element windkessel model.  American Journal of Physiology-Heart and Circulatory Physiology, 271(6):H2661–H2668, 1996.
  * [11] L. Formaggia, A. Quarteroni, and A. Veneziani, editors.  Cardiovascular Mathematics: Modeling and Simulation of the Circulatory System.  Number 1 in MS & A : Modeling, Simulation & Applications. Springer, Milano, 2009.
  * [12] L. Formaggia and A. Veneziani.  Reduced and multiscale models for the human cardiovascular system.  Technical report, Politecnico di Milano, October 2015.
  * [13] F. E. Fossan, J. Mariscal-Harana, J. Alastruey, and L. R. Hellevik.  Optimization of topological complexity for one-dimensional arterial blood flow models.  Journal of The Royal Society Interface, 15(149):20180546, 2018.
  * [14] S. Fujimoto, R. Mizuno, Y. Saito, and S. Nakamura.  Clinical application of wave intensity for the treatment of essential hypertension.  Heart and Vessels, 19(1):19–22, January 2004.
  * [15] Y. C. Fung.  Biomechanics: Mechanical Properties of Living Tissues.  Springer, second edition, 1993.  Second Edition.
  * [16] A. Ghigo.  Reduced-Order Models for Blood Flow in Networks of Large Arteries.  PhD thesis, Université Pierre et Marie Curie, Paris, September 2017\.
  * [17] J. K. Hale.  Ordinary Differential Equations.  John Wiley & Sons, Inc., 1969.
  * [18] A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy.  Uniformly High Order Accuracy Essentially Non–oscillatory Schemes, III.  Journal of Computational Physics, 71:231–303, 1987.
  * [19] A. Harten and S. Osher.  Uniformly High–Order Accurate Nonoscillatory Schemes. I.  SIAM Journal on Numerical Analysis, 24(2):279–309, 1987.
  * [20] P. J. Hunter.  Numerical simulation of arterial blood flow.  Master’s thesis, The University of Auckland, Auckland, 1972.
  * [21] J. K. J. Li, T. Cui, and G. M. Drzewiecki.  A nonlinear model of the arterial system incorporating a pressure-dependent compliance.  IEEE Transactions on Biomedical Engineering, 37(7):673–678, 1990\.
  * [22] K. S. Matthys, J. Alastruey, J. Peiró, A. W. Khir, P. Segers, P. R. Verdonck, K. H. Parker, and S. J. Sherwin.  Pulse wave propagation in a model human arterial network: Assessment of 1-D numerical simulations against in vitro measurements.  Journal of Biomechanics, 40(15):3476–3486, 2007.
  * [23] V. Milišić and A. Quarteroni.  Analysis of lumped parameter models for blood flow simulations and their relation with 1D models.  ESAIM: Mathematical Modelling and Numerical Analysis, 38(4):613–632, 2004.
  * [24] M. Mirramezani and S. C. Shadden.  A Distributed Lumped Parameter Model of Blood Flow.  Annals of Biomedical Engineering, 48(12):2870–2886, 2020.
  * [25] L. O. Müller and E. F. Toro.  A global multiscale mathematical model for the human circulation with emphasis on the venous system.  International Journal for Numerical Methods in Biomedical Engineering, 30(7):681–725, 2014.
  * [26] Lucas O. Müller and Eleuterio F. Toro.  Enhanced global mathematical model for studying cerebral venous blood flow.  Journal of Biomechanics, 47(13):3361–3372, October 2014.
  * [27] JP Murgo, N Westerhof, JP Giolma, and SA Altobelli.  Aortic input impedance in normal man: Relationship to pressure wave forms.  Circulation, 62:105–116, 1980.
  * [28] J. P. Mynard.  Computer Modelling and Wave Intensity Analysis of Perinatal Cardiovascular Function and Dysfunction.  PhD thesis, Department of Paediatrics, The University of Melbourne, August 2011.
  * [29] J. P. Mynard and J. J. Smolich.  One-Dimensional Haemodynamic Modeling and Wave Dynamics in the Entire Adult Circulation.  Annals of Biomedical Engineering, 43(6):1443–1460, 2015.
  * [30] S. Safaei, P. J. Blanco, L. O. Müller, F. E. Fossan, L. R. Hellevik, and P. J. Hunter.  Bond graph model of cerebral circulation: Toward clinically feasible systemic blood flow simulations.  Frontiers in Physiology, 9:148, 2018.
  * [31] K. Sagawa, R. K. Lie, and J. Schaefer.  Translation of Otto Frank’s paper “Die Grundform des Arteriellen Pulses” Zeitschrift für Biologie 37: 483–526 (1899).  Journal of Molecular and Cellular Cardiology, 22(3):253–254, 1990\.
  * [32] M. Saito, Y. Ikenaga, M. Matsukawa, Y. Watanabe, T. Asada, and P.-Y. Lagrée.  One-dimensional model for propagation of a pressure wave in a model of the human arterial network: Comparison of theoretical and experimental results.  Journal of biomechanical engineering, 133:121005, December 2011\.
  * [33] D. A. Sánchez.  Ordinary Differential Equations and Stability Theory.  Dover Publications, Inc., 1968.
  * [34] S. J. Sherwin, V. Franke, J. Peiró, and K. H. Parker.  One-dimensional modelling of a vascular network in space-time variables.  Journal of Engineering Mathematics, 47:217–250, January 2003.
  * [35] Y. Shi, P. Lawford, and R. Hose.  Review of Zero-D and 1-D Models of Blood Flow in the Cardiovascular System.  BioMedical Engineering OnLine, 10(1):33, 2011.
  * [36] A. Spilimbergo, E. F. Toro, and L. O. Müller.  One-Dimensional Blood Flow with Discontinuous Properties and Transport: Mathematical Analysis and Numerical Schemes.  Communications in Computational Physics, 29(3):649–697, 2021.
  * [37] E. F. Toro.  Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction.  Springer-Verlag Berlin Heidelberg, 2009.  Third Edition.
  * [38] E. F. Toro.  Brain venous haemodynamics, neurological diseases and mathematical modelling. A review.  Applied Mathematics and Computation, 272:542–579, 2016.
  * [39] E. F. Toro, R. C. Millington, and L. A. M. Nejad.  Towards Very High–Order Godunov Schemes.  In Godunov Methods: Theory and Applications. Edited Review, e. F. Toro (Editor), pages 905–937. Kluwer Academic/Plenum Publishers, 2001.
  * [40] M. Ursino.  A mathematical model of the carotid baroregulation in pulsating conditions.  IEEE Transactions on Biomedical Engineering, 46(4):382–392, April 1999.
  * [41] M. Ursino, A. Fiorenzi, and E. Belardinelli.  The role of pressure pulsatility in the carotid baroreflex control: A computer simulation study.  Computers in Biology and Medicine, 26(4):297–314, 1996.
  * [42] M Ursino and C. A. Lodi.  A simple mathematical model of the interaction between intracranial pressure and cerebral hemodynamics.  Journal of Applied Physiology, 82:1256–1269, 1997.
  * [43] B. van Leer.  On the Relation between the Upwind–Differencing Schemes of Godunov, Engquist-Osher and Roe.  SIAM Journal on Scientific and Statistical Computing, 5(1):1–20, 1985.
  * [44] N. Xiao, J. Alastruey, and C. A. Figueroa.  A systematic comparison between 1-D and 3-D hemodynamics in compliant arterial models.  International Journal for Numerical Methods in Biomedical Engineering, 30(2):204–231, February 2014.

##  Appendix A Stability analysis

###  A.1 First-order ordinary differential equations with periodic forcing
function.

It is a well-known result that, given a linear homogeneous system of ODEs with
constant coefficients, that is

$\frac{d\bm{x}(t)}{dt}=\bm{A}\bm{x}(t),$ (87)

the stability of the exact solution of such ODE system is determined by the
real part of the eigenvalues of the coefficient matrix $\bm{A}$. Precisely, a
necessary and sufficient condition for the ODE system (87) to be
asymptotically stable is that all eigenvalues of $\bm{A}$ have strictly
negative real part. If this is the case, then there exist positive constants
$\alpha$, $\beta$ such that

$\lVert e^{\bm{A}t}\rVert\leq\beta e^{-\alpha t},\quad t\geq 0.$ (88)

Proof of this result can be found in classical books on ordinary differential
equations, see for instance [33, 17].
Here, we want to prove that, when a periodic forcing function $\bm{b}(t)$ is
added to system (87) to obtain the inhomogeneous ODE system (52), if the
homogeneous part is asymptotically stable, then any solution of the complete
non-homogeneous system will converge to a periodic solution as $t$ increases.
Namely, the periodic solution of the inhomogeneous linear problem (52) is an
attracting stable solution for the system.
The exact solution of system (52) obtained with the method of variation of
parameters reads

$\bm{x}(t)=e^{\bm{A}t}\bm{x}_{0}+e^{\bm{A}t}\int_{0}^{t}e^{-\bm{A}s}\bm{b}(s)\,ds=e^{\bm{A}t}\bm{x}_{0}+\int_{0}^{t}e^{\bm{A}(t-s)}\bm{b}(s)\,ds,$
(89)

where $\bm{x}(0)=\bm{x}_{0}$ and $e^{\bm{A}t}$ is the matrix exponential.
Under the assumption of periodic forcing function $\bm{b}(t)$, i.e.
$\bm{b}(T_{0})=\bm{b}(0)$ for some period $T_{0}>0$, we are able to obtain a
particular solution of (52), that is its periodic solution satisfying
$\bm{x}(T_{0})=\bm{x}(0)$, by choosing appropriately the initial condition
$\bm{x}_{0}$. By enforcing this last condition, after straightforward
calculations, we find that the initial condition $\bm{x}_{0}^{P}$ which
provides a periodic solution of system (52) is given by

$\bm{x}_{0}^{P}=\left(\bm{I}-e^{\bm{A}T_{0}}\right)^{-1}e^{\bm{A}T_{0}}\int_{0}^{T_{0}}e^{-\bm{A}s}\bm{b}(s)\,ds.$
(90)

At this point, we want to show that, if we start from any initial position
different from the one which produces the periodic solution, the values of the
corresponding solution will converge to the values of the periodic solution as
$t$ increases. Let $\bm{x}^{P}$ be the periodic solution corresponding to the
initial condition $\bm{x}_{0}^{P}$ and let $\bm{x}^{NP}$ be any other solution
of system (52) associated to some initial condition $\bm{x}_{0}^{NP}$. Then,
by computing the norm of the difference between the two solutions and applying
(88), the following inequality holds

$\lVert\bm{x}^{P}(t)-\bm{x}^{NP}(t)\rVert=\lVert
e^{\bm{A}t}\left(\bm{x}_{0}^{P}-\bm{x}_{0}^{NP}\right)\rVert\leq\beta
e^{-\alpha t}\lVert\bm{x}_{0}^{P}-\bm{x}_{0}^{NP}\rVert,$ (91)

from which it clearly follows that $\beta e^{-\alpha t}\rightarrow 0$ as
$t\rightarrow+\infty$, that is the sought result. In conclusion, if the
forcing function $\bm{b}(t)$ is periodic and if the homogeneous part of system
(52) is asymptotically stable, then the exact solution of the original
inhomogeneous problem will converge to the periodic solution of system (52) as
$t$ increases, for any admissible choice of the initial condition
$\bm{x}(0)=\bm{x}_{0}=\left[V(0),Q(0)\right]^{T}$. As a consequence, the
stability of the exact solution of the complete ODE system (52) is determined
by the real part of the eigenvalues of the coefficient matrix $\bm{A}$.

###  A.2 First-order system of ordinary differential equations with null
eigenvalue.

Here, as discussed is Section 3, we derive the additional assumption on the
periodic forcing function $\bm{b}(t)$, which is needed to preserve the
stability of the solution of an inhomogeneous ODE system of the form (52)
whose coefficient matrix $\bm{A}$ has a null eigenvalue (and all other
eigenvalues with strictly negative real part). This assumption will ensure
that the stability of the homogeneous system (which is stable, but not
asymptotically stable) is inherited also by the corresponding inhomogeneous
system, namely that any solution of (52) for any admissible choice of the
initial condition will converge to the periodic one as $t\to+\infty$.
This condition on the periodic forcing function is first derived in the scalar
case for a single equation; then, the obtained result is extended to the case
of a system of ODEs, specifically to system (63) governing the linear
$(Q_{in},Q_{out})$-type 0D vessel.

####  A.2.1 Scalar case.

Let us start from the scalar case of a single ODE and consider the following
mass conservation equation

$\dfrac{dV}{dt}=Q_{in}-Q_{out}.$ (92)

For a scalar ODE, the coefficient matrix $\bm{A}$ in (52) is reduced to a
single constant coefficient $a$, which in the case of the above equation is
equal to zero and plays the role of the null eigenvalue of matrix $\bm{A}$. As
a consequence, the right-hand side of equation (92) only depends on the
scalar, periodic forcing function $b(t)=Q_{in}(t)-Q_{out}(t)$, where
$Q_{in}(t)$ and $Q_{out}(t)$ represent the inflow entering the vessel and the
outflow leaving the vessel, respectively, and are both periodic signals with
the same period $T_{0}>0$. The exact solution of equation (92) is

$V(t)=V(0)+\int_{0}^{t}b(s)\,ds,$ (93)

with initial condition $V(0)$. At this point, we are interested in studying
the asymptotic properties of such a solution, namely the behaviour of $V(t)$
as $t$ approaches $+\infty$.
Let $t=nT_{0}+\tau$, with $n\in\mathbb{N}$ sufficiently large and
$0\leq\tau<T_{0}$. Then, we have

$\begin{split}V(t)=V(nT_{0}+\tau)&=V(0)+\int_{0}^{nT_{0}+\tau}b(s)\,ds\\\
&=V(0)+\int_{0}^{nT_{0}}b(s)\,ds+\int_{nT_{0}}^{nT_{0}+\tau}b(s)\,ds\\\
&=V(0)+n\int_{0}^{T_{0}}b(s)\,ds+\int_{0}^{\tau}b(s)\,ds,\end{split}$ (94)

where, in the last equality, we have used the periodicity of the forcing term
$b(t)$. In the above expression of the exact solution $V(t)$, the initial
volume $V(0)$ is a fixed finite quantity, the second integral,
$\int_{0}^{\tau}b(s)\,ds,$ is also finite, since it is the integral of a
continuous periodic function over a fraction of its period $[0,T_{0}]$.
Therefore, studying the asymptotic behaviour of the solution $V(t)$, that is

$\lim_{t\to+\infty}V(t)$ (95)

is completely equivalent to study

$\lim_{n\to+\infty}n\int_{0}^{T_{0}}b(s)\,ds.$ (96)

We distinguish the following three cases:

$\lim_{n\to+\infty}n\int_{0}^{T_{0}}b(s)\,ds=\begin{cases}+\infty&\text{if
}\overline{b}>0,\\\ -\infty&\text{if }\overline{b}<0,\\\ 0&\text{if
}\overline{b}=0,\\\ \end{cases}$ (97)

where we have set

$\overline{b}:=\int_{0}^{T_{0}}b(s)\,ds.$ (98)

In other words, if the integral average of the periodic forcing function
$b(t)$ over its period $[0,T_{0}]$ is different from zero, either positive or
negative, then the solution of equation (92) will diverge either to $+\infty$
or to $-\infty$, depending on the sign of $\overline{b}$, as $t\to+\infty$. On
the other hand, if the integral average of $b(t)$ over its period is zero,
then any solution $V(t)$ of (92) will be bounded, namely

$\lim_{t\to+\infty}\big{|}V(t)\big{|}<+\infty\quad\Longleftrightarrow\quad\overline{b}=0\quad\Longleftrightarrow\quad\dfrac{1}{T_{0}}\int_{0}^{T_{0}}b(s)\,ds=0,$
(99)

and will converge to the periodic one, that is

$V(nT_{0}+\tau)\xrightarrow[]{n\to+\infty}V(0)+\int_{0}^{\tau}b(s)\,ds.$ (100)

We point out that this result is also physically consistent. Indeed, for
equation (92) the condition stated in (99) simply becomes

$\dfrac{1}{T_{0}}\int_{0}^{T_{0}}\left[Q_{in}(s)-Q_{out}(s)\right]\,ds=0\quad\Longleftrightarrow\quad\int_{0}^{T_{0}}Q_{in}(s)\,ds=\int_{0}^{T_{0}}Q_{out}(s)\,ds,$
(101)

which means that, in order for the volume $V(t)$ not to constantly
increase/decrease and asymptotically explode, the inflow entering the vessel
and the outflow leaving the vessel must be balanced over each period.

####  A.2.2 ODE system case.

Let us now extend the previous result to the case of a system of ODEs,
specifically to system (63). For any admissible choice of the initial
conditions $\bm{x}_{0}=\left[V(0),Q(0),V_{d}(0)\right]^{T}$, the general
solution of (63) is given by the formula in (89), where the constant
coefficient matrix $\bm{A}$ and the periodic forcing function $\bm{b}(t)$ are

$\bm{A}=\left[\begin{matrix}0&-1&0\\\\[7.74997pt]
\dfrac{2}{C_{0}L_{0}}&-\dfrac{R_{0}}{L_{0}}&-\dfrac{2}{C_{0}L_{0}}\\\\[7.74997pt]
0&1&0\end{matrix}\right],\quad\bm{b}(t)=\left[\begin{matrix}Q_{in}(t)\\\\[7.74997pt]
0\\\\[7.74997pt] -Q_{out}(t)\end{matrix}\right].$ (102)

In general, if an $n\times n$ matrix $\bm{A}$ is diagonalizable, with linearly
independent eigenvectors $\bm{v}_{1},\bm{v}_{2},\ldots,\bm{v}_{n}$ and
corresponding eigenvalues $\lambda_{1},\lambda_{2},\ldots,\lambda_{n}$, then,
denoted by $\bm{S}$ the matrix whose columns are the eigenvectors of $\bm{A}$,
we have

$\bm{S}^{-1}\bm{A}\bm{S}=\left[\begin{matrix}\lambda_{1}&0&\cdots&0\\\
0&\lambda_{2}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&\lambda_{n}\end{matrix}\right]=\bm{D}.$ (103)

The exponential of the diagonal matrix $\bm{D}$ is easily written as

$e^{\bm{D}}=\left[\begin{matrix}e^{\lambda_{1}}&0&\cdots&0\\\
0&e^{\lambda_{2}}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&e^{\lambda_{n}}\end{matrix}\right]$ (104)

and thus, it is straightforward to check that, under the assumption that
$\bm{A}$ is diagonalizable, the corresponding matrix exponential $e^{\bm{A}t}$
can be computed as follows

$e^{\bm{A}t}=\bm{S}\left[\begin{matrix}e^{\lambda_{1}t}&0&\cdots&0\\\
0&e^{\lambda_{2}t}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
0&0&\cdots&e^{\lambda_{n}t}\end{matrix}\right]\bm{S}^{-1}=\bm{S}e^{\bm{D}t}\bm{S}^{-1}.$
(105)

Then, by introducing the above result in (89), the exact solution of an
inhomogeneous ODE system of the form (52) with diagonalizable coefficient
matrix $\bm{A}$ can be reformulated as

$\bm{x}(t)=\left(\bm{S}e^{\bm{D}t}\bm{S}^{-1}\right)\bm{x}_{0}+\int_{0}^{t}\left(\bm{S}e^{\bm{D}(t-s)}\bm{S}^{-1}\right)\bm{b}(s)\,ds,$
(106)

As observed at the beginning of this paragraph, if the coefficient matrix
$\bm{A}$ has a zero eigenvalue and all other eigenvalues with strictly
negative real part, then the homogeneous part of system (52) is stable, but
not asymptotically stable, with an attractive line of equilibrium points.
Therefore, to study the asymptotic properties of the exact solution of the
inhomogeneous system (52) obtained by introducing a periodic forcing function
$\bm{b}(t)$, we have to focus on the integral term in (106), to study its
behaviour as $t$ tends to $+\infty$. We set

$\widetilde{\bm{S}}(t):=\bm{S}e^{\bm{D}t}\bm{S}^{-1}.$ (107)

By restricting to the order $n=3$ of the coefficient matrix of system (63),
after straightforward calculations, we get that the $i$th component of the
integrand vector function in (106), obtained by multiplying the $i$th row of
matrix $\widetilde{\bm{S}}$ by the vector forcing function $\bm{b}(t)$ and
using the fact that $\lambda_{1}=0$, is given by

$\begin{split}(\widetilde{\bm{S}}(t-s))_{i}\bm{b}(s)&=\sum_{j=1}^{3}(\widetilde{\bm{S}}(t-s))_{ij}\bm{b}_{j}(s)\\\
&=\dfrac{1}{|\bm{S}|}\left[\left(\bm{v}_{22}\bm{v}_{33}-\bm{v}_{23}\bm{v}_{32}\right)\bm{v}_{1i}+\left(\bm{v}_{23}\bm{v}_{31}-\bm{v}_{21}\bm{v}_{33}\right)\bm{v}_{2i}e^{\lambda_{2}(t-s)}+\left(\bm{v}_{21}\bm{v}_{32}-\bm{v}_{22}\bm{v}_{31}\right)\bm{v}_{3i}e^{\lambda_{3}(t-s)}\right]\bm{b}_{1}(s)\\\
&+\dfrac{1}{|\bm{S}|}\left[\left(\bm{v}_{13}\bm{v}_{32}-\bm{v}_{12}\bm{v}_{33}\right)\bm{v}_{1i}+\left(\bm{v}_{11}\bm{v}_{33}-\bm{v}_{13}\bm{v}_{31}\right)\bm{v}_{2i}e^{\lambda_{2}(t-s)}+\left(\bm{v}_{12}\bm{v}_{31}-\bm{v}_{11}\bm{v}_{32}\right)\bm{v}_{3i}e^{\lambda_{3}(t-s)}\right]\bm{b}_{2}(s)\\\
&+\dfrac{1}{|\bm{S}|}\left[\left(\bm{v}_{12}\bm{v}_{23}-\bm{v}_{13}\bm{v}_{22}\right)\bm{v}_{1i}+\left(\bm{v}_{13}\bm{v}_{21}-\bm{v}_{11}\bm{v}_{23}\right)\bm{v}_{2i}e^{\lambda_{2}(t-s)}+\left(\bm{v}_{11}\bm{v}_{22}-\bm{v}_{12}\bm{v}_{21}\right)\bm{v}_{3i}e^{\lambda_{3}(t-s)}\right]\bm{b}_{3}(s)\end{split}$
(108)

with $|\bm{S}|=\text{det}(\bm{S})$, for $i=1,\ldots,3$. As both eigenvalues
$\lambda_{2,3}$ have strictly negative real part, the presence of the
exponential factors $e^{\lambda_{2}t}$ and $e^{\lambda_{3}t}$ ensures the
convergence of the corresponding terms as $t\to+\infty$, independently of any
assumption on the periodic forcing function $\bm{b}(t)$. On the other hand,
since $e^{\lambda_{1}t}=e^{0}=1$ for any $t\geq 0$, when integrating the
corresponding terms without decaying exponential factor, namely

$\int_{0}^{t}\dfrac{\bm{v}_{1i}}{|\bm{S}|}\left[\left(\bm{v}_{22}\bm{v}_{33}-\bm{v}_{23}\bm{v}_{32}\right)\bm{b}_{1}(s)+\left(\bm{v}_{13}\bm{v}_{32}-\bm{v}_{12}\bm{v}_{33}\right)\bm{b}_{2}(s)+\left(\bm{v}_{12}\bm{v}_{23}-\bm{v}_{13}\bm{v}_{22}\right)\bm{b}_{3}(s)\right]\,ds$
(109)

for $i=1,\ldots,3$, we need to introduce additional assumptions on the
periodic forcing function $\bm{b}(t)$ in order for the above integral terms,
and hence the solution of system (63), not to explode as $t\to+\infty$.
As done in the scalar case, let $t=nT_{0}+\tau$, with $n\in\mathbb{N}$
sufficiently large and $0\leq\tau<T_{0}$. Then, using explicit expressions of
the components of $\bm{b}(t)$ given in (102), we have that the integral in
(109) can be rewritten as follows

$\dfrac{\bm{v}_{1i}}{|\bm{S}|}\int_{0}^{nT_{0}+\tau}\left[\alpha_{1}Q_{in}(s)-\alpha_{2}Q_{out}(s)\right]\,ds,\quad\text{for
}i=1,\dots,3,$ (110)

where we have set

$\alpha_{1}:=\left(\bm{v}_{22}\bm{v}_{33}-\bm{v}_{23}\bm{v}_{32}\right),\quad\alpha_{2}:=\left(\bm{v}_{12}\bm{v}_{23}-\bm{v}_{13}\bm{v}_{22}\right).$
(111)

If we discard the multiplicative factor $\frac{\bm{v}_{1i}}{|\bm{S}|}$, where
$\bm{v}_{1i}$ is the $i$th element of the eigenvector $\bm{v}_{1}$ associated
to the null eigenvalue $\lambda_{1}$, for all the three components of the
solution of system (63) we obtain the same integral, which, by exploiting the
periodicity of $Q_{in}(t)$ and $Q_{out}(t)$, becomes

$\begin{split}\int_{0}^{nT_{0}+\tau}\left[\alpha_{1}Q_{in}(s)-\alpha_{2}Q_{out}(s)\right]\,ds&=n\int_{0}^{T_{0}}\left[\alpha_{1}Q_{in}(s)-\alpha_{2}Q_{out}(s)\right]\,ds\\\
&+\int_{0}^{\tau}\left[\alpha_{1}Q_{in}(s)-\alpha_{2}Q_{out}(s)\right]\,ds.\end{split}$
(112)

The second integral in the above expression is finite, since it is the
integral of the linear combination of two continuous, periodic functions, with
the same period $[0,T_{0}]$, over a fraction of their period. Then, by means
of the same procedure used in the previous scalar case, also here we can
observe that, in order for each component of the solution of system (63) not
to diverge either to $+\infty$ or to $-\infty$ as $t\to+\infty$, namely

$\lim_{t\to+\infty}|\bm{x}_{i}(t)|<+\infty\quad\text{for }i=1,\ldots,3,$ (113)

and to converge to the periodic solution, a necessary and sufficient condition
is

$\int_{0}^{T_{0}}\left[\alpha_{1}Q_{in}(s)-\alpha_{2}Q_{out}(s)\right]\,ds=0\quad\Longleftrightarrow\quad\alpha_{1}\int_{0}^{T_{0}}Q_{in}(s)\,ds=\alpha_{2}\int_{0}^{T_{0}}Q_{out}(s)\,ds.$
(114)

We recall that this result is strictly valid only in the case where the
coefficient matrix $\bm{A}$ is diagonalizable and, from linear algebra, it is
well-known that an $n\times n$ matrix is diagonalizable if it has $n$
independent eigenvectors, which is true, for example, when the matrix has $n$
distinct eigenvalues.
It is easy to check that the matrix $\bm{A}$ of system (63) given in (102) is
diagonalizable in both cases $\Delta>0$ and $\Delta<0$, described in (64). If
$\Delta>0$, matrix $\bm{A}$ has three real and distinct eigenvalues

$\lambda_{1}=0,\quad\lambda_{2,3}=-\dfrac{R_{0}}{2L_{0}}\pm\dfrac{\sqrt{\Delta}}{2}$
(115)

and corresponding set of linearly independent eigenvectors

$\bm{v}_{1}=\left[\begin{matrix}1\\\\[5.81247pt] 0\\\\[5.81247pt]
1\end{matrix}\right],\quad\bm{v}_{2}=\left[\begin{matrix}1\\\
\left(\dfrac{R_{0}}{2L_{0}}+\dfrac{\sqrt{\Delta}}{2}\right)\\\
-1\end{matrix}\right],\quad\bm{v}_{3}=\left[\begin{matrix}1\\\
\left(\dfrac{R_{0}}{2L_{0}}-\dfrac{\sqrt{\Delta}}{2}\right)\\\
-1\end{matrix}\right];$ (116)

while, if $\Delta<0$, matrix $\bm{A}$ has three distinct eigenvalues, one
equal to zero and the other two complex and conjugate, given by

$\lambda_{1}=0,\quad\lambda_{2,3}=-\dfrac{R_{0}}{2L_{0}}\pm
i\dfrac{\sqrt{-\Delta}}{2}$ (117)

and corresponding set of linearly independent eigenvectors

$\bm{v}_{1}=\left[\begin{matrix}1\\\\[5.81247pt] 0\\\\[5.81247pt]
1\end{matrix}\right],\quad\bm{v}_{2}=\left[\begin{matrix}1\\\
\left(\dfrac{R_{0}}{2L_{0}}+i\dfrac{\sqrt{-\Delta}}{2}\right)\\\
-1\end{matrix}\right],\quad\bm{v}_{3}=\left[\begin{matrix}1\\\
\left(\dfrac{R_{0}}{2L_{0}}-i\dfrac{\sqrt{-\Delta}}{2}\right)\\\
-1\end{matrix}\right].$ (118)

Equipped with explicit expressions of the eigenvectors of the coefficient
matrix $\bm{A}$, we can then compute factors $\alpha_{1}$ and $\alpha_{2}$
given in (111) and defining the linear combination between the inflow function
$Q_{in}(t)$ and the outflow function $Q_{out}(t)$ in the integral (114). We
obtain

$\alpha_{1}=\alpha_{2}=\begin{cases}-\sqrt{\Delta},&\text{if }\Delta>0,\\\
-i\sqrt{-\Delta},&\text{if }\Delta<0,\\\ \end{cases}$ (119)

so that condition (114) on the periodic forcing function $\bm{b}(t)$ simply
becomes

$\int_{0}^{T_{0}}\left[Q_{in}(s)-Q_{out}(s)\right]\,ds=0\quad\Longleftrightarrow\quad\int_{0}^{T_{0}}Q_{in}(s)\,ds=\int_{0}^{T_{0}}Q_{out}(s)\,ds,$
(120)

which is exactly the same physically consistent condition found in the scalar
case to get asymptotically the periodic solution.