Fourth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (4th-CASAM-N): I. Mathematical Framework

Abstract

This work presents the fourth-order comprehensive sensitivity analysis methodology for nonlinear systems (abbreviated as “4th-CASAM-N”) for exactly and efficiently computing the first-, second-, third-, and fourth-order functional derivatives (customarily called “sensitivities”) of physical system responses (i.e., “system performance parameters”) to the system’s (or model) parameters. The qualifier “comprehensive” indicates that the 4th-CASAM-N methodology enables the exact and efficient computation not only of response sensitivities with respect to the customary model parameters (including computational input data, correlations, initial and/or boundary conditions) but also with respect to imprecisely known material boundaries, caused by manufacturing tolerances, of the system under consideration. The 4th-CASAM-N methodology presented in this work enables the hitherto very difficult, if not intractable, exact computation of all of the first-, second-, third-, and fourth-order response sensitivities for large-scale systems involving many parameters, as usually encountered in practice. Notably, the implementation of the 4th-CASAM-N requires very little additional effort beyond the construction of the adjoint sensitivity system needed for computing the first-order sensitivities. The application of the principles underlying the 4th-CASAM-N to an illustrative paradigm nonlinear heat conduction model will be presented in an accompanying work.

1. Introduction

The computational model of a physical system comprises the following conceptual components: (a) a well-posed system of that relate the system’s independent variables and parameters to the system’s state (i.e., dependent) variables; (b) probability distributions, moments thereof, inequality and/or equality constraints that define the range of variations of the system’s parameters; and (c) one or several quantities, customarily referred to as system responses (or objective functions, or indices of performance), which are computed using the mathematical model. This works presents a new, general-purpose methodology for computing exactly and efficiently functional derivatives (called “sensitivities”) of results (“system responses”), predicted by nonlinear mathematical models of systems (physical, engineering, biological) involving imprecisely known (i.e., uncertain) parameters, including input data, correlations, initial and/or boundary conditions, as well as manufacturing tolerances that affect domain of the model’s definition in phase space. This new method is called the “fourth-order comprehensive sensitivity analysis methodology for nonlinear systems” (abbreviated as “4th-CASAM-N”) since it enables the hitherto very difficult, if not intractable, exact computation of all of the first-, second-, third-, and fourth-order response sensitivities for large-scale systems involving many parameters, as is usually encountered in practice. The foundation for the material presented in this work is provided by the first-order adjoint sensitivity analysis procedure for nonlinear systems that was originally formulated in a general, functional analytic framework by Cacuci [1,2].

The aims and means of sensitivity theory/analysis are occasionally confused with the aims and means of optimization theory. The algorithms underlying optimization theory compute the values at which the first-order derivatives of a response with respect to the state functions and/or parameters vanish. In contradistinction, sensitivity theory/analysis aims at computing the response sensitivities to parameters at the nominal values for the model’s parameters and state functions. Therefore, although both sensitivity analysis and optimization algorithms evaluate first-order (and occasionally higher-order) response derivatives with respect to model parameters, these algorithms serve different purposes and conceptually differ from each other.

Response sensitivities to model parameters, which are computed using the methods of sensitivity analysis, are needed in many activities, including: (i) determining the effects of parameter variations on the system’s behavior; (ii) ranking the importance of model parameters in influencing the system response under consideration; (iii) quantifying uncertainties induced in responses by parameter uncertainties (e.g., by using the method of “propagation of uncertainties”); (iv) prioritizing possible improvements for the system under consideration and possibly reducing conservatism and redundancy; (v) validating the model under consideration by comparison to experiments, while taking into account both experimental and model uncertainties; and (vi) performing “predictive modeling” (which includes data assimilation and model calibration) for the purpose of obtaining best-estimate predicted results with reduced predicted uncertainties, using, e.g., the methodologies presented in [3,4,5] and in references therein.

First-order response sensitivities can be computed by using either statistical or deterministic methods; a comparative review of the most popular of these methods was presented in [6,7]. It is known that sensitivities cannot be computed exactly using statistical methods; this can only be performed with deterministic methods. Furthermore, for a system comprising $TP$ parameters, the computation of the first-order sensitivities by statistical methods requires $O\left(TP\right)$ large-scale computations while the adjoint sensitivity analysis originally developed by Cacuci [1,2] requires a single large-scale computation per response.

Since nonlinear operators do not admit bona-fide adjoint operators (only a linearized form of a nonlinear operator admits an adjoint operator), responses of nonlinear models can depend only on the forward functions. In contradistinction, model responses for linear systems may involve the solutions of both the forward and the adjoint linear models that correspond to the respective physical system. Hence, responses for linear systems cannot always be treated as particular cases of nonlinear systems but there is a need to develop a dedicated sensitivity analysis methodology for response-coupled forward and adjoint linear systems. The general methodology for computing arbitrarily high-order sensitivities for response-coupled linear forward/adjoint systems was developed by Cacuci [8,9,10,11,12]. The overwhelming impact of the higher-order (i.e., second-, third- and fourth-order) sensitivities on the model response was illustrated by Fang and Cacuci, in [13] and references therein, by means of an OECD/NEA reactor physics benchmark modeled by the linear neutron transport equation and comprising 21,976 uncertain model parameters.

The general mathematical framework for adjoint sensitivity analysis of nonlinear systems was extended by Cacuci to second-order [14] and third-order [15] sensitivities. The extension of the adjoint sensitivity analysis methodology to include the computation of sensitivities of model responses with respect to imprecisely known domain boundaries was presented in [16]. By extending all of the previous theoretical developments of the adjoint sensitivity analysis methodology, this work presents the theoretical framework for the fourth-order comprehensive sensitivity analysis methodology for nonlinear systems (abbreviated as “4th-CASAM-N”). This work is structured as follows. Section 2 presents the mathematical formulation of the that would define the computational model of a generic nonlinear physical system, including the definition of a generic response which depends on the model’s state variables and parameters, which are considered to be uncertain (i.e., not known precisely). Besides initial conditions and correlations, the model parameters are also considered to include geometrical parameters that describe the system’s boundaries and internal interfaces. Section 3 presents the 4th-CASAM-N methodology while Section 4 summarizes the salient features of this novel methodology. The application of the principles underlying the 4th-CASAM-N is illustrated in an accompanying work [17] by means of a paradigm nonlinear heat conduction model.

2. Generic Mathematical Modeling of a Nonlinear Physical System

As already mentioned, the computational model of a physical system comprises that relate the system’s independent variables and parameters to the system’s state variables. The model parameters usually stem from processes that are external to the system under consideration and are seldom, if ever, known precisely. The known characteristics of the model parameters may include their nominal (expected/mean) values and, possibly, higher-order moments or cumulants (i.e., variance/covariances, skewness, kurtosis), which are usually determined from experimental data and/or processes external to the physical system under consideration. Occasionally, only inequality and/or equality constraints that delimit the ranges of the system’s parameters are known. Without loss of generality, the imprecisely known model parameters can be considered to be real-valued scalar quantities. These model parameters will be denoted as ${\mathsf{\alpha }}_1$,…,${\mathsf{\alpha }}_{TP}$, where $TP$ denotes the “total number of imprecisely known parameters” underlying the model under consideration. For subsequent developments, it is convenient to consider that these parameters are components of a “vector of parameters” denoted as $\mathsf{\alpha }\triangleq {\left({\mathsf{\alpha }}_1,\dots ,{\mathsf{\alpha }}_{TP}\right)}^{†}\in {\mathrm{E}}_{\mathsf{\alpha }}\in {ℝ}^{TP}$, where ${\mathrm{E}}_{\mathsf{\alpha }}$ is also a normed linear space and where ${ℝ}^{TP}$ denotes the $TP$-dimensional subset of the set of real scalars. The components of the $TP$-dimensional column vector $\mathsf{\alpha }\in {ℝ}^{TP}$ are considered to include imprecisely known geometrical parameters that characterize the physical system’s boundaries in the phase space of the model’s independent variables. Matrices will be denoted using capital bold letters while vectors will be denoted using either capital or lower-case bold letters. The symbol “$\triangleq$” will be used to denote “is defined as” or “is by definition equal to”. Transposition will be indicated by a dagger $\left(†\right)$ superscript.

The generic nonlinear model is considered to comprise $TI$ independent variables which will be denoted as $x_i,i=1,\dots ,TI$, and which are considered to be components of a $TI$-dimensional column vector denoted as $x\triangleq {\left(x_1,\dots ,x_{TI}\right)}^{†}\in {ℝ}^{TI}$, where the sub/superscript “TI” denotes the “total number of independent variables”. The vector $x\in {ℝ}^{TI}$ of independent variables is considered to be defined on a phase-space domain which will be denoted as $\mathsf{\Omega }\left(\alpha \right)$ and which is defined as follows: $\mathsf{\Omega }\left(\alpha \right)\triangleq \left\{-\infty \le {\mathsf{\lambda }}_i\left(\alpha \right)\le x_i\le {\mathsf{\omega }}_i\left(\alpha \right)\le \infty ;i=1,\dots ,TI\right\}$. The lower boundary point of an independent variable is denoted as ${\mathsf{\lambda }}_i\left(\alpha \right)$ and the corresponding upper boundary point is denoted as ${\mathsf{\omega }}_i\left(\alpha \right)$. A typical example of boundaries that depend on both geometrical parameters and material properties are the “boundaries facing vacuum” in models based on diffusion theory, where conditions are imposed on the “extrapolated boundary” of the respective spatial domain. The “extrapolated boundary” depends both on the imprecisely known physical dimensions of the problem’s domain and also on the medium’s properties, such as atomic number densities and microscopic transport cross sections. The boundary of $\mathsf{\Omega }\left(\alpha \right)$, which will be denoted as $\partial \mathsf{\Omega }\left(\alpha \right)$, comprises the set of all of the endpoints ${\mathsf{\lambda }}_i\left(\alpha \right),{\mathsf{\omega }}_i\left(\alpha \right),i=1,\dots ,TI,$ of the respective intervals on which the components of $x$ are defined, i.e., $\partial \mathsf{\Omega }\left(\alpha \right)\triangleq \left\{{\mathsf{\lambda }}_i\left(\alpha \right)\cup {\mathsf{\omega }}_i\left(\alpha \right),i=1,\dots ,TI\right\}$.

A nonlinear physical system can be generally modeled by means of coupled which can be represented in operator form as follows:

$$N\left[u\left(x\right),\alpha \right]=Q\left(x,\alpha \right),x\in {\mathsf{\Omega }}_x\left(\alpha \right).$$

(1)

The quantities which appear in Equation (1) are defined as follows:

1. $u\left(x\right)\triangleq {\left[u_1\left(x\right),\dots ,u_{TD}\left(x\right)\right]}^{†}$ is a $TD$-dimensional column vector of dependent variables; the abbreviation “$TD$” denotes “total number of dependent variables”. The functions $u_i\left(x\right),i=1,\dots ,TD$ denote the system’s “dependent variables” (also called “state functions”); $u\left(x\right)\in {\mathrm{E}}_u$, where ${\mathrm{E}}_u$ is a normed linear space over the scalar field $\mathrm{F}$ of real numbers.

2. $N\left[u\left(x\right);\alpha \right]\triangleq {\left[N_1\left(u;\alpha \right),\dots ,N_{TD}\left(u;\alpha \right)\right]}^{†}$ denotes a $TD$-dimensional column vector. The components $N_i\left(u;\alpha \right),i=1,\dots ,TD$ are operators (including differential, difference, integral, distributions, and/or finite or infinite matrices) acting (usually) nonlinearly on the dependent variables $u\left(x\right)$, the independent variables $x$, and the model parameters $\alpha$. The mapping $N\left(u;\alpha \right)$ is defined on the combined domains of the model’s parameters and state functions, i.e., $\mathbf{N}:\mathrm{D}\subset \mathrm{E}\to {\mathrm{E}}_Q$, where $\mathrm{D}={\mathrm{D}}_u\oplus {\mathrm{D}}_{\mathsf{\alpha }}$, ${\mathrm{D}}_u\subset {\mathrm{E}}_u$, ${\mathrm{D}}_{\mathsf{\alpha }}\subset {\mathrm{E}}_{\mathsf{\alpha }}$, $\mathrm{E}={\mathrm{E}}_u\oplus {\mathrm{E}}_{\mathsf{\alpha }}$.

3. $Q\left(x,\alpha \right)\triangleq {\left[q_1\left(x;\alpha \right),\dots .,q_{TD}\left(x;\alpha \right)\right]}^{†}$ is a $TD$-dimensional column vector which represents inhomogeneous source terms, which usually depend nonlinearly on the uncertain parameters $\alpha$. The vector $Q\left(x,\alpha \right)$ is defined on a normed linear space denoted as ${\mathrm{E}}_Q$, i.e., $Q\in {\mathrm{E}}_Q$.

The equalities in this work are considered to hold in the weak (“distributional”) sense. The right sides of Equation (1) and of other various to be derived in this work may contain “generalized functions/functionals”, particularly Dirac distributions and derivatives thereof.

Boundary and/or initial conditions must also be provided if differential operators appear in Equation (1). In operator form, these boundaries and/or initial conditions are represented as follows:

$$\mathbf{B}\left[u\left(x\right);\alpha ;x\right]-C\left(x,\alpha \right)=0,x\in \partial {\mathsf{\Omega }}_x\left(\alpha \right).$$

(2)

where the column vector $0$ has $TD$ components, all of which are identically zero, i.e.,

$$0\triangleq \left[{\zeta }_1,\dots ,{\zeta }_{TD}\right];{\zeta }_i\equiv 0,i=1,\dots ,TD.$$

(3)

In Equation (2), the components $B_i\left(u;\alpha \right),i=1,\dots ,TD$ of $B\left(u;\alpha \right)\triangleq {\left[B_1\left(u;\alpha \right),\dots ,B_{TD}\left(u;\alpha \right)\right]}^{†}$ are nonlinear operators in $u\left(x\right)$ and $\alpha$, which are defined on the boundary $\partial {\mathsf{\Omega }}_x\left(\alpha \right)$ of the model’s domain ${\mathsf{\Omega }}_x\left(\alpha \right)$. The components $C_i\left(x;\alpha \right),i=1,\dots ,TD$ of $C\left(x,\alpha \right)\triangleq {\left[C_1\left(x;\alpha \right),\dots ,C_{TD}\left(x;\alpha \right)\right]}^{†}$ comprise inhomogeneous boundary sources which are nonlinear functions of $\alpha$.

Solving Equations (1) and (2) at the nominal parameter values, denoted as ${\alpha }^0\triangleq {\left[{\mathsf{\alpha }}_1^0,\dots ,{\mathsf{\alpha }}_i^0,..,{\mathsf{\alpha }}_{TP}^0\right]}^{†}$, provides the “nominal solution” $u^0\left(x\right)$, i.e., the vectors $u^0\left(x\right)$ and ${\alpha }^0$ satisfy the following:

$$N\left[u^0\left(x\right);{\alpha }^0\right]=Q\left(x,{\alpha }^0\right),x\in {\mathsf{\Omega }}_x,$$

(4)

$$B\left[u^0\left(x\right);{\alpha }^0;x\right]-C\left(x,{\alpha }^0\right)=0,x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right).$$

(5)

The superscript “0” will be used throughout this work to denote “nominal values”.

The results computed using a mathematical model are customarily called “model responses” (or “system responses” or “objective functions” or “indices of performance”). In general, a function-valued (i.e., operator-type) response $R\left[u\left(x\right);\alpha \right]$ can be represented by a spectral expansion in multidimensional orthogonal polynomials or Fourier series of the form:

$$R\left[u\left(x\right);\alpha \right]={\displaystyle \sum _{m_1}\dots }{\displaystyle \sum _{m_{TI}}{c}_{m_1\dots m_{TI}}\left[u\left(x\right);\alpha \right]P_{m_1}\left(x_1\right)P_{m_2}\left(x_2\right)}\dots P_{m_{TI}}\left(x_{TI}\right),$$

(6)

where the quantities $P_{m_i}\left(x_i\right)$, $i=1,\dots ,TI$ denote the corresponding spectral functions (e.g., orthogonal polynomials or Fourier exponential/trigonometric functions) and where the spectral Fourier coefficients $c_{m_1\dots m_{TI}}\left[u\left(x\right);\alpha \right]$ are defined as follows:

$$c_{m_1\dots m_{TI}}\left[u\left(x\right);\alpha \right]\triangleq {\displaystyle \underset{{\mathsf{\lambda }}_1\left(\alpha \right)}{\overset{{\mathsf{\omega }}_1\left(\alpha \right)}{\int }}\dots {\displaystyle \underset{{\mathsf{\lambda }}_i\left(\alpha \right)}{\overset{{\mathsf{\omega }}_i\left(\alpha \right)}{\int }}\dots {\displaystyle \underset{{\mathsf{\lambda }}_{TI}\left(\alpha \right)}{\overset{{\mathsf{\omega }}_{TI}\left(\alpha \right)}{\int }}R\left[u\left(x\right);\alpha \right]P_{m_1}\left(x_1\right)\dots P_{m_i}\left(x_i\right)\dots P_{m_{TI}}\left(x_{TI}\right)d{x}_1\dots d{x}_i\dots d{x}_{TI}}}}$$

(7)

The coefficients $c_{m_1\dots m_{TI}}\left[u\left(x\right);\alpha \right]$ can themselves be considered as “model responses” since the spectral polynomials $P_{m_i}\left(x_i\right)$ are perfectly well known while the expansion coefficients will contain all of the dependencies (directly or indirectly—through the state functions) of the respective response on the imprecisely known model parameters. This way, the sensitivity analysis of an operator-valued response $R\left[u\left(x\right);\alpha \right]$ can be reduced to the sensitivity analysis of the scalar-valued responses $c_{m_1\dots m_{TI}}\left[u\left(x\right);\alpha \right]$.

A measurement of a physical quantity that depends on the model’s state functions and parameters can be considered to be a response denoted as $R_p\left[u\left(x\right);\alpha \right]$, which is to be evaluated at $x=x^p\left(\alpha \right)$, where $x^p\left(\alpha \right)\triangleq {\left[x_1^p\left(\alpha \right),\dots x_k^p\left(\alpha \right),\dots ,x_{TI}^p\left(\alpha \right)\right]}^{†}$ denotes the location in the phase space of the specific “measurement point”. Such a measurement (or measurement-like) response can be represented mathematically as follows:

$$R_p\left[u\left(x\right);\alpha \right]\triangleq {\displaystyle \underset{{\mathsf{\lambda }}_1\left(\alpha \right)}{\overset{{\mathsf{\omega }}_1\left(\alpha \right)}{\int }}\dots {\displaystyle \underset{{\mathsf{\lambda }}_{TI}\left(\alpha \right)}{\overset{{\mathsf{\omega }}_{TI}\left(\alpha \right)}{\int }}F\left[u\left(x\right);\alpha ;x\right]\delta \left[x_1-x_1^p\left(\alpha \right)\right]\dots \delta \left[x_{TI}-x_{TI}^p\left(\alpha \right)\right]d{x}_1\dots d{x}_{TI}}}$$

(8)

where the function $F\left[u\left(x\right);\alpha ;x\right]$ denotes the mathematical dependence of the measurement device on the model’s dependent variable(s), and where the quantity $\delta \left[x_i-x_i^p\left(\alpha \right)\right]$ denotes the Dirac delta functional. The measurement’s location in phase space, $x^p\left(\alpha \right)$, may itself be afflicted by measurement (experimental) uncertainties. Hence, it is convenient to consider the components of $x^p\left(\alpha \right)$ to be included among the components of the vector $\alpha$ of model parameters, even though $x^p\left(\alpha \right)$ appears only in the definition of the response but does not appear in Equations (1) and (2), which mathematically define the physical model. Thus, the physical “system” is defined to comprise both the system’s computational model and the system’s response. In most cases, the coordinates $x_k^p\left(\alpha \right)$, $k=1,\dots ,TI$ will be independent (albeit uncertain) model parameters, in which case $\partial x_k^p\left(\alpha \right)/\partial {\mathsf{\alpha }}_n=1,if{\mathsf{\alpha }}_n\equiv x_k^p$ and $\partial x_k^p\left(\alpha \right)/\partial {\mathsf{\alpha }}_n=0,if{\mathsf{\alpha }}_n\ne x_k^p$.

The representations shown in Equations (6)–(8) indicate that model responses can be fundamentally analyzed by considering the following generic integral representation:

$$R\left[u\left(x\right);\alpha \right]\triangleq {\displaystyle \underset{{\mathsf{\lambda }}_1\left(\alpha \right)}{\overset{{\mathsf{\omega }}_1\left(\alpha \right)}{\int }}\dots {\displaystyle \underset{{\mathsf{\lambda }}_{TI}\left(\alpha \right)}{\overset{{\mathsf{\omega }}_{TI}\left(\alpha \right)}{\int }}S\left[u\left(x\right);\alpha ;x\right]d{x}_1\dots d{x}_{TI}}}$$

(9)

where $S\left[u\left(x\right);\alpha \right]$ is suitably differentiable nonlinear function of $u\left(x\right)$ and of $\alpha$. It is important to note that the components of $\alpha$ not only include parameters that appear in the defining the computational model per se, i.e., in Equations (1) and (2), but also include parameters that specifically occur only in the definition of the response under consideration.

It is also important to note that the system’s definition domain, $\mathsf{\Omega }\left(\alpha \right)$, in phase space is considered to be imprecisely known, subject to uncertainties in the components of the vector of model parameters $\alpha$. Therefore, the system domain’s boundary, $\partial \mathsf{\Omega }\left(\alpha \right)$, as well as the model response $R\left[u\left(x\right);\alpha \right]$, will be affected by the boundary uncertainties that affect the endpoints ${\mathsf{\lambda }}_i\left(\alpha \right),{\mathsf{\omega }}_i\left(\alpha \right),i=1,\dots ,TI$. Such boundary uncertainties stem most often from manufacturing uncertainties.

3. The Fourth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (4th-CASAM-N)

The starting point for building the mathematical framework of the novel 4th-CASAM-N is provided by the mathematical framework of the “first-order comprehensive adjoint sensitivity analysis methodology for nonlinear systems” presented in [16] which considered that the domain’s boundaries are also subject to uncertainties, thus generalizing all previous work on first-order adjoint sensitivity analysis. The 1st-CASAM-N is briefly reviewed in Section 3.1. Section 3.2 presents the 2nd-CASAM-N, which generalizes the material presented in [9,14,15]. Section 3.3 presents the 3rd-CASAM-N, comprising original material that provides the basis for the development of the 4th-CASAM, which is presented in Section 3.4.

3.1. The First-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (1st-CASAM-N)

The model and boundary parameters $\alpha$ are considered to be uncertain quantities, having unknown true values. The nominal (or mean) parameter vales ${\alpha }^0$ are considered to be known, and these will differ from the true values by quantities denoted as $\mathsf{\delta }\alpha \triangleq \left(\mathsf{\delta }{\mathsf{\alpha }}_1,\dots ,\mathsf{\delta }{\mathsf{\alpha }}_{TP}\right)$, where $\mathsf{\delta }{\mathsf{\alpha }}_i\triangleq {\mathsf{\alpha }}_i-{\mathsf{\alpha }}_i^0$. Since the forward state functions $u\left(x\right)$ are related to the model and boundary parameters $\alpha$ through Equations (1) and (2), it follows that the variations $\mathsf{\delta }\alpha$ in the model and boundary parameters will cause corresponding variations $v^{\left(1\right)}\left(x\right)\triangleq {\left[\mathsf{\delta }{u}_1\left(x\right),\dots ,\mathsf{\delta }{u}_{TD}\left(x\right)\right]}^{†}$ around the nominal solution $u^0\left(x\right)$ in the forward state functions. In turn, the variations $\mathsf{\delta }\alpha$ and $v^{\left(1\right)}\left(x\right)$ will induce variations in the system’s response. Cacuci [1,2] has shown that the most general definition of the sensitivity of an operator-valued model response $R\left(e\right)$, where $e\triangleq \left(\alpha ,u\right)\in \mathrm{E}$, to variations $h\triangleq \left(\mathsf{\delta }\alpha ,v^{\left(1\right)}\right)$ in the model parameters and state functions in a neighborhood around the nominal functions and parameter values $e^0\triangleq \left({\alpha }^0,u^0\right)\in \mathrm{E}$, is given by the first-order Gateaux (G) variation, which will be denoted as $\mathsf{\delta }R\left(e^0;h\right)$ and is defined as follows:

$$\mathsf{\delta }R\left(e^0;h\right)\triangleq {\left\{\frac{d}{d\mathsf{\epsilon }}\left[R\left(e^0+\mathsf{\epsilon }h\right)\right]\right\}}_{\mathsf{\epsilon }=0}\triangleq \underset{\mathsf{\epsilon }\to 0}{\lim }\frac{R\left(e^0+\mathsf{\epsilon }h\right)-R\left(e^0\right)}{\mathsf{\epsilon }}$$

(10)

for a scalar $\mathsf{\epsilon }\in \mathrm{F}$ and for all (i.e., arbitrary) vectors $h\in \mathrm{E}={\mathrm{E}}_{\mathsf{\alpha }}\times {\mathrm{E}}_u$ in a neighborhood $\left(e^0+\mathsf{\epsilon }h\right)$ around $e^0=\left({\alpha }^0,u^0\right)\in \mathrm{E}$. The G variation $\mathsf{\delta }R\left(e^0;h\right)$ is an operator defined on the same domain as $R\left(e\right)$ and has the same range as $R\left(e\right)$. The G variation $\mathsf{\delta }R\left(e^0;h\right)$ satisfies the relation: $R\left(e^0+\mathsf{\epsilon }h\right)-R\left(e^0\right)=\mathsf{\delta }R\left(e^0;h\right)+\mathsf{\Delta }\left(h\right)$, with $\underset{\mathsf{\epsilon }\to 0}{\lim }\left[\mathsf{\Delta }\left(\mathsf{\epsilon }h\right)\right]/\mathsf{\epsilon }=0$. The existence of the G variation $\mathsf{\delta }R\left(e^0;h\right)$ does not guarantee its numerical computability. Numerical methods most often require that $\mathsf{\delta }R\left(e^0;h\right)$ is linear in $h\triangleq \left(\mathsf{\delta }\alpha ;v^{\left(1\right)}\right)$ in a neighborhood $\left(e^0+\mathsf{\epsilon }h\right)$ around $e^0=\left(u^0,{\alpha }^0\right)\in \mathrm{E}$. Formally, the necessary and sufficient conditions for the G variation $\mathsf{\delta }R\left(e^0;h\right)$ of a nonlinear operator $R\left(e\right)$ to be linear and continuous in $h$ in a neighborhood $\left(e^0+\mathsf{\epsilon }h\right)$ around $e^0=\left({\alpha }^0,u^0\right)$ and therefore admit a total first-order G derivative, are as follows:

(i)

$R\left(e\right)$ satisfies a weak Lipschitz condition at $e^0$, namely:

$$‖R\left(e^0+\mathsf{\epsilon }h\right)-R\left(e^0\right)‖\le k‖\mathsf{\epsilon }\left(e^0\right)‖,k<\infty$$

(11)

(ii)

$R\left(e\right)$ satisfies the following condition

$$R\left(e^0+\mathsf{\epsilon }{h}_1+\mathsf{\epsilon }{h}_2\right)-R\left(e^0+\mathsf{\epsilon }{h}_1\right)-R\left(e^0+\mathsf{\epsilon }{h}_2\right)+R\left(e^0\right)=o\left(\mathsf{\epsilon }\right);h_1,h_2\in \mathrm{E};\mathsf{\epsilon }\in \mathrm{F}.$$

(12)

In practice, the relation provided in Equations (11) and (12) are seldom used directly since the computation of the expression on the right side of Equation (10) reveals if the respective expression is linear (or not) in $h\triangleq \left(\mathsf{\delta }\alpha ,v^{\left(1\right)}\right)$ and, hence, in $v^{\left(1\right)}\left(x\right)$. Numerical methods (e.g., Newton’s method, and variants thereof) for solving Equations (1) and (2) also require the existence of the first-order G derivatives of original model, in which case the components of the operators which appear in these must also satisfy the conditions provided in Equations (11) and (12). Therefore, the conditions provided in Equations (11) and (12) are henceforth considered to be satisfied by the operators which underly the physical system modeled by Equations (1) and (2), as well as by the model responses.

When the first-order G differential $\mathsf{\delta }R\left(e^0;h\right)$ satisfies the conditions provided in Equations (11) and (12), it can be written in the following form:

$$\mathsf{\delta }R\left(e^0;h\right)\equiv {\left\{\mathsf{\delta }R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right);\mathsf{\delta }\alpha \right]\right\}}_{{\alpha }^0}={\left\{\frac{\partial R\left(u;\alpha \right)}{\partial u}\right\}}_{{\alpha }^0}{v}^{\left(1\right)}\left(x\right)+{\left\{\frac{\partial R\left(u;\alpha \right)}{\partial \alpha }\right\}}_{{\alpha }^0}\mathsf{\delta }\alpha ,$$

(13)

where

$$\frac{\partial \left[\right]}{\partial u}{v}^{\left(1\right)}\left(x\right)\triangleq {\displaystyle \sum _{i=1}^{TD}\frac{\partial \left[\right]}{\partial u_i\left(x\right)}\mathsf{\delta }{u}_i\left(x\right)}$$

(14)

$$\frac{\partial \left[\right]}{\partial \alpha }\mathsf{\delta }\alpha \triangleq {\displaystyle \sum _{i=1}^{TP}\frac{\partial \left[\right]}{\partial {\mathsf{\alpha }}_i}\mathsf{\delta }{\mathsf{\alpha }}_i}$$

(15)

This, the quantities $\partial R\left(u;\alpha \right)/\partial u$ and $\partial R\left(u;\alpha \right)/\partial \alpha$ in Equation (13) denote the partial G derivatives of $R\left(e\right)$ with respect to $u$ and $\alpha$, evaluated at the nominal parameter values (and hence also nominal values of the state functions). The notation ${\left\{\right\}}_{{\alpha }^0}$ will be used in this work to indicate that the quantity enclosed within the bracket is to be evaluated at the respective nominal parameter and state functions values. The quantity ${\left\{\partial R\left(u;\alpha \right)/\partial \alpha \right\}}_{{\alpha }^0}\mathsf{\delta }\alpha$ is called the “direct-effect term” because it arises directly from parameter variations $\mathsf{\delta }\alpha$. The direct-effect term can be computed once the nominal values $e^0=\left(u^0,{\alpha }^0\right)$ are available. The quantity ${\left\{\partial R\left(u;\alpha \right)/\partial u\right\}}_{{\alpha }^0}{v}^{\left(1\right)}\left(x\right)$ is called the “indirect-effect term” because it arises indirectly, through the variations in the state functions (which are caused through the model by parameter variations). The indirect-effect term can be quantified only after having determined the variations $v^{\left(1\right)}\left(x\right)$ in terms of the variations $\mathsf{\delta }\alpha$.

In particular, the first-order Gateaux differential ${\left\{\mathsf{\delta }R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right);\mathsf{\delta }\alpha \right]\right\}}_{{\alpha }^0}$ of the generic response $R\left(u,\alpha \right)$ defined in Equation (9) has the following expression:

$$\begin{array}{l}{\left\{\mathsf{\delta }R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right);\mathsf{\delta }\alpha \right]\right\}}_{{\alpha }^0}\\ \triangleq \left\{\frac{d}{d\mathsf{\epsilon }}{\displaystyle \underset{{\mathsf{\lambda }}_1\left({\alpha }^0+\mathsf{\epsilon }\mathsf{\delta }\alpha \right)}{\overset{{\mathsf{\omega }}_1\left({\alpha }^0+\mathsf{\epsilon }\mathsf{\delta }\alpha \right)}{\int }}\dots {\displaystyle \underset{{\mathsf{\lambda }}_{TI}\left({\alpha }^0+\mathsf{\epsilon }\mathsf{\delta }\alpha \right)}{\overset{{\mathsf{\omega }}_{TI}\left({\alpha }^0+\mathsf{\epsilon }\mathsf{\delta }\alpha \right)}{\int }}S\left(u^0+\mathsf{\epsilon }{v}^{\left(1\right)},{\alpha }^0+\mathsf{\epsilon }\mathsf{\delta }\alpha ;x\right)d{x}_1\dots d{x}_{TI}}}\right\}\\ ={\left\{\mathsf{\delta }R\left[u\left(x\right);\alpha ;\mathsf{\delta }\alpha \right]\right\}}_{dir}+{\left\{\mathsf{\delta }R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind},\end{array}$$

(16)

where the “indirect-effect” term ${\left\{\mathsf{\delta }R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}$ comprises all dependencies on $v^{\left(1\right)}\left(x\right)$ and is defined as follows:

$${\left\{\mathsf{\delta }R\left[u\left(x\right);\alpha ;v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}\triangleq {\displaystyle \underset{{\mathsf{\lambda }}_1\left({\alpha }^0\right)}{\overset{{\mathsf{\omega }}_1\left({\alpha }^0\right)}{\int }}\dots {\displaystyle \underset{{\mathsf{\lambda }}_{TI}\left({\alpha }^0\right)}{\overset{{\mathsf{\omega }}_{TI}\left({\alpha }^0\right)}{\int }}{\left\{\frac{\partial S\left(u;\alpha ;x\right)}{\partial u}\right\}}_{{\alpha }^0}{v}^{\left(1\right)}\left(x\right)d{x}_1\dots d{x}_{TI}}},$$

(17)

and where the “direct-effect” term ${\left\{\mathsf{\delta }R\left[u\left(x\right);\alpha ;\mathsf{\delta }\alpha \right]\right\}}_{dir}$ comprises all dependencies on $\mathsf{\delta }\alpha$ and is defined as follows:

$${\left\{\mathsf{\delta }R\left[u\left(x\right);\alpha ;\mathsf{\delta }\alpha \right]\right\}}_{dir}\triangleq {\displaystyle \sum _{j_1=1}^{TP}{\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);\alpha \right]\right\}}_{dir}\mathsf{\delta }{\mathsf{\alpha }}_{j_1}},$$

(18)

with

$$\begin{array}{l}{\left\{R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);\mathsf{\alpha }\right]\right\}}_{dir}\triangleq {\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\frac{\partial S\left(\mathbf{u};\mathsf{\alpha };\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\mathrm{\dots }d{x}_{TI}\right\}}_{{\mathsf{\alpha }}^0}\\ +{\displaystyle {\displaystyle \sum }_{j=1}^{TI}}{\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\mathrm{\dots }\underset{{\mathsf{\lambda }}_{j-1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j-1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\underset{{\mathsf{\lambda }}_{j+1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j+1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\mathrm{\dots }\underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}S\left[\mathbf{u}\left(x_1,.,{\mathsf{\omega }}_j\left(\mathsf{\alpha }\right),.,x_{N_x}\right);\mathsf{\alpha }\right]\frac{\partial {\mathsf{\omega }}_j\left(\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\mathrm{\dots }d{x}_{TI}\right\}}_{{\mathsf{\alpha }}^0}\\ -{\displaystyle {\displaystyle \sum }_{j=1}^{TI}}{\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\mathrm{\dots }\underset{{\mathsf{\lambda }}_{j-1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j-1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\underset{{\mathsf{\lambda }}_{j+1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j+1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\mathrm{\dots }\underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}S\left[\mathbf{u}\left(x_1,.,{\mathsf{\lambda }}_j\left(\mathsf{\alpha }\right),.,x_{N_x}\right);\mathsf{\alpha }\right]\frac{\partial {\mathsf{\lambda }}_j\left(\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\mathrm{\dots }d{x}_{TI}\right\}}_{{\mathsf{\alpha }}^0}.\end{array}$$

(19)

The first-order relationship between the vectors $v^{\left(1\right)}\left(x\right)$ and $\mathsf{\delta }\mathsf{\alpha }$ is determined by taking the G differentials of Equations (1) and (2). Thus, applying the definition of the G differential to Equations (1) and (2) yields the following:

$${\left\{\frac{d}{d\mathsf{\epsilon }}N\left[u^0+\mathsf{\epsilon }{v}^{\left(1\right)}\left(x\right);{\mathsf{\alpha }}^0+\mathsf{\epsilon }\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{\mathsf{\epsilon }=0}={\left\{\frac{d}{d\mathsf{\epsilon }}Q\left(x,{\mathsf{\alpha }}^0+\mathsf{\epsilon }\mathsf{\delta }\mathsf{\alpha }\right)\right\}}_{\mathsf{\epsilon }=0},$$

(20)

$${\left\{\frac{d}{d\mathsf{\epsilon }}B\left[u^0+\mathsf{\epsilon }{v}^{\left(1\right)}\left(x\right);{\mathsf{\alpha }}^0+\mathsf{\epsilon }\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{\mathsf{\epsilon }=0}-{\left\{\frac{d}{d\mathsf{\epsilon }}C\left(x,{\mathsf{\alpha }}^0+\mathsf{\epsilon }\mathsf{\delta }\mathsf{\alpha }\right)\right\}}_{\mathsf{\epsilon }=0}=0.$$

(21)

Carrying out the differentiations with respect to $\mathsf{\epsilon }$ in Equations (20) and (21), and setting $\mathsf{\epsilon }=0$ in the resulting expressions yields the following:

$${\left\{V^{\left(1\right)}\left(u;\mathsf{\alpha }\right)v^{\left(1\right)}\left(x\right)\right\}}_{{\mathsf{\alpha }}^0}={\left\{q^{\left(1\right)}\left(u;\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\right\}}_{{\mathsf{\alpha }}^0},x\in {\mathsf{\Omega }}_x,$$

(22)

$${\left\{b_V^{\left(1\right)}\left(u;\mathsf{\alpha };v^{\left(1\right)};\mathsf{\delta }\mathsf{\alpha }\right)\right\}}_{{\mathsf{\alpha }}^0}=0,x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right).$$

(23)

In Equations (22) and (23), the superscript “(1)” indicates “1st-Level” and the various quantities which appear in these are defined as follows:

$$V^{\left(1\right)}\left(u;\mathsf{\alpha }\right)\triangleq \left\{\frac{\partial N\left(u;\mathsf{\alpha }\right)}{\partial u}\right\}\triangleq \left(\begin{array}{ccc}\frac{\partial N_1}{\partial u_1}& \dots & \frac{\partial N_1}{\partial u_{TD}}\\ ⋮& \ddots & ⋮\\ \frac{\partial N_{TD}}{\partial u_1}& \cdots & \frac{\partial N_{TD}}{\partial u_{TD}}\end{array}\right);$$

(24)

$$q^{\left(1\right)}\left(u;\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\triangleq \frac{\partial \left[Q\left(\mathsf{\alpha }\right)-N\left(u;\mathsf{\alpha }\right)\right]}{\partial \mathsf{\alpha }}\mathsf{\delta }\mathsf{\alpha }\triangleq {\displaystyle \sum _{j_1=1}^{TP}{s}_V^{\left(1\right)}\left(j_1;u;\mathsf{\alpha }\right)\mathsf{\delta }{\mathsf{\alpha }}_{j_1}};$$

(25)

$$s_V^{\left(1\right)}\left(j_1;u;\mathsf{\alpha }\right)\triangleq \frac{\partial \left[Q\left(\mathsf{\alpha }\right)-N\left(u;\mathsf{\alpha }\right)\right]}{\partial {\mathsf{\alpha }}_{j_1}};$$

(26)

$${\left\{b_V^{\left(1\right)}\left(u;\mathsf{\alpha };v^{\left(1\right)};\mathsf{\delta }\mathsf{\alpha }\right)\right\}}_{{\mathsf{\alpha }}^0}\triangleq {\left\{\frac{\partial B\left(u;\mathsf{\alpha }\right)}{\partial u}{v}^{\left(1\right)}\right\}}_{{\mathsf{\alpha }}^0}+{\displaystyle \sum _{j_1=1}^{TP}{b}_V^{\left(1\right)}\left(j_1;u;\mathsf{\alpha }\right)\mathsf{\delta }{\mathsf{\alpha }}_{j_1}};$$

(27)

$$b_V^{\left(1\right)}\left(j_1;u;\mathsf{\alpha }\right)\triangleq \frac{\partial \left[B\left(u;\mathsf{\alpha }\right)-C\left(\mathsf{\alpha }\right)\right]}{\partial {\mathsf{\alpha }}_{j_1}}.$$

(28)

The system comprising Equations (22) and (23) is called the “1st-Level Variational Sensitivity System” (1st-LVSS). In order to determine the solutions of the 1st-LVSS that would correspond to every parameter variation $\mathsf{\delta }{\mathsf{\alpha }}_{j_1}$, $j_1=1,\dots ,TP$, the 1st-LVSS would need to be solved $TP$ times, with distinct right sides for each $\mathsf{\delta }{\mathsf{\alpha }}_{j_1}$, as follows:

$${\left\{V^{\left(1\right)}\left(u;\mathsf{\alpha }\right)v^{\left(1\right)}\left(j_1;x\right)\right\}}_{{\mathsf{\alpha }}^0}={\left\{s_V^{\left(1\right)}\left(j_1;u;\mathsf{\alpha }\right)\right\}}_{{\mathsf{\alpha }}^0};j_1=1,\dots ,TP;x\in {\mathsf{\Omega }}_x,$$

(29)

$${\left\{\frac{\partial B\left(u;\mathsf{\alpha }\right)}{\partial u}{v}^{\left(1\right)}\left(j_1;x\right)\right\}}_{{\mathsf{\alpha }}^0}=b_V^{\left(1\right)}\left(j_1;\mathsf{\alpha };u\right);j_1=1,\dots ,TP;x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right).$$

(30)

Subsequently, the solutions $v^{\left(1\right)}\left(j_1;x\right)$ could be used, in turn, in Equation (17) to compute the indirect-effect term corresponding to each parameter variation $\mathsf{\delta }{\mathsf{\alpha }}_{j_1}$, $j_1=1,\dots ,TP$, to obtain the following contribution from the indirect-effect term to the respective partial sensitivity of the response with respect to the parameter ${\mathsf{\alpha }}_{j_1}$:

$${\left\{\mathsf{\delta }R\left[j_1;u\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(j_1;x\right)\right]\right\}}_{ind}\triangleq {\left\{{\displaystyle \underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{\int }}\dots {\displaystyle \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{\int }}\left[\frac{\partial S\left(u;\mathsf{\alpha };x\right)}{\partial u}{v}^{\left(1\right)}\left(j_1;x\right)\right]d{x}_1\dots d{x}_{TI}}}\right\}}_{{\mathsf{\alpha }}^0}.$$

(31)

Adding the contribution from the indirect-effect term obtained in Equation (31) to the contribution from the direct-effect term provided in Equation (18) yields the following expression for the sensitivity (i.e., partial G derivative) $\partial R\left[u\left(x\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_1}$ of the response $R\left[u\left(x\right);\mathsf{\alpha }\right]$ to the parameter ${\mathsf{\alpha }}_{j_1}$, $j_1=1,\dots ,TP$:

$$\frac{\partial R\left[u\left(x\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_1}}={\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);\mathsf{\alpha }\right]\right\}}_{dir}+{\left\{\mathsf{\delta }R\left[j_1;u\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(j_1;x\right)\right]\right\}}_{ind}.$$

(32)

The quantities $\partial R\left[u\left(x\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_1}$ are independent of the parameter variations $\mathsf{\delta }{\mathsf{\alpha }}_{i_1}$ and represent the first-order partial sensitivities (first-order partial G derivatives) of the response $R\left(e\right)$ with respect to each of the model parameters ${\mathsf{\alpha }}_{j_1}$, $j_1=1,\dots ,TP$, evaluated at the nominal values $e^0=\left({\mathsf{\alpha }}^0,u^0\right)$. Computing the response sensitivities by using the solutions $v^{\left(1\right)}\left(j_1;x\right)$, $j_1=1,\dots ,TP$, of the 1st-LVSS requires $TP$ large-scale forward computations in order to determine the functions $v^{\left(1\right)}\left(j_1;x\right)$, $j_1=1,\dots ,TP$. Since most problems of practical interest are characterized by many parameters (i.e., $\mathsf{\alpha }$ has many components) and comparatively few responses, it becomes prohibitively expensive to solve repeatedly the 1st-LVSS in order to determine the functions $v^{\left(1\right)}\left(j_1;x\right)$, $j_1=1,\dots ,TP$. Even though the 1st-LVSS contains first-order parameter and state-functions variations, it is called “first-level” (rather than “first-order”) in anticipation of determining second-order sensitivities, which will use “second-level” forward and adjoint systems. These “second-level” systems will not be called “second-order” because they will not contain second-order parameter and/or state-function variations, but will also contain only first-order variations, even though they will be used for determining second-order sensitivities. Similar terminology, i.e., “third-level” (as opposed to “third-order”) forward/adjoint systems, will be used for determining the third-order sensitivities, and so on.

In most practical situations, the number of model parameters significantly exceeds the number of functional responses of interest, i.e., $TR\ll TP$, so it would be advantageous to perform just $TR$ (rather than $TP$) computations. The goal of the “1st-order comprehensive adjoint sensitivity analysis methodology for nonlinear systems (1st-CASAM-N)” is to compute exactly and efficiently the “indirect effect term” defined in Equation (17) without needing to compute explicitly the vectors $v^{\left(1\right)}\left(j_1;x\right)$, $j_1=1,\dots ,TP$. The qualifier “comprehensive” is meant to indicate that that the 1st-CASAM-N considers that the internal and external boundaries $\partial {\mathsf{\Omega }}_x\left(\mathsf{\alpha }\right)$ of the phase-space domain depend on the uncertain model parameters $\mathsf{\alpha }$ and are thereby imprecisely known, subject to uncertainties. Thus, the 1st-CASAM-N represents a generalization of the pioneering works by Cacuci [1,2] that conceived the “adjoint sensitivity analysis methodology”, in which the domain boundary was considered to be perfectly well known, free of uncertainties. The fundamental ideas underlying the 1st-CASAM-N are as conceived by Cacuci [1,2], aiming at eliminating the appearance of the vectors $v^{\left(1\right)}\left(j_1;x\right)$ from the expression of the indirect-effect term defined in Equation (17). This elimination is achieved by expressing the right side of Equation (17) in terms of the solutions of the “1st-Level Adjoint Sensitivity System (1st-LASS)”, the construction of which requires the introduction of adjoint operators. Adjoint operators can be defined in Banach spaces but are most useful in Hilbert spaces. Since real Hilbert spaces provide the natural mathematical setting for computational purposes, the derivations presented in this section are set in real (as opposed to complex) Hilbert spaces, without affecting the generality of the concepts presented herein. Thus, the spaces ${\mathrm{E}}_u$ and ${\mathrm{E}}_Q$ are henceforth considered to be self-dual Hilbert spaces and will be denoted as ${\mathrm{H}}_1\left({\mathsf{\Omega }}_x\right)$. The inner product of two vectors $u^{(a)}\left(x\right)\in {\mathrm{H}}_1$ and $u^{(b)}\left(x\right)\in {\mathrm{H}}_1$ will be denoted as ${〈u^{(a)},u^{(b)}〉}_1$, and is defined as follows:

$${〈u^{(a)},u^{(b)}〉}_1\triangleq {\left\{{\displaystyle \underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{\int }}\dots {\displaystyle \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{\int }}\left[u^{(a)}\left(x\right)·u^{(b)}\left(x\right)\right]d{x}_1\dots d{x}_{TI}}}\right\}}_{{\mathsf{\alpha }}^0}$$

(33)

where the dot indicates the “scalar product of two vectors” defined as follows: $u^{(a)}\left(x\right)·u^{(b)}\left(x\right)\triangleq {\displaystyle \sum _{i=1}^{TD}{u}_i^{(a)}\left(x\right)u_i^{(b)}}\left(x\right)$. It is important to note that the inner product defined in Equation (33) is continuous in $\mathsf{\alpha }$, i.e., it holds at any value particular value of $\mathsf{\alpha }$, including at the nominal parameter values ${\mathsf{\alpha }}^0$.

The construction of the 1st-LASS commences by noting that the vector $v^{\left(1\right)}\left(x\right)$ itself is independent of the index $j_1=1,\dots ,TP$, since

$$v^{\left(1\right)}\left(x\right)={\displaystyle \sum _{j_1=1}^{TP}{v}^{\left(1\right)}\left(j_1;x\right)}.$$

(34)

The next step is to form the inner product of Equation (22) with a vector $a^{(1)}\left(x\right)\triangleq \left[a_1^{\left(1\right)}\left(x\right),\dots ,a_{TD}^{\left(1\right)}\left(x\right)\right]\in {\mathrm{H}}_1$, where the superscript “(1)” indicates “1st-Level”:

$${\left\{{〈a^{(1)}\left(x\right),V^{\left(1\right)}\left(u;\mathsf{\alpha }\right)v^{\left(1\right)}\left(x\right)〉}_1\right\}}_{{\mathsf{\alpha }}^0}={\left\{{〈a^{(1)}\left(x\right),q^{\left(1\right)}\left(u;\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)〉}_1\right\}}_{{\mathsf{\alpha }}^0}.$$

(35)

Using the definition of the adjoint operator in ${\mathrm{H}}_1\left({\mathsf{\Omega }}_x\right)$, the left side of Equation (35) is transformed as follows:

$${\left\{{〈a^{(1)},V^{\left(1\right)}\left(u;\mathsf{\alpha }\right)v^{\left(1\right)}〉}_1\right\}}_{{\mathsf{\alpha }}^0}={\left\{{〈A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)a^{(1)},v^{\left(1\right)}〉}_1\right\}}_{{\mathsf{\alpha }}^0}+{\left\{{\left[P^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};v^{\left(1\right)}\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0},$$

(36)

where ${\left[P^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};v^{\left(1\right)}\right)\right]}_{\partial {\mathsf{\Omega }}_x}$ denotes the associated bilinear concomitant evaluated on the space/time domain’s boundary $\partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right)$ and where $A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)$ is the operator adjoint to $V^{\left(1\right)}\left(u;\mathsf{\alpha }\right)$, i.e., $A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)\triangleq {\left[V^{\left(1\right)}\left(u;\mathsf{\alpha }\right)\right]}^{\ast }$. The symbol ${\left[\right]}^{\ast }$ will be used in this work to indicate “adjoint” operator. In certain situations, it might be computationally advantageous to include certain boundary components of ${\left[P^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};v^{\left(1\right)}\right)\right]}_{\partial {\mathsf{\Omega }}_x}$ into the components of $A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)$.

The domain of $A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)$ is determined by selecting appropriate adjoint boundary and/or initial conditions, which will be denoted in operator form as:

$${\left\{b_A^{\left(1\right)}\left(u;a^{\left(1\right)};\mathsf{\alpha }\right)\right\}}_{{\mathsf{\alpha }}^0}=0,x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right).$$

(37)

The above boundary conditions for $A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)$ are usually inhomogeneous, i.e., $b_A^{\left(1\right)}\left(0;0;\mathsf{\alpha }\right)\ne 0$, and are obtained by imposing the following requirements:

1. They must be independent of unknown values of $v^{\left(1\right)}\left(x\right)$ and $\mathsf{\delta }\mathsf{\alpha }$;

2. The substitution of the boundary and/or initial conditions represented by Equations (23) and (37) into the expression of ${\left\{{\left[P^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};v^{\left(1\right)}\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ must cause all terms containing unknown values of $v^{\left(1\right)}\left(x\right)$ to vanish.

Constructing the adjoint initial and/or boundary conditions for $A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)$ as described above and implementing them together with the variational boundary and/or initial conditions represented by Equation (23) into Equation (35) reduces the bilinear concomitant ${\left\{{\left[P^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};v^{\left(1\right)}\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ to a quantity that will contain boundary terms involving only known values of $\mathsf{\delta }\mathsf{\alpha }$, ${\mathsf{\alpha }}^0$, $u^0$, and ${\mathsf{\psi }}^{(1)}$; this quantity will be denoted by ${\left\{{\left[{\widehat{P}}^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$. In general, the boundary terms represented by ${\left\{{\left[{\widehat{P}}^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ do not vanish automatically. In certain cases, however, ${\left\{{\left[{\widehat{P}}^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ may vanish automatically or it may be forced to vanish by considering appropriately constructed extensions of the adjoint operator $A^{\left(1\right)}\left(\mathsf{\alpha },u\right)$; however, such extensions are seldom needed in practice. Since ${\left\{{\left[{\widehat{P}}^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ is linear in $\mathsf{\delta }\mathsf{\alpha }$, it can be expressed in the following form: ${\left[{\widehat{P}}^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}={\displaystyle \sum _{j_1=1}^{TP}\left[\partial {\widehat{P}}^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)}\right)/\partial {\mathsf{\alpha }}_{j_1}\right]\mathsf{\delta }{\mathsf{\alpha }}_{j_1}}.$

Implementing the forward and adjoint boundary and/or initial conditions given in Equations (23) and (37) into Equation (36) transforms the later into the following relation:

$${\left\{{〈A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)a^{(1)},v^{\left(1\right)}〉}_1\right\}}_{{\mathsf{\alpha }}^0}={\left\{{〈a^{(1)},V^{\left(1\right)}\left(u;\mathsf{\alpha }\right)v^{\left(1\right)}〉}_1\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0},$$

(38)

Replacing the quantity $V^{\left(1\right)}\left(\mathsf{\alpha };u\right)v^{\left(1\right)}$ in the first term on the right side of Equation (38) by the right side of Equation (22) yields the following relation:

$${\left\{{〈A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)a^{(1)},v^{\left(1\right)}〉}_1\right\}}_{{\mathsf{\alpha }}^0}={\left\{{〈a^{(1)},q^{\left(1\right)}\left(u;\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)〉}_1\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}.$$

(39)

The definition of the function $a^{(1)}\left(x\right)$ will now be completed by requiring that the left side of Equation (39) is the same as the indirect-effect term defined in Equation (17), which is achieved by imposing the following relationship:

$${\left\{A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)a^{\left(1\right)}\left(x\right)\right\}}_{{\mathsf{\alpha }}^0}={\left\{\partial S\left(u;\mathsf{\alpha }\right)/\partial u\right\}}_{{\mathsf{\alpha }}^0}\triangleq s_A^{\left(1\right)}\left[u\left(x\right);\mathsf{\alpha }\right],x\in {\mathsf{\Omega }}_x,$$

(40)

while satisfying the adjoint boundary conditions represented by Equation (37). The subscript “A” attached to the source term on the right side of Equation (40) indicates “adjoint”. Since the source $s_A^{\left(1\right)}\left[u\left(x\right);\mathsf{\alpha }\right]$ may contain distributions (e.g., Dirac delta functions and derivatives thereof), the equality in Equation (40) is considered to hold in the weak sense. The well-known Riesz representation theorem ensures that the relationship in Equation (40) holds uniquely.

The results obtained in Equations (38)–(40) are now replaced in Equation (17) to obtain the following expression of the indirect-effect term as a function of $a^{(1)}\left(x\right)$:

$$\begin{array}{l}{\left\{\mathsf{\delta }R\left[u\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}={\left\{{〈a^{(1)},q^{\left(1\right)}\left(u;\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)〉}_1\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}\\ \equiv {\left\{\mathsf{\delta }R\left[u\left(x\right);\mathsf{\alpha };a^{(1)}\left(x\right);\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{ind}\triangleq {\displaystyle \sum _{j_1=1}^{TP}{\left\{R^{\left(1\right)}\left[j_1;u\left(x\right);a^{(1)}\left(x\right);\mathsf{\alpha }\right]\right\}}_{ind}\mathsf{\delta }{\mathsf{\alpha }}_{j_1}},\end{array}$$

(41)

where, for each $j_1=1,\dots ,TP$, the contribution of the indirect-effect term to the sensitivity of the response with respect to the parameter ${\mathsf{\alpha }}_{j_1}$ is given by

$$\begin{array}{l}{\left\{R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]\right\}}_{ind}\triangleq {\left\{\langle {\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{s}}_V^{\left(1\right)}{\left(j_1;\mathbf{u};\mathsf{\alpha }\right)\rangle }_1-\partial {\widehat{P}}^{\left(1\right)}\left(\mathbf{u};\mathsf{\alpha };{\mathbf{a}}^{\left(1\right)}\right)/\partial {\mathsf{\alpha }}_{j_1}\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\mathrm{\dots }\underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)·\frac{\partial \left[\mathbf{Q}\left(\mathsf{\alpha }\right)-\mathbf{N}\left(\mathbf{u};\mathsf{\alpha }\right)\right]}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\mathrm{\dots }d{x}_{TI}-\frac{\partial {\widehat{P}}^{\left(1\right)}\left(\mathbf{u};\mathsf{\alpha };{\mathbf{a}}^{\left(1\right)}\right)}{\partial {\mathsf{\alpha }}_{j_1}}\right\}}_{{\mathsf{\alpha }}^0}.\end{array}$$

(42)

As the identity on the right side of Equation (41) indicates, the desired elimination of all unknown values of $v^{\left(1\right)}\left(x\right)$ from the expression of the indirect-effect term ${\left\{\mathsf{\delta }R\left[u\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}$ is accomplished. Instead of depending on $v^{\left(1\right)}\left(x\right)$, the indirect-effect term ${\left\{\mathsf{\delta }R\left[u\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}$ now depends on the adjoint function $a^{(1)}\left(x\right)\in {\mathrm{H}}_1$.

Replacing in Equation (16) the result obtained in Equation (41) together with the expression provided (18) for the direct-effect term yields the following expression

$$\begin{array}{l}{\left\{\mathsf{\delta }R\left[u\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right);\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}={\left\{\mathsf{\delta }R\left[u\left(x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}+{\left\{{〈a^{(1)},q^{\left(1\right)}\left(u;\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)〉}_1\right\}}_{{\mathsf{\alpha }}^0}\\ -{\left\{{\left[{\widehat{P}}^{\left(1\right)}\left(u;\mathsf{\alpha };a^{(1)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}.\end{array}$$

(43)

The expressions of the first-order response sensitivities ${\left\{\partial R\left[u\left(x\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_1}\right\}}_{{\mathsf{\alpha }}^0}$ of the response $R\left[u\left(x\right);\mathsf{\alpha }\right]$ with respect to the parameters ${\mathsf{\alpha }}_{j_1}$ are obtained identifying in Equation (43) the quantities that multiply the respective parameter variations $\mathsf{\delta }{\mathsf{\alpha }}_{j_1}$, $j_1=1,\dots ,TP$. This identification yields the following expressions for the first-order response sensitivities ${\left\{\partial R\left[u\left(x\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_1}\right\}}_{{\mathsf{\alpha }}^0}$, $j_1=1,\dots ,TP$, computed at the model’s nominal parameter and state function values:

$$\begin{array}{l}{\left\{\partial R\left[\mathbf{u}\left(\mathbf{x}\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_1}\right\}}_{{\mathsf{\alpha }}^0}\triangleq {\left\{R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]\right\}}_{ind}+{\left\{R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);\mathbf{\alpha }\right]\right\}}_{dir}\\ ={\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\mathrm{\dots }\underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)·\frac{\partial \left[\mathbf{Q}\left(\mathsf{\alpha }\right)-\mathbf{N}\left(\mathbf{u};\mathsf{\alpha }\right)\right]}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\mathrm{\dots }d{x}_{TI}-\frac{\partial {\widehat{P}}^{\left(1\right)}\left(\mathbf{u};\mathsf{\alpha };{\mathbf{a}}^{\left(1\right)}\right)}{\partial {\mathsf{\alpha }}_{j_1}}\right\}}_{{\mathsf{\alpha }}^0}\\ +{\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\frac{\partial S\left(\mathbf{u};\mathsf{\alpha };\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI}\right\}}_{{\mathsf{\alpha }}^0}\\ +{\displaystyle {\displaystyle \sum }_{j=1}^{TI}}{\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{j-1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j-1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\underset{{\mathsf{\lambda }}_{j+1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j+1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}S\left[\mathbf{u}\left(x_1,.,{\mathsf{\omega }}_j\left(\mathsf{\alpha }\right),.,x_{N_x}\right);\mathsf{\alpha }\right]\frac{\partial {\mathsf{\omega }}_j\left(\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI}\right\}}_{{\mathsf{\alpha }}^0}\\ -{\displaystyle {\displaystyle \sum }_{j=1}^{TI}}{\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{j-1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j-1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\underset{{\mathsf{\lambda }}_{j+1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j+1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}S\left[\mathbf{u}\left(x_1,.,{\mathsf{\lambda }}_j\left(\mathsf{\alpha }\right),.,x_{N_x}\right);\mathsf{\alpha }\right]\frac{\partial {\mathsf{\lambda }}_j\left(\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI}\right\}}_{{\mathsf{\alpha }}^0}.\\ \triangleq R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right];j_1=1,\dots ,TP.\end{array}$$

(44)

As indicated by Equation (44), each of the first-order sensitivities $R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha }\right]$ of the response $R\left[u\left(x\right);\mathsf{\alpha }\right]$ with respect to the model parameters ${\mathsf{\alpha }}_{j_1}$ (including boundary and initial conditions) can be computed inexpensively after having obtained the function $a^{(1)}\left(x\right)\in {\mathrm{H}}_1$, using just quadrature formulas to evaluate the various inner products involving $a^{(1)}\left(x\right)\in {\mathrm{H}}_1$ in the expression of the indirect-effect term obtained in Equation (42). The function $a^{(1)}\left(x\right)\in {\mathrm{H}}_1$ is obtained by solving numerically Equations (37) and (40), which is the only large-scale computation needed for obtaining all of the first-order sensitivities. Equations (37) and (40) will be called the first-level adjoint sensitivity system (1st-LASS), and its solution, $a^{(1)}\left(x\right)\in {\mathrm{H}}_1\left({\mathsf{\Omega }}_x\right)$, will be called the first-level adjoint function. It is very important to note that the 1st-LASS is independent of parameter variation $\mathsf{\delta }{\mathsf{\alpha }}_{j_1}$, $j_1=1,\dots ,TP$, and therefore needs to be solved only once, regardless of the number of model parameters under consideration. Furthermore, since Equation (40) is linear in $a^{(1)}\left(x\right)$${\mathsf{\psi }}_{1,i_1}^{\left(2\right)}\left(x\right)$, solving it requires less computational effort than solving the original Equation (1), which is nonlinear in $u\left(x\right)$.

3.2. The Second-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (2nd-CASAM-N)

The 2nd-CASAM-N relies on the same fundamental concepts as introduced in [14], but the 2nd-CASAM-N also enables the computation of response sensitivities with respect to imprecisely known domain boundaries, thus including all possible types of uncertain parameters. Fundamentally, the second-order sensitivities are defined as the “1st-order sensitivities of the 1st-order sensitivities”. This definition stems from the inductive definition of the second-order total G differential of correspondingly differentiable function, which is also defined inductively as “the total 1st-order differential of the 1st-order total differential” of a function.

The computation of the second-order sensitivities $R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha }\right]$ requires the prior computation of the original forward function $u\left(x\right)$ and the first-level adjoint function $a^{\left(1\right)}\left(x\right)$. For each $j_1=1,\dots ,TP$, the first-order sensitivities $R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha }\right]$ will be assumed to satisfy the conditions stated in Equations (11) and (12). Under these conditions, the first-order total G differential ${\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right);\mathsf{\delta }{a}^{\left(1\right)}\left(x\right);\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}$ of $R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha }\right]$ will exist in a neighborhood around the nominal values of the parameters and state functions and will have, by definition, the following expression, for each $j_1=1,\dots ,TP$:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha };{\mathbf{v}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\delta }{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}\\ \triangleq {\left\{\frac{d}{d\mathsf{\epsilon }}\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right)+\mathsf{\epsilon }{\mathbf{v}}^{\left(1\right)}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)+\mathsf{\epsilon }\mathsf{\delta }{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }+\mathsf{\epsilon }\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{\mathsf{\epsilon }=0}\\ \triangleq {\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}+{\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha };{\mathbf{v}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\delta }{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)\right]\right\}}_{ind},\end{array}$$

(45)

where

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}\triangleq {\left\{\frac{\partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\mathsf{\delta }\mathsf{\alpha }\right\}}_{{\mathsf{\alpha }}^0}\\ ={\displaystyle \sum _{j_2=1}^{TP}\frac{\partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_2}}}\mathsf{\delta }{\mathsf{\alpha }}_{j_2},\end{array}$$

(46)

and

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha };{\mathbf{v}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\delta }{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)\right]\right\}}_{ind}\\ \triangleq {\displaystyle {\displaystyle \sum }_{i=1}^{TD}}\frac{\partial R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial u_i\left(\mathbf{x}\right)}\mathbf{\delta }{u}_i\left(\mathbf{x}\right)+{\displaystyle {\displaystyle \sum }_{i=1}^{TD}}\frac{\partial R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial a_i^{\left(1\right)}\left(\mathbf{x}\right)}\mathsf{\delta }{a}_i^{\left(1\right)}\left(\mathbf{x}\right)\\ \triangleq {\left\{\partial R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]/\partial u\right\}}_{{\mathsf{\alpha }}^0}{\mathbf{v}}^{\left(1\right)}\left(\mathbf{x}\right)\\ +{\left\{\partial R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]/\partial {\mathbf{a}}^{\left(1\right)}\right\}}_{{\mathsf{\alpha }}^0}\mathsf{\delta }{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right).\end{array}$$

(47)

In Equation (45), the quantities ${\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right);\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}$ and ${\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}$ denote, respectively, the indirect-effect term and the direct-effect term. The direct-effect term ${\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}$ comprises all dependencies on the vector $\mathsf{\delta }\mathsf{\alpha }$ of parameter variations and, in view of Equation (44), has the following expression obtained by using Equation (44):

$$\begin{array}{l}\frac{\partial R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_2}}=\frac{\partial }{\partial {\mathsf{\alpha }}_{j_2}}\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\frac{\partial S\left(\mathbf{u};\mathsf{\alpha };\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\mathrm{\dots }d{x}_{TI}\right\}\\ +\frac{\partial }{\partial {\mathsf{\alpha }}_{j_2}}\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)·\frac{\partial \left[\mathbf{Q}\left(\mathsf{\alpha }\right)-\mathbf{N}\left(\mathbf{u};\mathsf{\alpha }\right)\right]}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI}-\frac{\partial {\widehat{P}}^{\left(1\right)}\left(\mathbf{u};\mathsf{\alpha };{\mathbf{a}}^{\left(1\right)}\right)}{\partial {\mathsf{\alpha }}_{j_1}}\right\}\\ +{\displaystyle {\displaystyle \sum }_{j=1}^{TI}}\frac{\partial }{\partial {\mathsf{\alpha }}_{j_2}}\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{j-1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j-1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\underset{{\mathsf{\lambda }}_{j+1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j+1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}S\left[\mathbf{u}\left(x_1,.,{\mathsf{\omega }}_j\left(\mathsf{\alpha }\right),.,x_{N_x}\right);\mathsf{\alpha }\right]\frac{\partial {\mathsf{\omega }}_j\left(\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI}\right\}\\ -{\displaystyle {\displaystyle \sum }_{j=1}^{TI}}\frac{\partial }{\partial {\mathsf{\alpha }}_{j_2}}\left\{\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{j-1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j-1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\underset{{\mathsf{\lambda }}_{j+1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j+1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}S\left[\mathbf{u}\left(x_1,.,{\mathsf{\lambda }}_j\left(\mathsf{\alpha }\right),.,x_{N_x}\right);\mathsf{\alpha }\right]\frac{\partial {\mathsf{\lambda }}_j\left(\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI}\right\},\\ for{j}_1,j_2=1,\dots ,TP.\end{array}$$

(48)

The indirect-effect term ${\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right);\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}$ comprises all dependencies on the vectors $v^{\left(1\right)}\left(x\right)$ and $\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)$ of variations in the state functions $u\left(x\right)$ and $a^{\left(1\right)}\left(x\right)$, respectively. The components of the indirect-effect term have the following expressions obtained by using Equation (44):

$$\begin{array}{l}\frac{\partial R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial u_i}=\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\frac{{\partial }^{2}S\left(\mathbf{u};\mathsf{\alpha };\mathsf{\alpha }\right)}{\partial u_i\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI}\\ +\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots {{\displaystyle \int }}_{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}^{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)·\frac{\partial \left[\mathbf{Q}\left(\mathsf{\alpha }\right)-\partial \mathbf{N}\left(\mathbf{u};\mathsf{\alpha }\right)/\partial u_i\right]}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI}-\frac{{\partial }^2{\widehat{P}}^{\left(1\right)}\left(\mathbf{u};\mathsf{\alpha };{\mathbf{a}}^{\left(1\right)}\right)}{\partial u_i\partial {\mathsf{\alpha }}_{j_1}}\\ +{\displaystyle {\displaystyle \sum }_{j=1}^{TI}}\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{j-1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j-1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\underset{{\mathsf{\lambda }}_{j+1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j+1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\frac{\partial S\left[\mathbf{u}\left(x_1,.,{\mathsf{\omega }}_j\left(\mathsf{\alpha }\right),.,x_{N_x}\right);\mathsf{\alpha }\right]}{\partial u_i}\frac{\partial {\mathsf{\omega }}_j\left(\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI}\\ -{\displaystyle {\displaystyle \sum }_{j=1}^{TI}}\underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\mathrm{\dots }\underset{{\mathsf{\lambda }}_{j-1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j-1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\underset{{\mathsf{\lambda }}_{j+1}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{j+1}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\dots \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{{\displaystyle \int }}}\frac{\partial S\left[\mathbf{u}\left(x_1,.,{\mathsf{\lambda }}_j\left(\mathsf{\alpha }\right),.,x_{N_x}\right);\mathsf{\alpha }\right]}{\partial u_i}\frac{\partial {\mathsf{\lambda }}_j\left(\mathsf{\alpha }\right)}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI},\\ for{j}_1=1,\dots ,TP;i=1,\dots ,TD;\end{array}$$

(49)

and

$$\begin{array}{l}\frac{\partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha }\right]}{\partial a_i^{\left(1\right)}}=\frac{\partial }{\partial a_i^{\left(1\right)}}\left\{{\displaystyle \underset{{\mathsf{\lambda }}_1\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_1\left(\mathsf{\alpha }\right)}{\int }}\dots {\displaystyle \underset{{\mathsf{\lambda }}_{TI}\left(\mathsf{\alpha }\right)}{\overset{{\mathsf{\omega }}_{TI}\left(\mathsf{\alpha }\right)}{\int }}{a}^{(1)}\left(x\right)·\frac{\partial \left[Q\left(\mathsf{\alpha }\right)-N\left(u;\mathsf{\alpha }\right)\right]}{\partial {\mathsf{\alpha }}_{j_1}}d{x}_1\dots d{x}_{TI}}}\right\},\\ for{j}_1=1,\dots ,TP;i=1,\dots ,TD.\end{array}$$

(50)

The direct-effect term ${\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}$ can be computed immediately while the indirect-effect term ${\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right);\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}$ can be computed only after having determined the vectors $v^{\left(1\right)}\left(x\right)$ and $\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)$. The vector $v^{\left(1\right)}\left(x\right)$ is the solution of the 1st-LVSS defined by Equations (22) and (23). On the other hand, the vector $\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)$ is the solution of the G differentiated 1st-LASS. Thus, taking the G differential of Equations (37) and (40) yields the following system of equations for $\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)$:

$${\left\{V_{21}^{\left(2\right)}\left(u;a^{\left(1\right)};\mathsf{\alpha }\right)v^{\left(1\right)}\left(x\right)+V_{22}^{\left(2\right)}\left(u;\mathsf{\alpha }\right)\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)\right\}}_{{\mathsf{\alpha }}^0}={\left\{q_2^{\left(2\right)}\left(u;a^{\left(1\right)};\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\right\}}_{{\mathsf{\alpha }}^0},x\in {\mathsf{\Omega }}_x,$$

(51)

$${\left\{\mathsf{\delta }{b}_A^{\left(1\right)}\left(u;a^{\left(1\right)};\mathsf{\alpha };v^{\left(1\right)};\mathsf{\delta }{a}^{\left(1\right)};\mathsf{\delta }\mathsf{\alpha }\right)\right\}}_{{\mathsf{\alpha }}^0}=0,x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right),$$

(52)

where

$$V_{21}^{\left(2\right)}\left(u;a^{\left(1\right)};\mathsf{\alpha }\right)\triangleq \frac{\partial \left[A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)a^{\left(1\right)}\right]}{\partial u}-\frac{s^{\left(1\right)}\left[u\left(x\right);\mathsf{\alpha }\right]}{\partial u};$$

(53)

$$V_{22}^{\left(2\right)}\left(u;\mathsf{\alpha }\right)\triangleq A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)$$

(54)

$$q_2^{\left(2\right)}\left(u;\mathsf{\alpha };a^{\left(1\right)};\mathsf{\delta }\mathsf{\alpha }\right)\triangleq \frac{s^{\left(1\right)}\left[u\left(x\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\mathsf{\delta }\mathsf{\alpha }-\frac{\partial \left[A^{\left(1\right)}\left(u;\mathsf{\alpha }\right)a^{\left(1\right)}\left(x\right)\right]}{\partial \mathsf{\alpha }}\mathsf{\delta }\mathsf{\alpha };$$

(55)

$$\mathsf{\delta }{b}_A^{\left(1\right)}\left(u;a^{\left(1\right)};\mathsf{\alpha }\right)\triangleq \frac{\partial b_A^{\left(1\right)}\left(u;a^{\left(1\right)};\mathsf{\alpha }\right)}{\partial u}{v}^{\left(1\right)}\left(x\right)+\frac{\partial b_A^{\left(1\right)}\left(u;a^{\left(1\right)};\mathsf{\alpha }\right)}{\partial a^{\left(1\right)}}\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)+\frac{\partial b_A^{\left(1\right)}\left(u;a^{\left(1\right)};\mathsf{\alpha }\right)}{\partial \mathsf{\alpha }}\mathsf{\delta }\mathsf{\alpha }$$

(56)

Although Equations (51) and (52) are coupled to the 1st-LVSS, they can be solved sequentially, after having obtained the solution $v^{\left(1\right)}\left(x\right)$ of the 1st-LVSS. Formally, the functions $v^{\left(1\right)}\left(x\right)$ and $\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)$ are obtained by solving the following second-level variational sensitivity system (2nd-LVSS), which is obtained by concatenating the 1st-LVSS to Equations (51) and (52):

$${\left\{V{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]V^{\left(2\right)}\left(2;x\right)\right\}}_{{\mathsf{\alpha }}^0}={\left\{Q_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0},x\in {\mathsf{\Omega }}_x,$$

(57)

$${\left\{B_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[2\right],0\left[2\right]\triangleq {\left[0,0\right]}^{†},x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right)$$

(58)

The argument “2” which appears in the list of arguments of the vector $U^{\left(2\right)}\left(2;x\right)$ and the “variational vector” $V^{\left(2\right)}\left(2;x\right)$ in Equations (57) and (58) indicates that each of these vectors is a two-block column vector, with each block comprising a column vector of dimension $TD$, defined as follows:

$$U^{\left(2\right)}\left(2;x\right)\triangleq \left(\begin{array}{c}{u}^{\left(1\right)}\left(x\right)\\ a^{\left(1\right)}\left(x\right)\end{array}\right);V^{\left(2\right)}\left(2;x\right)\triangleq \mathsf{\delta }{U}^{\left(2\right)}\left(2;x\right)\triangleq \left(\begin{array}{c}{v}^{\left(2\right)}\left(1;x\right)\\ v^{\left(2\right)}\left(2;x\right)\end{array}\right)\triangleq \left(\begin{array}{c}{v}^{\left(1\right)}\left(x\right)\\ \mathsf{\delta }{a}^{\left(1\right)}\left(x\right)\end{array}\right).$$

(59)

Thus, the (column) block vector $U^{\left(2\right)}\left(2;x\right)$ has a total of $2\times TD$ components; evidently, the (column) block vector $V^{\left(2\right)}\left(2;x\right)$ also has a total of $2\times TD$ components. In the relatively simple case regarding the components of either the vector $U^{\left(2\right)}\left(2;x\right)$ or the vector $V^{\left(2\right)}\left(2;x\right)$, the numbers “1” and “2” could also be used as subscripts, but such a subscript notation would become unwieldy for the higher-level (adjoint) functions, which will be introduced in the sections to follow below. The superscript “(2)” which appears in the notation of the vectors $U^{\left(2\right)}\left(2;x\right)$ and $V^{\left(2\right)}\left(2;x\right)$ indicates “2nd-level”. Henceforth, such “higher-level” (i.e., level higher than first) variational and adjoint functions/vectors will be denoted using bold capital letters. The argument “2” in the expression $0\left[2\right]\triangleq {\left[0,0\right]}^{†}$ indicates that the quantity $0\left[2\right]$ is a two-block column vector comprising two vectors, each of which has $TD$ components, all of which are zero-valued, as defined in Equation (3). Thus, the column vector $0\left[2\right]$ has a total of $2\times TD$ components, all of which are identically zero.

To distinguish block vectors from block matrices, two capital bold letter are used (and will henceforth be used) to denote block matrices, as in the case of “the second-level” “variational matrix” $V{M}^{\left(2\right)}\left[2\times 2;u^{\left(2\right)}\left(x\right);\mathsf{\alpha }\right]$. The “2nd-level” is indicated by the superscript “(2)”. Subsequently in this work, levels higher than second will also be indicated by a corresponding superscript attached to the appropriate block vectors and/or block matrices. The argument “$2\times 2$”, which appears in the list of arguments of $V{M}^{\left(2\right)}\left[2\times 2;u^{\left(2\right)}\left(x\right);\mathsf{\alpha }\right]$, indicates that this matrix is a $2\times 2$-dimensional block matrix comprising four matrices, each with of dimensions $TD\times TD$, having the following structure:

$$V{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]\triangleq \left(\begin{array}{cc}{V}^{\left(1\right)}& 0\\ V_{21}^{\left(2\right)}& V_{22}^{\left(2\right)}\end{array}\right).$$

(60)

Thus, the matrix $V{M}^{\left(2\right)}\left[2\times 2;u^{\left(2\right)}\left(x\right);\mathsf{\alpha }\right]$ has a total of ${\left(2\times TD\right)}^2$ components (or elements). The other quantities which appear in Equations (57) and (58) are also two-block vectors, with the same structure as $V^{\left(2\right)}\left(2;x\right)$, and are defined as follows:

$$Q_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\triangleq \left(\begin{array}{c}{q}_V^{\left(2\right)}\left(1;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\\ q_V^{\left(2\right)}\left(2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\end{array}\right)\triangleq \left(\begin{array}{c}{q}^{\left(1\right)}\left(u;\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\\ q_2^{\left(2\right)}\left(u;a^{\left(1\right)};\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\end{array}\right);$$

(61)

$$\begin{array}{l}{B}_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\triangleq \left(\begin{array}{c}{b}_V^{\left(2\right)}\left[1;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\\ b_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\end{array}\right)\\ \triangleq \left(\begin{array}{c}{b}_V^{\left(1\right)}\left(u^{\left(1\right)};\mathsf{\alpha };\mathsf{\delta }{u}^{\left(1\right)};\mathsf{\delta }\mathsf{\alpha }\right)\\ \mathsf{\delta }{b}_A^{\left(1\right)}\left[U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\end{array}\right).\end{array}$$

(62)

Solving the 2nd-LVSS requires $T{P}^2$ large-scale computations, which is unrealistic to perform for large-scale systems comprising many parameters. The 2nd-CASAM-N circumvents the need for solving the 2nd-LVSS by deriving an alternative expression for the indirect-effect term defined in Equation (47), in which the function $V^{\left(2\right)}\left(2;x\right)$ is replaced by a second-level adjoint function which is independent of variations in the model parameter and state functions. This second-level adjoint function will satisfy a second-level adjoint sensitivity system (2nd-LASS), which will be constructed by using the 2nd-LVSS as the starting point, following the same principles outlined in Section 3.1 The 2nd-LASS will be constructed in a Hilbert space which will be denoted as ${\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$ and which comprises as elements block vectors of the same form as $V^{\left(2\right)}\left(2;x\right)$. Thus, a generic vector in ${\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$, denoted as ${\Psi }^{(2)}\left(2;x\right)\triangleq {\left[{\mathsf{\psi }}^{(2)}\left(1;x\right),{\mathsf{\psi }}^{(2)}\left(2;x\right)\right]}^{†}\in {\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$, comprises two components of the form ${\mathsf{\psi }}^{\left(2\right)}\left(1;x\right)\triangleq {\left[{\mathsf{\psi }}_1^{\left(2\right)}\left(1;x\right),\dots ,{\mathsf{\psi }}_{TD}^{\left(2\right)}\left(1;x\right)\right]}^{†}\in {\mathrm{H}}_1\left({\mathsf{\Omega }}_x\right)$ and ${\mathsf{\psi }}^{\left(2\right)}\left(2;x\right)\triangleq {\left[{\mathsf{\psi }}_1^{\left(2\right)}\left(2;x\right),\dots ,{\mathsf{\psi }}_{TD}^{\left(2\right)}\left(2;x\right)\right]}^{†}\in {\mathrm{H}}_1\left({\mathsf{\Omega }}_x\right)$, each of which are $TD$-dimensional column vectors; hence, ${\Psi }^{(2)}\left(2;x\right)$ is a $\left(2\times TD\right)$-dimensional column vector.

The inner product of two vectors ${\Psi }^{(2)}\left(2;x\right)\triangleq {\left[{\mathsf{\psi }}^{(2)}\left(1;x\right),{\mathsf{\psi }}^{(2)}\left(2;x\right)\right]}^{†}\in {\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$ and ${\Phi }^{(2)}\left(x\right)\triangleq {\left[{\mathsf{\phi }}^{(2)}\left(1;x\right),{\mathsf{\phi }}^{(2)}\left(2;x\right)\right]}^{†}\in {\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$ in the Hilbert space ${\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$ will be denoted as ${〈{\Psi }^{(2)}\left(2;x\right),{\Phi }^{(2)}\left(x\right)〉}_2$ and defined as follows:

$${〈{\Psi }^{(2)}\left(2;x\right),{\Phi }^{(2)}\left(x\right)〉}_2\triangleq {\displaystyle \sum _{i=1}^2{〈{\mathsf{\psi }}^{(2)}\left(i;x\right),{\mathsf{\phi }}^{(2)}\left(i;x\right)〉}_1}$$

(63)

The inner product defined in Equation (63) is continuous in $\mathsf{\alpha }$ in a neighborhood of ${\mathsf{\alpha }}^0$. Using the definition of the inner product defined in Equation (63), construct the inner product of Equation (57) with a vector $A^{(2)}\left(2;x\right)\triangleq {\left[a^{(2)}\left(1;x\right),a^{(2)}\left(2;x\right)\right]}^{†}\in {\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$ to obtain the following relation:

$$\begin{array}{l}{\left\{{〈A^{(2)}\left(2;x\right),V{M}^{\left(2\right)}\left[2\times 2;u^{\left(2\right)}\left(x\right);\mathsf{\alpha }\right]V^{\left(2\right)}\left(2;x\right)〉}_2\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{{〈A^{(2)}\left(2;x\right),Q_V^{\left(2\right)}\left[2;u^{\left(2\right)}\left(x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]〉}_2\right\}}_{{\mathsf{\alpha }}^0}.\end{array}$$

(64)

The inner product on the left side of Equation (64) is now further transformed by using the definition of the adjoint operator to obtain the following relation:

$$\begin{array}{l}{\left\{{〈A^{(2)}\left(2;x\right),V{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]V^{\left(2\right)}\left(2;x\right)〉}_2\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{{〈V^{\left(2\right)}\left(2;x\right),A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]A^{(2)}\left(2;x\right)〉}_2\right\}}_{{\mathsf{\alpha }}^0}\\ +{\left\{{\left[P^{\left(2\right)}\left(U^{\left(2\right)};A^{(2)};V^{\left(2\right)};\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0},\end{array}$$

(65)

where the adjoint matrix-valued operator $A{M}^{\left(2\right)}\left[2\times 2;u^{\left(2\right)}\left(x\right);\mathsf{\alpha }\right]$ is defined as follows:

$$A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]\triangleq {\left[V{M}^{\left(2\right)}\left(2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right)\right]}^{\ast }=\left(\begin{array}{cc}{\left[V^{\left(1\right)\ast }\right]}^{†}& {\left[V_{21}^{\left(2\right)\ast }\right]}^{†}\\ 0& {\left[V_{22}^{\left(2\right)\ast }\right]}^{†}\end{array}\right)$$

(66)

The matrix $A{M}^{\left(2\right)}\left[2\times 2;u^{\left(2\right)}\left(x\right);\mathsf{\alpha }\right]$ comprises $\left(2\times 2\right)$ block matrices, each with dimensions $T{D}^2$, thus comprising a total of $\left(2\times 2\right)T{D}^2$ components (or elements).

In Equation (66), the quantity ${\left\{{\left[P^{\left(2\right)}\left(U^{\left(2\right)};A^{(2)};V^{\left(2\right)};\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ denotes the corresponding bilinear concomitant on the domain’s boundary, evaluated at the nominal values for the parameters and respective state functions. The definition domain of the adjoint (matrix-valued) operator $A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]$ is specified by requiring the function $A^{(2)}\left(2;x\right)\triangleq {\left[a^{(2)}\left(1;x\right),a^{(2)}\left(2;x\right)\right]}^{†}\in {\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$ to satisfy adjoint boundary/initial conditions denoted as follows:

$${\left\{B_A^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[2\right],x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right).$$

(67)

The second-level adjoint boundary/initial conditions represented by Equation (67) are determined by requiring that: (a) they must be independent of unknown values of $V^{\left(2\right)}\left(2;x\right)$; and (b) the substitution of the boundary and/or initial conditions represented by Equations (58) and (67) into the expression of ${\left\{{\left[P^{\left(2\right)}\left(U^{\left(2\right)};A^{(2)};V^{\left(2\right)};\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ must cause all terms containing unknown values of $V^{\left(2\right)}\left(2;x\right)$ to vanish.

Implementing the second-level (forward and adjoint) boundary/initial conditions, namely Equations (58) and (67) into Equation (65), will transform the later into the following form:

$$\begin{array}{l}{\left\{{〈V^{\left(2\right)}\left(2;x\right),A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]A^{(2)}\left(2;x\right)〉}_2\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{{〈A^{(2)}\left(2;x\right),V{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]V^{\left(2\right)}\left(2;x\right)〉}_2\right\}}_{{\mathsf{\alpha }}^0}\\ -{\left\{{\left[{\widehat{P}}^{\left(2\right)}\left(U^{\left(2\right)};A^{(2)};\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0},\end{array}$$

(68)

where ${\left\{{\left[{\widehat{P}}^{\left(2\right)}\left(U^{\left(2\right)};A^{(2)};\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ denotes residual boundary terms which do not depend on $V^{\left(2\right)}\left(2;x\right)$ but may not have vanished automatically. The right side of Equation (64) is now used to replace the vector $V{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]V^{\left(2\right)}\left(2;x\right)$ in the first term on the right side of Equation (68), thereby obtaining the following relation:

$$\begin{array}{l}{\left\{{〈V^{\left(2\right)}\left(2;x\right),A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]A^{(2)}\left(2;x\right)〉}_2\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{{〈A^{(2)}\left(2;x\right),Q_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]〉}_2\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(2\right)}\left(U^{\left(2\right)};A^{(2)};\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}.\end{array}$$

(69)

The definition of the second-level adjoint function $A^{(2)}\left(2;x\right)\triangleq {\left[a^{(2)}\left(1;x\right),a^{(2)}\left(2;x\right)\right]}^{†}\in {\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$ is now completed by requiring that both the left side of Equation (69) and the right side of Equation (47) represent the “indirect-effect term” ${\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right)\right]\right\}}_{ind}$. As shown in Equation (45), there are a total of $TP$ indirect-effect terms, each one corresponding to a first-order sensitivity $R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha }\right]$, $j_1=1,\dots ,TP$. Hence, there will be a total of $TP$ second-level adjoint functions of the form $A^{(2)}\left(2;j_1;x\right)\triangleq {\left[a^{(2)}\left(1;j_1;x\right),a^{(2)}\left(2;j_1;x\right)\right]}^{†}\in {\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$, $j_1=1,\dots ,TP$, with each such adjoint function corresponding to a specific $j_1$-dependent indirect-effect term. The left side of Equation (69) will be identical to the right side of Equation (47) by requiring that the following relation is satisfied for each $j_1=1,\dots ,TP$ by the second-level adjoint functions (block vectors) $A^{(2)}\left(2;j_1;x\right)\triangleq {\left[a^{(2)}\left(1;j_1;x\right),a^{(2)}\left(2;j_1;x\right)\right]}^{†}\in {\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$:

$$\begin{array}{l}{\left\{A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]A^{(2)}\left(2;j_1;x\right)\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{Q_A^{\left(2\right)}\left[2;j_1;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0},j_1=1,\dots ,TP,x\in {\mathsf{\Omega }}_x,\end{array}$$

(70)

where

$$\begin{array}{l}{Q}_A^{\left(2\right)}\left[2;j_1;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]\triangleq \left(\begin{array}{c}{q}_A^{\left(2\right)}\left(1;j_1;U^{\left(2\right)};\mathsf{\alpha }\right)\\ q_A^{\left(2\right)}\left(2;j_1;U^{\left(2\right)};\mathsf{\alpha }\right)\end{array}\right)\\ \triangleq \left(\begin{array}{c}\partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right)\right]/\partial u\\ \partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha };v^{\left(1\right)}\left(x\right)\right]/\partial a^{\left(1\right)}\end{array}\right),j_1=1,\dots ,TP.\end{array}$$

(71)

The boundary conditions to be satisfied by each of the second-level adjoint functions $A^{(2)}\left(2;j_1;x\right)\triangleq {\left[a^{(2)}\left(1;j_1;x\right),a^{(2)}\left(2;j_1;x\right)\right]}^{†}\in {\mathrm{H}}_2\left({\mathsf{\Omega }}_x\right)$ are as represented by Equation (67), namely:

$${\left\{B_A^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[2\right];j_1=1,\dots ,TP;x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right).$$

(72)

Since the source $Q_A^{\left(2\right)}\left[2;j_1;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]$ may contain distributions (e.g., Dirac delta functions and derivatives thereof), the equality in Equation (70) is considered to hold in the weak sense. The well-known Riesz representation theorem ensures that the relationship in Equation (70) holds uniquely.

The system of equations represented by Equations (70)–(72) will be called the second-level adjoint sensitivity system (2nd-LASS) and its solution, $A^{(2)}\left(2;j_1;x\right)\triangleq {\left[a^{(2)}\left(1;j_1;x\right),a^{(2)}\left(2;j_1;x\right)\right]}^{†}$, will be called the second-level adjoint function. The 2nd-LASS is independent of parameter variations $\mathsf{\delta }{\alpha }$ and variations $V^{\left(2\right)}\left(2;x\right)$ in the respective state functions. It is also important to note that the ${\left(2\times TD\right)}^2$-dimensional matrix $A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]$ is independent of the index $j_1$. Only the source term $Q_A^{\left(2\right)}\left[2;j_1;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]$ depends on the index $j_1$. Therefore, the same solver can be used to invert the $A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]$ and numerically solve the 2nd-LASS for each $j_1$-dependent source $Q_A^{\left(2\right)}\left[2;j_1;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]$ in order to obtain the corresponding $j_1$-dependent $2\times TD$-dimensional second-level adjoint function (column vector) $A^{(2)}\left(2;j_1;x\right)\triangleq {\left[a^{(2)}\left(1;j_1;x\right),a^{(2)}\left(2;j_1;x\right)\right]}^{†}$. Computationally, it would be efficient to store, if possible, the inverse matrix ${\left\{A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]\right\}}^{-1}$, in order to multiply the inverse matrix ${\left\{A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]\right\}}^{-1}$ directly with the corresponding source term $Q_A^{\left(2\right)}\left[2;j_1;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]$, for each index $j_1$, in order to obtain the corresponding $j_1$-dependent $2\times TD$-dimensional second-level adjoint function $A^{(2)}\left(2;j_1;x\right)\triangleq {\left[a^{(2)}\left(1;j_1;x\right),a^{(2)}\left(2;j_1;x\right)\right]}^{†}$. The two components $a^{(2)}\left(1;j_1;x\right)$ and $a^{(2)}\left(2;j_1;x\right)$ of the second-level adjoint function $A^{(2)}\left(2;j_1;x\right)\triangleq {\left[a^{(2)}\left(1;j_1;x\right),a^{(2)}\left(2;j_1;x\right)\right]}^{†}$ are distinguished from each other by the use of the numbers “1” and, respectively, “2” in the respective list of arguments. In this particularly simple case, the numbers “1” and “2” could also be used as subscripts, in the customary notation for vector components, but such a use would not lend itself to generalizations because the subscript notation would become unwieldy for the higher-level adjoint functions, which will be introduced in the sections that follow below.

Using the underlying the 2nd-LASS and Equation (68) in Equation (47) yields the following expression for the indirect-effect term, for each $j_1=1,\dots ,TP$:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha };{\mathbf{v}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\delta }{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)\right]\right\}}_{ind}\\ ={\left\{{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right),{\mathbf{Q}}_V^{\left(2\right)}{\left[2;U^{\left(2\right)}\left(2;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]}_2\right\}}_{{\mathsf{\alpha }}^0}\\ -{\left\{{\left[{\widehat{P}}^{\left(2\right)}\left[{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right);{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}\\ \equiv {\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;{\mathbf{u}}^{\left(2\right)}\left(\mathbf{x}\right);{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{ind}.\end{array}$$

(73)

As the identity in Equation (73) indicates, the dependence of the indirect-effect term on the functions $v^{\left(1\right)}\left(x\right)$ and $\mathsf{\delta }{a}^{\left(1\right)}\left(x\right)$ is replaced by the dependence on the second-level adjoint function $A^{(2)}\left(2;j_1;x\right)$, for each $j_1=1,\dots ,TP$.

Replacing the expression obtained in Equation (73) for the indirect-effect term together with the expression for the direct-effect term provided in Equation (46) yields the following expression for the total differential defined by Equation (45):

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}={\left\{\frac{\partial R^{\left(1\right)}\left[j_1;u\left(x\right);a^{\left(1\right)}\left(x\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\mathsf{\delta }\mathsf{\alpha }\right\}}_{{\mathsf{\alpha }}^0}\\ +{\left\{{〈A^{(2)}\left(2;j_1;x\right),Q_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]〉}_2\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(2\right)}\left(U^{\left(2\right)};A^{(2)};\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}.\end{array}$$

(74)

Using Equations (25), (46), and (55) in Equation (74) yields the following component form for the total differential expressed by Equation (74):

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}\\ ={\displaystyle \sum _{j_2=1}^{TP}{\left\{R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}}\mathsf{\delta }{\mathsf{\alpha }}_{j_2},j_1=1,\dots ,TP,\end{array}$$

(75)

where the quantity ${\left\{R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}$ denotes the second-order sensitivity of the generic scalar-valued response $R\left[u\left(x\right);\mathsf{\alpha }\right]$ with respect to the parameters ${\mathsf{\alpha }}_{j_1}$ and ${\mathsf{\alpha }}_{j_2}$ computed at the nominal values of the parameters and respective state functions. The expression of the second-order sensitivity $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]$ of the response $R\left[u\left(x\right);\mathsf{\alpha }\right]$ with respect to two model parameters ${\mathsf{\alpha }}_{j_1}$ and ${\mathsf{\alpha }}_{j_2}$ has the following expression:

$$\begin{array}{l}For{j}_1,j_2=1,\dots ,TP:R^{\left(2\right)}\left[j_2;j_1;{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right);{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]\triangleq \frac{\partial R^{\left(1\right)}\left[j_1;\mathbf{u}\left(\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_2}}\\ -\frac{{\left\{\partial {\widehat{P}}^{\left(2\right)}\left[{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right);{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]\right\}}_{\partial {\mathsf{\Omega }}_x}}{\partial {\mathsf{\alpha }}_{j_2}}+{〈{\mathbf{a}}^{\left(2\right)}\left(1;j_1;\mathbf{x}\right),\frac{\partial \left[\mathbf{Q}\left(\mathsf{\alpha }\right)-\mathbf{N}\left(\mathbf{u};\mathsf{\alpha }\right)\right]}{\partial {\mathsf{\alpha }}_{j_2}}〉}_1\\ +{〈{\mathbf{a}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right),\frac{{\mathbf{s}}^{\left(1\right)}\left[\mathbf{u}\left(\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_2}}-\frac{\partial \left[{\mathbf{A}}^{\left(1\right)}\left(\mathbf{u};\mathsf{\alpha }\right){\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)\right]}{\partial {\mathsf{\alpha }}_{j_2}}〉}_1\triangleq \frac{{\partial }^{2}R\left[\mathbf{u}\left(\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_2}\partial {\mathsf{\alpha }}_{j_1}}.\end{array}$$

(76)

As Equations (70) and (72) indicate, solving the 2nd-LASS once provides the second-level adjoint function $a^{(2)}\left(j_1;x\right)$, for each index $j_1=1,\dots ,TP$, which enables the exact and efficient computation of the G differential $\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$ for each index $j_1=1,\dots ,TP$. Notably, the G differential $\mathsf{\delta }{R}^{\left(1\right)}\left[j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$ comprises one complete row (running on the index $j_1=1,\dots ,TP$) of mixed second-order partial sensitivities $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]={\partial }^{2}R/\partial {\mathsf{\alpha }}_{j_2}\partial {\mathsf{\alpha }}_{j_1}$. Thus, the exact computation of all of the partial second-order sensitivities, ${\left\{{\partial }^{2}R\left(\mathsf{\alpha },u\right)/\partial {\mathsf{\alpha }}_j\partial {\mathsf{\alpha }}_i\right\}}_{\left({\mathsf{\alpha }}^0\right)},$ $i,j=1,\dots ,TP,$ requires at most $TP$ large-scale (adjoint) computations using the 2nd-LASS, rather than at least $O\left(T{P}^2\right)$ large-scale computations, which would be required by forward methods.

Since the adjoint matrix $A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]$ is block diagonal, solving the 2nd-LASS is equivalent to solving two 1st-LASS, with two different source terms. Thus, the “solvers” and the computer program used for solving the 1st-LASS can also be used for solving the 2nd-LASS. The 2nd-LASS was designated as the “second-level” rather than the “second-order” adjoint sensitivity system, since the 2nd-LASS does not involve any explicit second-order G derivatives of the operators underlying the original system but involves the inversion of the same operators that needed to be inverted for solving the 1st-LASS.

It is important to note that if the 2nd-LASS is solved $TP$-times, the second-order mixed sensitivities $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]\equiv {\partial }^{2}R/\partial {\mathsf{\alpha }}_{j_2}\partial {\mathsf{\alpha }}_{j_1}$ will be computed twice, in two different ways, in terms of two distinct second-level adjoint functions. Consequently, the symmetry property ${\partial }^{2}R\left[u\left(x\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_2}\partial {\mathsf{\alpha }}_{j_1}={\partial }^{2}R\left[u\left(x\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_1}\partial {\mathsf{\alpha }}_{j_2}$ enjoyed by the second-order sensitivities provides an intrinsic (numerical) verification that the components of the second-level adjoint function $A^{(2)}\left(2;j_1;x\right)$, as well as the first-level adjoint function $a^{(1)}\left(x\right)$ are computed accurately.

The structure of the 2nd-LASS enables full flexibility for prioritizing the computation of the second-order sensitivities. The computation of the second-order sensitivities would logically be prioritized based on the relative magnitudes of the first-order sensitivities: the largest relative first-order response sensitivity should have the highest priority for computing the corresponding second-order mixed sensitivities; then, the second largest relative first-order response sensitivity should be considered next, and so on. The unimportant second-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them. Computing second-order sensitivities that correspond to vanishing first-order sensitivities may also be of interest, since vanishing first-order sensitivities may indicate critical points of the response in the phase space of model parameters.

3.3. The Third-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (3rd-CASAM-N)

The second-order sensitivities $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]\equiv {\partial }^{2}R/\partial {\mathsf{\alpha }}_{j_2}\partial {\mathsf{\alpha }}_{j_1}$ will be assumed to satisfy the conditions stated in Equations (11) and (12) for each $j_1,j_2=1,\dots ,TP$, so that the first-order total G differential of $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]$ will exist and will be linear in the variations $V^{\left(2\right)}\left(2;x\right)$ and $\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)$ in a neighborhood around the nominal values of the parameters and the respective state functions. By definition, the first-order total G differential of $R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]$, which will be denoted as ${\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };V^{\left(2\right)}\left(2;x\right);\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right);\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}$, is given by the following expression:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right);{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right);\mathsf{\alpha };{\mathbf{V}}^{\left(2\right)}\left(2;\mathbf{x}\right);\mathsf{\delta }{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right);\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}\\ \triangleq {\left\{\frac{d}{d\mathsf{\epsilon }}\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right)+\mathsf{\epsilon }{\mathbf{V}}^{\left(2\right)}\left(2;\mathbf{x}\right);{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right)+\mathsf{\epsilon }\mathsf{\delta }{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right);\mathsf{\alpha }+\mathsf{\epsilon }\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{\mathsf{\epsilon }=0}\\ \triangleq {\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right);{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}\\ +{\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right);{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right);\mathsf{\alpha };{\mathbf{V}}^{\left(2\right)}\left(2;\mathbf{x}\right);\mathsf{\delta }{A}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right)\right]\right\}}_{ind}.\end{array}$$

(77)

In Equation (77), the quantity ${\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}$ denotes the “direct-effect term”, which comprises all of the dependencies on the vector $\mathsf{\delta }\mathsf{\alpha }$ of parameter variations and has the following expression:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}\\ \triangleq {\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\mathsf{\delta }\mathsf{\alpha }\right\}}_{{\mathsf{\alpha }}^0}\\ ={\displaystyle \sum _{j_3=1}^{TP}{\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_3}}\right\}}_{{\mathsf{\alpha }}^0}}\mathsf{\delta }{\mathsf{\alpha }}_{j_3}.\end{array}$$

(78)

In Equation (77), the quantity ${\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };V^{\left(2\right)}\left(2;x\right);\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)\right]\right\}}_{ind}$ denotes the “indirect-effect term” which comprises all of the dependencies on the vectors $V^{\left(2\right)}\left(2;x\right)$ and $\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)$ of variations in the state functions $U^{\left(2\right)}\left(2;x\right)$ and $A^{(2)}\left(2;j_1;x\right)$; this indirect-effect term is defined as follows:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };V^{\left(2\right)}\left(2;x\right);\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)\right]\right\}}_{ind}\\ \triangleq {\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)};A^{(2)};\mathsf{\alpha }\right]}{\partial U^{\left(2\right)}\left(2;x\right)}\right\}}_{{\mathsf{\alpha }}^0}{V}^{\left(2\right)}\left(2;x\right)+{\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)};A^{(2)};\mathsf{\alpha }\right]}{\partial A^{(2)}\left(2;j_1;x\right)}\right\}}_{{\mathsf{\alpha }}^0}\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right),\end{array}$$

(79)

where

$$\begin{array}{l}{\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)};A^{(2)};\mathsf{\alpha }\right]}{\partial U^{\left(2\right)}\left(2;x\right)}\right\}}_{{\mathsf{\alpha }}^0}{V}^{\left(2\right)}\left(2;x\right)\triangleq {\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)};A^{(2)};\mathsf{\alpha }\right]}{\partial u\left(x\right)}\right\}}_{{\mathsf{\alpha }}^0}{v}^{\left(1\right)}\left(x\right)\\ +{\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)};A^{(2)};\mathsf{\alpha }\right]}{\partial a^{\left(1\right)}\left(x\right)}\right\}}_{{\mathsf{\alpha }}^0}\mathsf{\delta }{a}^{\left(1\right)}\left(x\right),\end{array}$$

(80)

and

$$\frac{\partial R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)};a^{\left(2\right)};\mathsf{\alpha }\right]}{\partial A^{(2)}\left(2;j_1;x\right)}\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)\triangleq {\displaystyle \sum _{i=1}^2\frac{\partial R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)};a^{\left(2\right)};\mathsf{\alpha }\right]}{\partial a^{(2)}\left(i;j_1;x\right)}\mathsf{\delta }{a}^{(2)}\left(i;j_1;x\right)}.$$

(81)

The direct-effect term ${\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}$ can be computed immediately; however, ${\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };V^{\left(2\right)}\left(2;x\right);\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)\right]\right\}}_{ind}$ can be computed only after having determined the vectors $V^{\left(2\right)}\left(2;x\right)$ and $\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)$. The vector $V^{\left(2\right)}\left(2;x\right)$ is the solution of the 2nd-LVSS defined by Equations (57) and (58). The vector $\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)$ is the solution of the G-differentiated 2nd-LASS. Taking the G differential of Equations (70) and (72) yields the following system of equations for $\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)$:

$$\begin{array}{l}\frac{\partial \left\{\mathbf{A}{\mathbf{M}}^{\left(2\right)}\left[2\times 2;{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right);\mathsf{\alpha }\right]{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right)\right\}}{\partial {\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right)}{\mathbf{V}}^{\left(2\right)}\left(2;\mathbf{x}\right)-\frac{\partial {\mathbf{Q}}_A^{\left(2\right)}\left[2;j_1;{\mathbf{u}}^{\left(2\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial {\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right)}{\mathbf{V}}^{\left(2\right)}\left(2;\mathbf{x}\right)\\ +\mathbf{A}{\mathbf{M}}^{\left(2\right)}\left[2\times 2;{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right);\mathsf{\alpha }\right]\mathsf{\delta }{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right)=\frac{\partial {\mathbf{Q}}_A^{\left(2\right)}\left[2;j_1;{\mathbf{u}}^{\left(2\right)}\left(\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\partial \mathsf{\alpha }\\ -\frac{\partial \left\{\mathbf{A}{\mathbf{M}}^{\left(2\right)}\left[2\times 2;{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right);\mathsf{\alpha }\right]{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right)\right\}}{\partial \mathsf{\alpha }}\partial \mathsf{\alpha },j_1=1,\dots ,TP;\mathbf{x}\in {\mathsf{\Omega }}_x,\end{array}$$

(82)

$${\left\{\mathsf{\delta }{B}_A^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[2\right];j_1=1,\dots ,TP;x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right).$$

(83)

The quantities which appear in Equation (83) are evaluated at the nominal values of the parameters and respective state functions, but the notation ${\left\{\right\}}_{{\mathsf{\alpha }}^0}$, which indicates this evaluation, is omitted in order to simplify the notation.

Concatenating Equations (82) and (83) with the 2nd-LVSS represented by Equations (57) and (58) yields the following system of, which will be called the “3rd-order variational sensitivity system” (3rd-LVSS), for determining the vectors $V^{\left(2\right)}\left(2;x\right)$ and $\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)$:

$$\begin{array}{l}{\left\{V{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]V^{\left(3\right)}\left(4;j_1;x\right)\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{Q_V^{\left(3\right)}\left[4;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0},x\in {\mathsf{\Omega }}_x,\end{array}$$

(84)

$${\left\{B_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;j_1;x\right);V^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[4\right];0\left[4\right]\triangleq {\left[0,0,0,0\right]}^{†},x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right)$$

(85)

where

$$\begin{array}{l}V{M}^{\left(3\right)}\left(4\times 4;U^{\left(3\right)};\mathsf{\alpha }\right)\triangleq \left(\begin{array}{cc}V{M}^{\left(2\right)}\left(2\times 2\right)& 0\left[2\times 2\right]\\ V{M}_{21}^{\left(3\right)}\left(2\times 2\right)& V{M}_{22}^{\left(3\right)}\left(2\times 2\right)\end{array}\right);\\ U^{\left(3\right)}\left(4;j_1;x\right)\triangleq \left(\begin{array}{c}{U}^{\left(2\right)}\left(2;x\right)\\ A^{(2)}\left(2;j_1;x\right)\end{array}\right);\end{array}$$

(86)

$$\begin{array}{l}{V}^{\left(3\right)}\left(4;j_1;x\right)\triangleq \mathsf{\delta }{U}^{\left(3\right)}\left(4;j_1;x\right)=\left(\begin{array}{c}{V}^{\left(2\right)}\left(2;x\right)\\ \mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)\end{array}\right)\\ ={\left[v^{\left(1\right)}\left(x\right),\mathsf{\delta }{a}^{\left(1\right)}\left(x\right),\mathsf{\delta }{a}^{\left(2\right)}\left(1;j_1;x\right),\mathsf{\delta }{a}^{\left(2\right)}\left(2;j_1;x\right)\right]}^{†};\end{array}$$

(87)

$$\begin{array}{l}V{M}_{21}^{\left(3\right)}\left(2\times 2;x\right)\triangleq \frac{\partial \left\{A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]A^{(2)}\left(2;j_1;x\right)\right\}}{\partial U^{\left(2\right)}\left(2;x\right)}\\ -\frac{\partial Q_A^{\left(2\right)}\left[2;j_1;u^{\left(2\right)}\left(x\right);\mathsf{\alpha }\right]}{\partial U^{\left(2\right)}\left(2;x\right)};\end{array}$$

(88)

$$V{M}_{22}^{\left(3\right)}\left(2\times 2;x\right)\triangleq A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right];0\left[2\times 2\right]\triangleq \left(\begin{array}{cc}0& 0\\ 0& 0\end{array}\right);$$

(89)

$$\begin{array}{l}{Q}_V^{\left(3\right)}\left[4;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\triangleq \left(\begin{array}{c}{Q}_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\\ Q_2^{\left(3\right)}\left[2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\end{array}\right)\\ \triangleq {\left\{q_V^{\left(3\right)}\left[1;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right],\dots ,q_V^{\left(3\right)}\left[4;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}^{†};\end{array}$$

(90)

$$\begin{array}{l}{Q}_2^{\left(3\right)}\left[2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\triangleq \frac{\partial Q_A^{\left(2\right)}\left[2;j_1;u^{\left(2\right)}\left(x\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\partial \mathsf{\alpha }\\ -\frac{\partial \left\{A{M}^{\left(2\right)}\left[2\times 2;U^{\left(2\right)}\left(2;x\right);\mathsf{\alpha }\right]A^{(2)}\left(2;j_1;x\right)\right\}}{\partial \mathsf{\alpha }}\partial \mathsf{\alpha };\end{array}$$

(91)

$$B_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;j_1;x\right);V^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\triangleq \left(\begin{array}{c}{B}_V^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);V^{\left(2\right)}\left(2;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\\ \mathsf{\delta }{B}_A^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]\end{array}\right).$$

(92)

For subsequent reference, it is important to note that the quantities $q_V^{\left(3\right)}\left[i;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$ are linear in the parameter variations $\mathsf{\delta }{\mathsf{\alpha }}_i$, $i=1,\dots ,TP$, and can therefore be written in the following form:

$$q_V^{\left(3\right)}\left[i;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\triangleq {\displaystyle \sum _{j_3=1}^{TP}{s}_V^{\left(3\right)}\left[i;j_3;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]}\mathsf{\delta }{\mathsf{\alpha }}_{j_3}$$

(93)

The variational matrix $V{M}^{\left(3\right)}\left(4\times 4;U^{\left(3\right)};\mathsf{\alpha }\right)$ comprises $\left(4\times 4\right)$ block matrices, each comprising $T{D}^2$ components/elements; thus, the matrix $V{M}^{\left(3\right)}\left(4\times 4;U^{\left(3\right)};\mathsf{\alpha }\right)$ comprises a total of $\left(4\times 4\right)T{D}^2$ components/elements. Each of the vectors $V^{\left(3\right)}\left(4;j_1;x\right)$, $Q_V^{\left(3\right)}\left[4;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$ and $B_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(2;j_1;x\right);V^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$ comprise four $TD$-dimensional vectors, as shown in their respective definitions; thus, each of the vectors $V^{\left(3\right)}\left(4;j_1;x\right)$, $Q_V^{\left(3\right)}\left[4;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$, and $B_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(2;j_1;x\right);V^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$ comprise $\left(4\times TD\right)$ components/elements.

Solving the 3rd-LVSS would require $T{P}^3$ large-scale computations, which is unrealistic for large-scale systems comprising many parameters. The 3rd-CASAM-N circumvents the need for solving the 3rd-LVSS by deriving an alternative expression for the indirect-effect term defined in Equation (81), in which the function $V^{\left(3\right)}\left(4;j_1;x\right)$ is replaced by a third-level adjoint function which is independent of parameter variations. This third-level adjoint function will be the solution of a third-level adjoint sensitivity system (3rd-LASS) which will be constructed by applying the same principles as those used for constructing the 1st-LASS and the 2nd-LASS. The Hilbert space appropriate for constructing the 3rd-LASS will be denoted as ${\mathrm{H}}_3\left({\mathsf{\Omega }}_x\right)$ and comprise element block vectors of the same form as $V^{\left(3\right)}\left(4;j_1;x\right)$. Thus, a generic block vector in ${\mathrm{H}}_3\left({\mathsf{\Omega }}_x\right)$, denoted as ${\Psi }^{(3)}\left(4;x\right)\triangleq {\left[{\mathsf{\psi }}^{(3)}\left(1;x\right),{\mathsf{\psi }}^{(3)}\left(2;x\right),{\mathsf{\psi }}^{(3)}\left(3;x\right),{\mathsf{\psi }}^{(3)}\left(4;x\right)\right]}^{†}\in {\mathrm{H}}_3\left({\mathsf{\Omega }}_x\right)$, comprises four $TD$-dimensional vector components of the form ${\mathsf{\psi }}^{\left(3\right)}\left(i;x\right)\triangleq {\left[{\mathsf{\psi }}_1^{\left(3\right)}\left(i;x\right),\dots ,{\mathsf{\psi }}_{TD}^{\left(3\right)}\left(i;x\right)\right]}^{†}\in {\mathrm{H}}_1\left({\mathsf{\Omega }}_x\right)$, $i=1,2,3,4$, where each of these four components is a $TD$-dimensional column vector. The inner product of two vectors ${\Psi }^{(3)}\left(4;x\right)\in {\mathrm{H}}_3\left({\mathsf{\Omega }}_x\right)$ and ${\Phi }^{(3)}\left(4;x\right)\in {\mathrm{H}}_3\left({\mathsf{\Omega }}_x\right)$ in the Hilbert space ${\mathrm{H}}_3\left({\mathsf{\Omega }}_x\right)$ will be denoted as ${〈{\Psi }^{(3)}\left(4;x\right),{\Phi }^{(3)}\left(4;x\right)〉}_3$ and defined as follows:

$${〈{\Psi }^{(3)}\left(4;x\right),{\Phi }^{(3)}\left(4;x\right)〉}_3\triangleq {\displaystyle \sum _{i=1}^4{〈{\mathsf{\psi }}^{(3)}\left(i;x\right),{\mathsf{\phi }}^{(3)}\left(i;x\right)〉}_1}$$

(94)

The inner product defined in Equation (94) is continuous in $\mathsf{\alpha }$, in a neighborhood around ${\mathsf{\alpha }}^0$. Using the definition of the inner product defined in Equation (94), the inner product of Equation (84) is constructed with a vector $A^{(3)}\left(4;x\right)\triangleq {\left[a^{(3)}\left(1;x\right),a^{(3)}\left(2;x\right),a^{(3)}\left(3;x\right),a^{(3)}\left(4;x\right)\right]}^{†}\in {\mathrm{H}}_3\left({\mathsf{\Omega }}_x\right)$ to obtain the following relation:

$$\begin{array}{l}{\left\{{〈A^{(3)}\left(4;x\right),V{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]V^{\left(3\right)}\left(4;j_1;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{{〈A^{(3)}\left(4;x\right),Q_V^{\left(3\right)}\left[4;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]〉}_3\right\}}_{{\mathsf{\alpha }}^0},x\in {\mathsf{\Omega }}_x.\end{array}$$

(95)

The inner product on the left side of Equation (95) is further transformed by using the definition of the adjoint operator to obtain the following relation:

$$\begin{array}{l}{\left\{{〈A^{(3)}\left(4;x\right),V{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]V^{\left(3\right)}\left(4;j_1;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{{〈V^{\left(3\right)}\left(4;j_1;x\right),A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]A^{(3)}\left(4;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}\\ +{\left\{{\left[P^{\left(3\right)}\left(U^{\left(3\right)};A^{(3)};V^{\left(3\right)};\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0},\end{array}$$

(96)

where

$$\begin{array}{l}A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]\triangleq {\left[V{M}^{\left(3\right)}\left(4\times 4;U^{\left(3\right)};\mathsf{\alpha }\right)\right]}^{\ast }\\ =\left(\begin{array}{cc}{\left\{{\left[V{M}^{\left(2\right)}\left(2\times 2\right)\right]}^{\ast }\right\}}^{†}& {\left\{{\left[V{M}_{21}^{\left(3\right)}\left(2\times 2\right)\right]}^{\ast }\right\}}^{†}\\ 0\left[2\times 2\right]& {\left\{{\left[V{M}_{22}^{\left(3\right)}\left(2\times 2\right)\right]}^{\ast }\right\}}^{†}\end{array}\right),\end{array}$$

(97)

and where ${\left\{{\left[P^{\left(3\right)}\left(U^{\left(3\right)};A^{(3)};V^{\left(3\right)};{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ denotes the corresponding bilinear concomitant on the domain’s boundary, evaluated at the nominal values for the parameters and respective state functions. The adjoint matrix $A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]$ comprises ${\left(4\times TD\right)}^2$ components/elements, while the adjoint function $A^{(3)}\left(4;x\right)\in {\mathrm{H}}_3\left({\mathsf{\Omega }}_x\right)$ comprises $\left(4\times TD\right)$ components/elements.

The domain of the adjoint matrix operator $A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]$ is specified by requiring that the function $A^{(3)}\left(4;x\right)\triangleq {\left[a^{(3)}\left(1;x\right),a^{(3)}\left(2;x\right),a^{(3)}\left(3;x\right),a^{(3)}\left(4;x\right)\right]}^{†}\in {\mathrm{H}}_3\left({\mathsf{\Omega }}_x\right)$ satisfies adjoint boundary/initial conditions denoted as follows:

$${\left\{B_A^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;j_1;x\right);A^{(3)}\left(4;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[4\right],x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right).$$

(98)

The third-level adjoint boundary/initial conditions represented by Equation (98) are determined by requiring that: (a) they must be independent of unknown values of $V^{\left(3\right)}\left(4;j_1;x\right)$; and (b) the substitution of the boundary and/or initial conditions represented by Equations (85) and (98) into the expression of ${\left\{{\left[P^{\left(3\right)}\left(U^{\left(3\right)};A^{(3)};V^{\left(3\right)};\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ must cause all terms containing unknown values of $V^{\left(3\right)}\left(4;j_1;x\right)$ to vanish.

Implementing the boundary/initial conditions represented by Equations (85) and (98) into Equation (96) will transform the latter into the following form:

$$\begin{array}{l}{\left\{{〈V^{\left(3\right)}\left(4;j_1;x\right),A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]A^{(3)}\left(4;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{{〈A^{(3)}\left(4;x\right),V{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]V^{\left(3\right)}\left(4;j_1;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}\\ -{\left\{{\left[{\widehat{P}}^{\left(3\right)}\left(U^{\left(3\right)};A^{(3)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0},\end{array}$$

(99)

where ${\left\{{\left[{\widehat{P}}^{\left(3\right)}\left(U^{\left(3\right)};A^{(3)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ denotes residual boundary terms which may not vanish automatically. The right side of Equation (95) is now used in the first term on the right side of Equation (68) to obtain the following relation:

$$\begin{array}{l}{\left\{{〈V^{\left(3\right)}\left(4;j_1;x\right),A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]A^{(3)}\left(4;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{{〈A^{(3)}\left(4;x\right),Q_V^{\left(3\right)}\left[4;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]〉}_3\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(3\right)}\left(U^{\left(3\right)};A^{(3)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}.\end{array}$$

(100)

The definition of the third-level adjoint function $A^{(3)}\left(4;x\right)\in {\mathrm{H}}_3\left({\mathsf{\Omega }}_x\right)$ is now completed by requiring that the left side of Equation (100) and the right side of Equation (81) represent the “indirect-effect term” ${\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha };V^{\left(2\right)}\left(2;x\right);\mathsf{\delta }{A}^{(2)}\left(2;j_1;x\right)\right]\right\}}_{ind}$ for each of the indices $j_1=1,\dots ,TP;j_2=1,....,j_1$. Hence, there will be $TP\left(TP+1\right)/2$ distinct third-level adjoint functions $A^{(3)}\left(4;j_2;j_1;x\right)\triangleq {\left[a^{(3)}\left(1;j_2;j_1;x\right),a^{(3)}\left(2;j_2;j_1;x\right),a^{(3)}\left(3;j_2;j_1;x\right),a^{(3)}\left(4;j_2;j_1;x\right)\right]}^{†}$, corresponding to the indices $j_1=1,\dots ,TP;j_2=1,....,j_1$. Each of these distinct third-level adjoint functions will correspond to a specific $\left(j_1,j_2\right)$-dependent indirect-effect term.

The left side of Equation (100) will be identical to the right side of Equation (81) by requiring that the following relation be satisfied by the third-level adjoint functions $A^{(3)}\left(4;j_2;j_1;x\right)\triangleq {\left[a^{(3)}\left(1;j_2;j_1;x\right),a^{(3)}\left(2;j_2;j_1;x\right),a^{(3)}\left(3;j_2;j_1;x\right),a^{(3)}\left(4;j_2;j_1;x\right)\right]}^{†}$:

$$\begin{array}{l}{\left\{A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]A^{(3)}\left(4;j_2;j_1;x\right)\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{Q_A^{\left(3\right)}\left[4;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0},j_1=1,\dots ,TP;j_2=1,....,j_1,\end{array}$$

(101)

$$Q_A^{\left(3\right)}\left[4;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]\triangleq \left[q_A^{\left(3\right)}\left(1;j_2;j_1;U^{\left(3\right)};\mathsf{\alpha }\right),\dots ,q_A^{\left(3\right)}\left(1;j_2;j_1;U^{\left(3\right)};\mathsf{\alpha }\right)\right];$$

(102)

$$q_A^{\left(3\right)}\left(1;j_2;j_1;U^{\left(3\right)};\mathsf{\alpha }\right)\triangleq \partial R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)};a^{\left(2\right)};\mathsf{\alpha }\right]/\partial u^{\left(1\right)};$$

(103)

$$q_A^{\left(3\right)}\left(2;j_2;j_1;U^{\left(3\right)};\mathsf{\alpha }\right)\triangleq \partial R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)};a^{\left(2\right)};\mathsf{\alpha }\right]/\partial a^{\left(1\right)};$$

(104)

$$q_A^{\left(3\right)}\left(3;j_2;j_1;U^{\left(3\right)};\mathsf{\alpha }\right)\triangleq \partial R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)};a^{\left(2\right)};\mathsf{\alpha }\right]/\partial a^{(2)}\left(1;j_1;x\right);$$

(105)

$$q_A^{\left(3\right)}\left(3;j_2;j_1;U^{\left(3\right)};\mathsf{\alpha }\right)\triangleq \partial R^{\left(2\right)}\left[j_2;j_1;u^{\left(2\right)};a^{\left(2\right)};\mathsf{\alpha }\right]/\partial a^{(2)}\left(2;j_1;x\right).$$

(106)

The boundary conditions to be satisfied by each of the third-level adjoint functions $A^{(3)}\left(4;j_2;j_1;x\right)\triangleq {\left[a^{(3)}\left(1;j_2;j_1;x\right),a^{(3)}\left(2;j_2;j_1;x\right),a^{(3)}\left(3;j_2;j_1;x\right),a^{(3)}\left(4;j_2;j_1;x\right)\right]}^{†}$ are those represented by Equation (98), namely

$$\begin{array}{l}{\left\{B_A^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;j_1;x\right);A^{(3)}\left(4;j_2;j_1;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[4\right];\\ j_1=1,\dots ,TP;j_2=1,\dots ,j_1;x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right).\end{array}$$

(107)

Since the source $Q_A^{\left(3\right)}\left[4;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]$ may contain distributions, the equality in Equation (101) is considered to hold in the weak sense. The Riesz representation theorem ensures that the weak equality in Equation (101) holds uniquely.

The system of represented by Equations (101) and (107) will be called the third-level adjoint sensitivity system (3rd-LASS); its solution, $A^{(3)}\left(4;j_2;j_1;x\right)$, will be called the third-level adjoint function. Using the underlying the 3rd-LASS and Equation (100) in Equation (81) yields the following expression for the indirect-effect term:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;{\mathbf{U}}^{\left(2\right)}\left(2;\mathbf{x}\right);{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right);\mathsf{\alpha };{\mathbf{V}}^{\left(2\right)}\left(2;\mathbf{x}\right);\mathsf{\delta }{\mathbf{A}}^{\left(2\right)}\left(2;j_1;\mathbf{x}\right)\right]\right\}}_{ind}\\ ={\left\{{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right),{\mathbf{Q}}_V^{\left(3\right)}{\left[4;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]}_3\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(3\right)}\left({\mathbf{U}}^{\left(3\right)};{\mathbf{A}}^{\left(3\right)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}\\ \equiv {\left\{{\delta }{R}^{\left(2\right)}\left[j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]\right\}}_{ind};j_1=1,\dots ,TP;j_2=1,\dots ,j_1.\end{array}$$

(108)

As the identity in Equation (108) indicates, the dependence of the indirect-effect term on the function $V^{\left(3\right)}\left(4;j_1;x\right)$ is replaced by the dependence on the adjoint function $A^{(3)}\left(4;j_2;j_1;x\right)$, for each $j_1=1,\dots ,TP;j_2=1,\dots ,j_1$.

Replacing the expression obtained in Equation (108) for the indirect-effect term together with the expression for the direct-effect term provided in Equation (78) yields the following expression for the total differential defined by Equation (77):

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{\frac{\partial R^{\left(2\right)}\left[j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\mathsf{\delta }\mathsf{\alpha }\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(3\right)}\left({\mathbf{U}}^{\left(3\right)};{\mathbf{A}}^{\left(3\right)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}\\ +{\left\{{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right),{\mathbf{Q}}_V^{\left(3\right)}{\left[4;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]}_3\right\}}_{{\mathsf{\alpha }}^0}.\end{array}$$

(109)

In component form, the total differential expressed by Equation (109) has the following expression:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{(3)}\left(4;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}\\ ={\displaystyle \sum _{j_3=1}^{TP}{\left\{R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{(3)}\left(4;j_2;j_1;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}}\mathsf{\delta }{\mathsf{\alpha }}_{j_3},j_1;j_2=1,\dots ,TP,\end{array}$$

(110)

where the quantity $R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{(3)}\left(4;j_2;j_1;x\right);\mathsf{\alpha }\right]$ denotes the third-order sensitivity of the generic scalar-valued response $R\left[u\left(x\right);\mathsf{\alpha }\right]$ with respect to any three model parameters ${\mathsf{\alpha }}_{j_1}$, ${\mathsf{\alpha }}_{j_2}$, ${\mathsf{\alpha }}_{j_3}$, and has the following expression:

$$\begin{array}{l}{R}^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]\triangleq \frac{\partial R^{\left(2\right)}\left[j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_3}}\\ -\frac{{\left[\partial {\widehat{P}}^{\left(3\right)}\left({\mathbf{U}}^{\left(3\right)};{\mathbf{A}}^{\left(3\right)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}}{\partial {\mathsf{\alpha }}_{j_3}}+{\displaystyle {\displaystyle \sum }_{i=1}^4}{〈{\mathbf{a}}^{\left(3\right)}\left(i;j_2;j_1;\mathbf{x}\right),{\mathbf{s}}_V^{\left(3\right)}\left[i;j_3;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]〉}_1\\ \triangleq \frac{{\partial }^{3}R\left[\mathbf{u}\left(\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_3}\partial {\mathsf{\alpha }}_{j_2}\partial {\mathsf{\alpha }}_{j_1}},for{j}_1,j_2,j_3=1,\dots ,TP.\end{array}$$

(111)

As Equations (101)–(107) indicate, solving the 3rd-LASS once provides the third-level adjoint function $A^{(3)}\left(4;j_2;j_1;x\right)$, for each of the indices $j_1=1,\dots ,TP;j_2=1,\dots ,j_1$. In turn, the availability of $A^{(3)}\left(4;j_2;j_1;x\right)$ enables the exact and efficient computation of the G differential $\mathsf{\delta }{R}^{\left(2\right)}\left[j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{(3)}\left(4;j_2;j_1;x\right);\mathsf{\alpha }\right]$. Thus, the exact computation of all of the partial third-order sensitivities ${\left\{{\partial }^{3}R\left(\mathsf{\alpha },u\right)/\partial {\mathsf{\alpha }}_{j_3}\partial {\mathsf{\alpha }}_{j_2}\partial {\mathsf{\alpha }}_{j_1}\right\}}_{\left({\mathsf{\alpha }}^0\right)}$ requires at most $TP\left(TP+1\right)/2$ large-scale (adjoint) computations using the 3rd-LASS, rather than at least $O\left(T{P}^3\right)$ large-scale computations, which would be required by forward methods.

The matrix $A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]$ is block diagonal; therefore, solving the 3rd-LASS is equivalent to solving three 1st-LASS, with different source terms. The 3rd-LASS was designated as “the third-level” rather than “third-order” adjoint sensitivity system since the 3rd-LASS does not involve any explicit second-order and/or third-order G derivatives of the operators underlying the original system, but involves only the inversion of the same operators that needed to be inverted for solving the 1st-LASS.

By solving the 3rd-LASS $TP\left(TP+1\right)/2$ times, the third-order mixed sensitivities ${\left\{{\partial }^{3}R\left(\mathsf{\alpha },u\right)/\partial {\mathsf{\alpha }}_{j_3}\partial {\mathsf{\alpha }}_{j_2}\partial {\mathsf{\alpha }}_{j_1}\right\}}_{\left({\mathsf{\alpha }}^0\right)}$ will be computed three times, in three different ways. Consequently, the multiple symmetries intrinsic to the third-order sensitivities provide an intrinsic numerical verification that the components of the first-, second-, and third-level adjoint functions are computed accurately.

The structure of the 3rd-LASS enables full flexibility for prioritizing the computation of the third-order sensitivities. The computation of the third-order sensitivities would logically be prioritized based on the relative magnitudes of the second-order sensitivities, so that the unimportant third-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them.

3.4. The Fourth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (4th-CASAM-N)

The third-order sensitivities $R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha }\right]\equiv {\partial }^{3}R\left[u\left(x\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_3}\partial {\mathsf{\alpha }}_{j_2}\partial {\mathsf{\alpha }}_{j_1}$ will be assumed to satisfy the conditions stated in Equations (11) and (12) for each $j_1,j_2,j_3=1,\dots ,TP$, so that the first-order total G differential of $R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha }\right]$ will exist and will be linear in the variations $V^{\left(3\right)}\left(4;j_1;x\right)\triangleq \mathsf{\delta }{U}^{\left(3\right)}\left(4;j_1;x\right)$ and $\mathsf{\delta }{A}^{(3)}\left(4;j_2;j_1;x\right)$ in a neighborhood around the nominal values of the parameters and the respective state functions. By definition, the first-order total G differential of $R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{(3)}\left(4;j_2;j_1;x\right);\mathsf{\alpha }\right]$, which will be denoted as ${\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);A^{(3)}\left(4;j_2;j_1;x\right);\mathsf{\alpha };V^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\delta }{A}^{(3)}\left(4;j_2;j_1;x\right);\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}$, is given by the following expression:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha };{\mathbf{V}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\delta }{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}\\ \triangleq \left\{\frac{d}{d\mathsf{\epsilon }}\right.R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right)+\mathsf{\epsilon }{\mathbf{V}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\right.\\ {\left.{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right)+\mathsf{\epsilon }\mathsf{\delta }{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha }+\mathsf{\epsilon }\mathsf{\delta }\mathsf{\alpha }\right\}}_{\mathsf{\epsilon }=0}\\ \triangleq {\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}\\ +{\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha };{\mathbf{V}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\delta }{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right)\right]\right\}}_{ind}.\end{array}$$

(112)

In Equation (112), the quantity ${\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}$ denotes the “direct-effect term” which comprises all of the dependencies on the vector $\mathsf{\delta }\mathsf{\alpha }$ of parameter variations and has the following expression:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}\triangleq \\ {\left\{\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\mathsf{\delta }\mathsf{\alpha }\right\}}_{{\mathsf{\alpha }}^0}\\ ={\displaystyle {\displaystyle \sum }_{j_4=1}^{TP}}{\left\{\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_4}}\right\}}_{{\mathsf{\alpha }}^0}\mathsf{\delta }{\mathsf{\alpha }}_{j_4}.\end{array}$$

(113)

In Equation (112), the quantity ${\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha };V^{\left(3\right)};\mathsf{\delta }{A}^{(3)}\right]\right\}}_{ind}$ denotes the “indirect-effect term” which comprises all of the dependencies on the vectors $V^{\left(3\right)}\left(4;j_1;x\right)$ and $\mathsf{\delta }{A}^{(3)}\left(4;j_2;j_1;x\right)$ of variations in the state functions $U^{\left(3\right)}\left(4;j_1;x\right)$ and $A^{(3)}\left(4;j_2;j_1;x\right)$; this indirect-effect term is defined as follows:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha };{\mathbf{V}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\delta }{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right)\right]\right\}}_{ind}\\ \triangleq {\left\{\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)};{\mathbf{A}}^{\left(3\right)};\mathsf{\alpha }\right]}{\partial {\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right)}\right\}}_{{\mathsf{\alpha }}^0}{\mathbf{V}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right)\\ +{\left\{\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)};{\mathbf{A}}^{\left(3\right)};\mathsf{\alpha }\right]}{\partial {\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right)}\right\}}_{{\mathsf{\alpha }}^0}\mathsf{\delta }{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right),\end{array}$$

(114)

where

$$\begin{array}{l}\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)};{\mathbf{A}}^{\left(3\right)};\mathsf{\alpha }\right]}{\partial {\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right)}{\mathbf{V}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right)\\ \triangleq \frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)};{\mathbf{A}}^{\left(3\right)};\mathsf{\alpha }\right]}{\partial \mathbf{u}\left(\mathbf{x}\right)}{\mathbf{v}}^{\left(1\right)}\left(\mathbf{x}\right)+\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)};{\mathbf{A}}^{\left(3\right)};\mathsf{\alpha }\right]}{\partial {\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)}\mathsf{\delta }{\mathbf{a}}^{\left(1\right)}\left(\mathbf{x}\right)\\ +{\displaystyle {\displaystyle \sum }_{i=1}^2}\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)};{\mathbf{A}}^{\left(3\right)};\mathsf{\alpha }\right]}{\partial {\mathbf{a}}^{\left(2\right)}\left(i;j_1;\mathbf{x}\right)}\mathsf{\delta }{\mathbf{a}}^{\left(2\right)}\left(i;j_1;\mathbf{x}\right),\end{array}$$

(115)

and

$$\begin{array}{l}\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha }\right]}{\partial A^{(3)}\left(4;j_2;j_1;x\right)}\mathsf{\delta }{A}^{(3)}\left(4;j_2;j_1;x\right)\\ \triangleq {\displaystyle \sum _{i=1}^4\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha }\right]}{\partial a^{(3)}\left(i;j_2;j_1;x\right)}\mathsf{\delta }{a}^{(3)}\left(i;j_2;j_1;x\right)}.\end{array}$$

(116)

The direct-effect term ${\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{dir}$ can be computed immediately; however, the indirect-effect term ${\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha };V^{\left(3\right)};\mathsf{\delta }{A}^{(3)}\right]\right\}}_{ind}$ can be computed only after having determined the vectors $V^{\left(3\right)}\left(4;j_1;x\right)$ and $\mathsf{\delta }{A}^{(3)}\left(4;j_2;j_1;x\right)$. The vector $V^{\left(3\right)}\left(4;j_1;x\right)$ is the solution of the 3rd-LVSS defined by Equations (84) and (85). The vector $\mathsf{\delta }{A}^{(3)}\left(4;j_2;j_1;x\right)$ is the solution of the G-differentiated 3rd-LASS. Taking the G differential of Equations (101) and (107) yields the following system of for $\mathsf{\delta }{A}^{(3)}\left(4;j_2;j_1;x\right)$ and for $j_1=1,\dots ,TP;j_2=1,\dots ,j_1$:

$$\begin{array}{l}\frac{\partial \left\{\mathbf{A}{\mathbf{M}}^{\left(3\right)}\left[4\times 4;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right)\right\}}{\partial {\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right)}{\mathbf{V}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right)\\ -\frac{\partial {\mathbf{Q}}_A^{\left(3\right)}\left[4;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial {\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right)}{\mathbf{V}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right)\\ +\mathbf{A}{\mathbf{M}}^{\left(3\right)}\left[4\times 4;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]\mathsf{\delta }{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right)\\ =\frac{\partial {\mathbf{Q}}_A^{\left(3\right)}\left[4;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\partial \mathsf{\alpha }\\ -\frac{\partial \left\{\mathbf{A}{\mathbf{M}}^{\left(3\right)}\left[4\times 4;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right)\right\}}{\partial \mathsf{\alpha }}\partial \mathsf{\alpha },\mathbf{x}\in {\mathsf{\Omega }}_x,\end{array}$$

(117)

$${\left\{\mathsf{\delta }{B}_A^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;j_1;x\right);A^{(3)}\left(4;j_2;j_1;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[4\right];x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right).$$

(118)

The quantities which appear in Equation (117) are evaluated at the nominal values of the parameters and respective state functions, but the notation ${\left\{\right\}}_{{\mathsf{\alpha }}^0}$, which indicates this evaluation, is omitted in order to simplify the notation.

Concatenating Equations (117) and (118) with the 3rd-LVSS represented by Equations (84) and (85) yields the following system of, which will be called the “4th-order variational sensitivity system” (4th-LVSS), for determining the vectors $V^{\left(3\right)}\left(4;j_1;x\right)$ and $\mathsf{\delta }{A}^{(3)}\left(4;j_2;j_1;x\right)$, for $j_1=1,\dots ,TP;j_2=1,\dots ,j_1$:

$$\begin{array}{l}{\left\{V{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]V^{\left(4\right)}\left(8;j_2;j_1;x\right)\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{Q_V^{\left(4\right)}\left[8;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0},x\in {\mathsf{\Omega }}_x,\end{array}$$

(119)

$${\left\{B_V^{\left(4\right)}\left[8;U^{\left(4\right)}\left(8;j_1;x\right);V^{\left(4\right)}\left(8;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[8\right];x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right);$$

(120)

where $0\left[8\right]\triangleq {\left[0,0,0,0,0,0,0,0\right]}^{†}$, and

$$V{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]\triangleq \left(\begin{array}{cc}V{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]& 0\left[4\times 4\right]\\ V{M}_{21}^{\left(4\right)}\left(4\times 4;x\right)& V{M}_{22}^{\left(4\right)}\left(4\times 4;x\right)\end{array}\right);$$

(121)

$$U^{\left(4\right)}\left(8;j_2;j_1;x\right)\triangleq \left(\begin{array}{c}{U}^{\left(3\right)}\left(4;j_1;x\right)\\ A^{(3)}\left(4;j_2;j_1;x\right)\end{array}\right);$$

(122)

$$\begin{array}{l}{V}^{\left(4\right)}\left(8;j_2;j_1;x\right)\triangleq \mathsf{\delta }{U}^{\left(4\right)}\left(8;j_2;j_1;x\right)=\left(\begin{array}{c}{V}^{\left(3\right)}\left(4;j_1;x\right)\\ \mathsf{\delta }{A}^{(3)}\left(4;j_2;j_1;x\right)\end{array}\right)\\ =\left[v^{\left(1\right)}\left(x\right),\mathsf{\delta }{a}^{\left(1\right)}\left(x\right),\mathsf{\delta }{a}^{\left(2\right)}\left(1;j_1;x\right),\mathsf{\delta }{a}^{\left(2\right)}\left(2;j_1;x\right)\right.,\\ {\left.\mathsf{\delta }{a}^{(3)}\left(1;j_2;j_1;x\right),\dots ,\mathsf{\delta }{a}^{(3)}\left(4;j_2;j_1;x\right)\right]}^{†};\end{array}$$

(123)

$$\begin{array}{l}V{M}_{21}^{\left(4\right)}\left(4\times 4;x\right)\triangleq -\frac{\partial Q_A^{\left(3\right)}\left[4;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]}{\partial U^{\left(3\right)}\left(4;j_1;x\right)}\\ +\frac{\partial \left\{A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]A^{(3)}\left(4;j_2;j_1;x\right)\right\}}{\partial U^{\left(3\right)}\left(4;j_1;x\right)};\end{array}$$

(124)

$$V{M}_{22}^{\left(4\right)}\left(4\times 4;x\right)\triangleq A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right];$$

(125)

$$\begin{array}{l}{Q}_V^{\left(4\right)}\left[8;j_1;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\triangleq \left(\begin{array}{c}{Q}_V^{\left(3\right)}\left[4;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\\ Q_2^{\left(4\right)}\left[4;j_1;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\end{array}\right)\\ \equiv {\left\{q_V^{\left(4\right)}\left[1;j_1;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right],\dots ,q_V^{\left(4\right)}\left[8;j_1;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\right\}}^{†};\end{array}$$

(126)

$$\begin{array}{l}{Q}_2^{\left(4\right)}\left[4;j_1;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\triangleq \frac{\partial Q_A^{\left(3\right)}\left[4;j_2;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\partial \mathsf{\alpha }\\ -\frac{\partial \left\{A{M}^{\left(3\right)}\left[4\times 4;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]A^{(3)}\left(4;j_2;j_1;x\right)\right\}}{\partial \mathsf{\alpha }}\partial \mathsf{\alpha };\end{array}$$

(127)

$$\begin{array}{l}{B}_V^{\left(4\right)}\left[8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);V^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\\ \triangleq \left(\begin{array}{c}{B}_V^{\left(3\right)}\left[4;U^{\left(3\right)}\left(4;j_1;x\right);V^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\\ \mathsf{\delta }{B}_A^{\left(2\right)}\left[2;U^{\left(2\right)}\left(2;x\right);A^{(2)}\left(2;j_1;x\right);\mathsf{\alpha }\right]\end{array}\right).\end{array}$$

(128)

For subsequent reference, it is noted that the quantities $q_V^{\left(4\right)}\left[i;j_2;j_1;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$ are linear in the parameter variations $\mathsf{\delta }{\mathsf{\alpha }}_i$, $i=1,\dots ,TP$, and can therefore be written in the following form:

$$q_V^{\left(4\right)}\left[i;j_2;j_1;U^{\left(4\right)}\left(4;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]\triangleq {\displaystyle \sum _{j_4=1}^{TP}{s}_V^{\left(3\right)}\left[i;j_4;j_2;j_1;U^{\left(4\right)}\left(4;x\right);\mathsf{\alpha }\right]}\mathsf{\delta }{\mathsf{\alpha }}_{j_4};i=1,\dots ,8.$$

(129)

The variational matrix $V{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]$ comprises $\left(8\times 8\right)$ block matrices, each comprising $T{D}^2$ components/elements; thus, the matrix $V{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]$ comprises a total of $\left(8\times 8\right)T{D}^2$ components/elements. Each of the vectors $V^{\left(4\right)}\left(8;j_2;j_1;x\right)$, $Q_V^{\left(4\right)}\left[8;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$, and $B_V^{\left(4\right)}\left[8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);V^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$ comprise eight $TD$-dimensional vectors, as shown in their respective definitions; thus, each of the vectors $B_V^{\left(4\right)}\left[8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);V^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$, $V^{\left(4\right)}\left(8;j_2;j_1;x\right)$, $Q_V^{\left(4\right)}\left[8;j_1;U^{\left(3\right)}\left(4;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]$ comprise $\left(8\times TD\right)$ components/elements.

Solving the 4th-LVSS would require $T{P}^4$ large-scale computations, which is unrealistic for large-scale systems comprising many parameters. The 4th-CASAM-N circumvents the need for solving the 3rd-LVSS by deriving an alternative expression for the indirect-effect term defined in Equation (114), in which the function $V^{\left(4\right)}\left(8;j_2;j_1;x\right)$ is replaced by a fourth-level adjoint function which is independent of parameter variations. This fourth-level adjoint function will be the solution of a fourth-level adjoint sensitivity system (4th-LASS) which will be constructed by applying the same principles as those used for constructing the 1st-LASS, the 2nd-LASS, and the 3rd-LASS. The Hilbert space appropriate for constructing the 4th-LASS will be denoted as ${\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$ and comprises block vectors elements of the same form as $V^{\left(4\right)}\left(8;j_2;j_1;x\right)$. Thus, a generic block vector in ${\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$ will have the structure ${\Psi }^{(4)}\left(8;x\right)\triangleq {\left[{\mathsf{\psi }}^{(4)}\left(1;x\right),\dots ,{\mathsf{\psi }}^{(4)}\left(8;x\right)\right]}^{†}\in {\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$, comprising eight $TD$-dimensional vectors of the form ${\mathsf{\psi }}^{\left(4\right)}\left(i;x\right)\triangleq {\left[{\mathsf{\psi }}_1^{\left(4\right)}\left(i;x\right),\dots ,{\mathsf{\psi }}_{TD}^{\left(4\right)}\left(i;x\right)\right]}^{†}\in {\mathrm{H}}_1\left({\mathsf{\Omega }}_x\right)$, $i=1,\dots ,8$. The inner product of two vectors ${\Psi }^{(4)}\left(8;x\right)\in {\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$ and ${\Phi }^{(4)}\left(8;x\right)\in {\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$ in the Hilbert space ${\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$ will be denoted as ${〈{\Psi }^{(4)}\left(8;x\right),{\Phi }^{(4)}\left(8;x\right)〉}_4$ and defined as follows:

$${〈{\Psi }^{(4)}\left(8;x\right),{\Phi }^{(4)}\left(8;x\right)〉}_4\triangleq {\displaystyle \sum _{i=1}^8{〈{\mathsf{\psi }}^{(3)}\left(i;x\right),{\mathsf{\phi }}^{(3)}\left(i;x\right)〉}_1}$$

(130)

The inner product defined in Equation (130) is continuous in $\mathsf{\alpha }$, in a neighborhood around ${\mathsf{\alpha }}^0$. Using the definition of the inner product defined in Equation (130), construct the inner product of Equation (119) with a vector $A^{(4)}\left(8;x\right)\triangleq {\left[a^{(4)}\left(1;x\right),\dots ,a^{(4)}\left(8;x\right)\right]}^{†}\in {\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$ to obtain the following relation:

$$\begin{array}{l}{\left\{{〈A^{(4)}\left(8;x\right),V{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]V^{\left(4\right)}\left(8;j_2;j_1;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{{〈A^{(4)}\left(8;x\right),Q_V^{\left(4\right)}\left[8;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]〉}_3\right\}}_{{\mathsf{\alpha }}^0},x\in {\mathsf{\Omega }}_x.\end{array}$$

(131)

The inner product on the left side of Equation (131) is further transformed by using the definition of the adjoint operator to obtain the following relation:

$$\begin{array}{l}{\left\{{〈A^{(4)}\left(8;x\right),V{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]V^{\left(4\right)}\left(8;j_2;j_1;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}=\\ {\left\{{〈V^{\left(4\right)}\left(8;j_2;j_1;x\right),A{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]A^{(4)}\left(8;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}\\ +{\left\{{\left[P^{\left(4\right)}\left(U^{\left(4\right)};A^{(4)};V^{\left(4\right)};\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0},\end{array}$$

(132)

where

$$\begin{array}{l}A{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]\triangleq {\left[V{M}^{\left(4\right)}\left(8\times 8;U^{\left(4\right)};\mathsf{\alpha }\right)\right]}^{\ast }\\ =\left(\begin{array}{cc}{\left\{{\left[V{M}^{\left(3\right)}\left(4\times 4\right)\right]}^{\ast }\right\}}^{†}& {\left\{{\left[V{M}_{21}^{\left(4\right)}\left(4\times 4\right)\right]}^{\ast }\right\}}^{†}\\ 0\left[4\times 4\right]& {\left\{{\left[V{M}_{22}^{\left(3\right)}\left(4\times 4\right)\right]}^{\ast }\right\}}^{†}\end{array}\right),\end{array}$$

(133)

and where ${\left\{{\left[P^{\left(4\right)}\left(U^{\left(4\right)};A^{(4)};V^{\left(4\right)};\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ denotes the corresponding bilinear concomitant on the domain’s boundary, evaluated at the nominal values for the parameters and respective state functions. The adjoint matrix $A{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]$ comprises $\left(8\times 8\right)T{D}^2$ components/elements, while the adjoint function $A^{(4)}\left(8;x\right)\in {\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$ comprises $\left(8\times TD\right)$ components/elements.

The domain of the adjoint matrix operator $A{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]$ is specified by requiring that the function $A^{(4)}\left(8;x\right)\triangleq {\left[a^{(4)}\left(1;x\right),\dots ,a^{(4)}\left(8;x\right)\right]}^{†}\in {\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$ satisfies adjoint boundary/initial conditions denoted as follows:

$${\left\{B_A^{\left(4\right)}\left[8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{(4)}\left(8;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[8\right],x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right).$$

(134)

The fourth-level adjoint boundary/initial conditions represented by Equation (134) are determined by requiring that: (a) they must be independent of unknown values of $V^{\left(4\right)}\left(8;j_2;j_1;x\right)$; and (b) the substitution of the boundary and/or initial conditions represented by Equations (85) and (98) into the expression of ${\left\{{\left[P^{\left(4\right)}\left(U^{\left(4\right)};A^{(4)};V^{\left(4\right)};\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ must cause all terms containing unknown values of $V^{\left(4\right)}\left(8;j_2;j_1;x\right)$ to vanish.

Implementing the boundary/initial conditions represented by Equations (128) and (134) into Equation (132) will transform the later relation into the following form:

$$\begin{array}{l}{\left\{{〈V^{\left(4\right)}\left(8;j_2;j_1;x\right),A{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]A^{(4)}\left(8;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}\\ =-{\left\{{\left[{\widehat{P}}^{\left(4\right)}\left(U^{\left(4\right)};A^{(4)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}\\ +{\left\{{〈A^{(4)}\left(8;x\right),V{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]V^{\left(4\right)}\left(8;j_2;j_1;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0},\end{array}$$

(135)

where ${\left\{{\left[{\widehat{P}}^{\left(4\right)}\left(U^{\left(4\right)};A^{(4)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}$ denotes residual boundary terms which may have not vanish automatically. The right side of Equation (131) is now used in the first term on the right side of Equation (135) to obtain the following relation:

$$\begin{array}{l}{\left\{{〈V^{\left(4\right)}\left(8;j_2;j_1;x\right),A{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]A^{(4)}\left(8;x\right)〉}_3\right\}}_{{\mathsf{\alpha }}^0}=\\ {\left\{{〈A^{(4)}\left(8;x\right),Q_V^{\left(4\right)}\left[8;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]〉}_3\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(4\right)}\left(U^{\left(4\right)};A^{(4)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}.\end{array}$$

(136)

The definition of the fourth-level adjoint function $A^{(4)}\left(8;x\right)\triangleq {\left[a^{(1)}\left(1;x\right),\dots ,a^{(4)}\left(8;x\right)\right]}^{†}\in {\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$ is now completed by requiring that the left side of Equation (136) and the right side of Equation (114) represent the “indirect-effect term” ${\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha };V^{\left(3\right)};\mathsf{\delta }{A}^{(3)}\right]\right\}}_{ind}$, for each of the indices $j_1=1,\dots ,TP$; $j_2=1,....,j_1$; $j_3=1,....,j_2$. Hence, there will be $TP\left(TP+1\right)\left(TP+2\right)/6$ distinct fourth-level adjoint functions $A^{(4)}\left(8;x\right)\in {\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$, each corresponding to one combination of the indices $j_1=1,\dots ,TP$; $j_2=1,....,j_1$; $j_3=1,....,j_2$. Each of these distinct fourth-level adjoint functions will correspond to a specific $\left(j_1,j_2,j_3\right)$-dependent indirect-effect term.

The left side of Equation (136) will be identical to the right side of Equation (114) by requiring that the following relation be satisfied by the fourth-level adjoint functions $A^{(4)}\left(8;j_3;j_2;j_1;x\right)\triangleq {\left[a^{(4)}\left(1;j_3;j_2;j_1;x\right),\dots ,a^{(4)}\left(8;j_3;j_2;j_1;x\right)\right]}^{†}\in {\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$, for each value of the indices $j_1=1,\dots ,TP$; $j_2=1,....,j_1$; $j_3=1,....,j_2$:

$$\begin{array}{l}A{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]A^{(4)}\left(8;j_3;j_2;j_1;x\right)\\ =Q_A^{\left(4\right)}\left[8;j_3;j_2;j_1;U^{\left(4\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right],\end{array}$$

(137)

where

$$\begin{array}{l}{Q}_A^{\left(4\right)}\left[8;j_3;j_2;j_1;U^{\left(4\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]\triangleq \left\{q_A^{\left(4\right)}\left[1;j_3;j_2;j_1;U^{\left(4\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]\right.,\dots \\ {\left.\dots ,q_A^{\left(4\right)}\left[8;j_3;j_2;j_1;U^{\left(4\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]\right\}}^{†},j_1=1,\dots ,TP;j_2=1,\dots ,j_1;j_3=1,\dots ,j_2;\end{array}$$

(138)

$$q_A^{\left(4\right)}\left[1;j_3;j_2;j_1;U^{\left(4\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]\triangleq \partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha }\right]/\partial u\left(x\right);$$

(139)

$$q_A^{\left(4\right)}\left[2;j_3;j_2;j_1;U^{\left(4\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]\triangleq \partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha }\right]/\partial a^{\left(1\right)}\left(x\right);$$

(140)

$$q_A^{\left(4\right)}\left[2+i;j_3;j_2;j_1;U^{\left(4\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]\triangleq \frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha }\right]}{\partial a^{\left(2\right)}\left(i;j_1;x\right)};i=1,2;$$

(141)

$$q_A^{\left(4\right)}\left[4+i;j_3;j_2;j_1;U^{\left(4\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]\triangleq \frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(3\right)};A^{(3)};\mathsf{\alpha }\right]}{\partial a^{(3)}\left(i;j_2;j_1;x\right)};i=1,2,3,4;$$

(142)

The boundary conditions to be satisfied by each of the fourth-level adjoint functions $A^{(4)}\left(8;j_3;j_2;j_1;x\right)\triangleq {\left[a^{(4)}\left(1;j_3;j_2;j_1;x\right),\dots ,a^{(4)}\left(8;j_3;j_2;j_1;x\right)\right]}^{†}\in {\mathrm{H}}_4\left({\mathsf{\Omega }}_x\right)$ are those represented by Equation (134), namely

$$\begin{array}{l}{\left\{B_A^{\left(4\right)}\left[8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{(4)}\left(8;j_3;j_2;j_1;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}=0\left[8\right],x\in \partial {\mathsf{\Omega }}_x\left({\mathsf{\alpha }}^0\right);\\ j_1=1,\dots ,TP;j_2=1,\dots ,j_1;j_3=1,\dots ,j_2.\end{array}$$

(143)

Since the source $Q_A^{\left(4\right)}\left[8;j_3;j_2;j_1;U^{\left(4\right)}\left(4;j_1;x\right);\mathsf{\alpha }\right]$ may contain distributions, the equality in Equation (137) is considered to hold in the weak sense, and the Riesz representation theorem ensures that this weak equality holds uniquely.

The system of represented by Equations (101) and (107) will be called the fourth-level adjoint sensitivity system (4th-LASS); its solution, $A^{(4)}\left(8;j_3;j_2;j_1;x\right)$, will be called the fourth-level adjoint function. Using the underlying the 4th-LASS and Equation (136) in Equation (114) yields the following expression for the indirect-effect term:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha };{\mathbf{V}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);\mathsf{\delta }{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right)\right]\right\}}_{ind}\\ \triangleq {\left\{{〈{\mathbf{A}}^{\left(4\right)}\left(8;\mathbf{x}\right),{\mathbf{Q}}_V^{\left(4\right)}\left[8;j_1;{\mathbf{U}}^{\left(4\right)}\left(8;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]〉}_3\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(4\right)}\left({\mathbf{U}}^{\left(4\right)};{\mathbf{A}}^{\left(4\right)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}\\ \equiv {\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(4\right)}\left(8;j_2;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(4\right)}\left(8;j_3;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]\right\}}_{ind}.\end{array}$$

(144)

As the identity in Equation (144) indicates, the dependence of the indirect-effect term on the function $V^{\left(4\right)}\left(8;j_2;j_1;x\right)$ is replaced by the dependence on the adjoint function $A^{(4)}\left(8;j_3;j_2;j_1;x\right)$ for each $j_1=1,\dots ,TP;$ $j_2=1,\dots ,j_1$; $j_3=1,\dots ,j_2$.

Replacing the expression obtained in Equation (144) for the indirect-effect term together with the expression for the direct-effect term provided in Equation (113) yields the following expression for the total differential defined by Equation (112):

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(4\right)}\left(8;j_2;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(4\right)}\left(8;j_3;j_2;j_1;\mathbf{x}\right);\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}\\ ={\left\{\frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial \mathsf{\alpha }}\mathsf{\delta }\mathsf{\alpha }\right\}}_{{\mathsf{\alpha }}^0}-{\left\{{\left[{\widehat{P}}^{\left(4\right)}\left({\mathbf{U}}^{\left(4\right)};{\mathbf{A}}^{\left(4\right)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}\right\}}_{{\mathsf{\alpha }}^0}\\ +{\left\{{〈{\mathbf{A}}^{\left(4\right)}\left(8;j_3;j_2;j_1;\mathbf{x}\right),{\mathbf{Q}}_V^{\left(4\right)}\left[8;j_1;{\mathbf{U}}^{\left(4\right)}\left(8;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha };\mathsf{\delta }\mathsf{\alpha }\right]〉}_3\right\}}_{{\mathsf{\alpha }}^0}\end{array}$$

(145)

In component form, the total differential expressed by Equation (145) can be written in the following form, for each $j_1=1,\dots ,TP;$ $j_2=1,\dots ,j_1$; $j_3=1,\dots ,j_2$:

$$\begin{array}{l}{\left\{\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{(4)}\left(8;j_3;j_2;j_1;x\right);\mathsf{\delta }\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}\\ ={\displaystyle \sum _{j_4=1}^{TP}{\left\{R^{\left(4\right)}\left[j_4;j_3;j_2;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{(4)}\left(8;j_3;j_2;j_1;x\right);\mathsf{\alpha }\right]\right\}}_{{\mathsf{\alpha }}^0}}\mathsf{\delta }{\mathsf{\alpha }}_{j_4},\end{array}$$

(146)

where the quantity $R^{\left(4\right)}\left[j_4;j_3;j_2;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{(4)}\left(8;j_3;j_2;j_1;x\right);\mathsf{\alpha }\right]$ denotes the fourth-order sensitivity of the generic scalar-valued response $R\left[u\left(x\right);\mathsf{\alpha }\right]$ with respect to any four model parameters ${\mathsf{\alpha }}_{j_1}$, ${\mathsf{\alpha }}_{j_2}$, ${\mathsf{\alpha }}_{j_3}$, ${\mathsf{\alpha }}_{j_3}$, and has the following expression:

$$\begin{array}{l}{R}^{\left(4\right)}\left[j_4;j_3;j_2;j_1;{\mathbf{U}}^{\left(4\right)}\left(8;j_2;j_1;\mathbf{x}\right);A^{\left(4\right)}\left(8;j_3;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]\\ \triangleq \frac{\partial R^{\left(3\right)}\left[j_3;j_2;j_1;{\mathbf{U}}^{\left(3\right)}\left(4;j_1;\mathbf{x}\right);{\mathbf{A}}^{\left(3\right)}\left(4;j_2;j_1;\mathbf{x}\right);\mathsf{\alpha }\right]}{\partial {\mathsf{\alpha }}_{j_4}}-\frac{{\left[\partial {\widehat{P}}^{\left(4\right)}\left({\mathbf{U}}^{\left(4\right)};{\mathbf{A}}^{\left(4\right)};\mathsf{\delta }\mathsf{\alpha }\right)\right]}_{\partial {\mathsf{\Omega }}_x}}{\partial {\mathsf{\alpha }}_{j_4}}\\ +{\displaystyle {\displaystyle \sum }_{i=1}^8}{〈{\mathbf{a}}^{\left(4\right)}\left(i;j_3;j_2;j_1;\mathbf{x}\right),{\mathbf{s}}_V^{\left(3\right)}\left[i;j_4;j_2;j_1;{\mathbf{U}}^{\left(4\right)}\left(4;\mathbf{x}\right);\mathsf{\alpha }\right]〉}_1\\ \equiv {\partial }^{4}R\left[\mathbf{u}\left(\mathbf{x}\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_1}\partial {\mathsf{\alpha }}_{j_3}\partial {\mathsf{\alpha }}_{j_3}\partial {\mathsf{\alpha }}_{j_4};j_1,j_2,j_3,j_4=1,\dots ,j_3=1,\dots ,TP.\end{array}$$

(147)

As Equations (137) and (143) indicate, solving the 4th-LASS once provides the fourth-level adjoint function $A^{(4)}\left(8;j_3;j_2;j_1;x\right)$ for each of the indices $j_1=1,\dots ,TP;$ $j_2=1,\dots ,j_1$, and $j_3=1,\dots ,j_2$. In turn, the availability of $A^{(4)}\left(8;j_3;j_2;j_1;x\right)$ enables the exact and efficient computation of the G differential $\mathsf{\delta }{R}^{\left(3\right)}\left[j_3;j_2;j_1;U^{\left(4\right)}\left(8;j_2;j_1;x\right);A^{(4)}\left(8;j_3;j_2;j_1;x\right);\mathsf{\delta }\mathsf{\alpha }\right]$. Thus, the exact computation of all of the partial fourth-order sensitivities, ${\partial }^{4}R\left[u\left(x\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_1}\partial {\mathsf{\alpha }}_{j_3}\partial {\mathsf{\alpha }}_{j_3}\partial {\mathsf{\alpha }}_{j_4}$ requires at most $TP\left(TP+1\right)\left(TP+2\right)/6$ large-scale (adjoint) computations using the 4th-LASS, rather than at least $O\left(T{P}^4\right)$ large-scale computations, which would be required by forward methods.

The adjoint matrix $A{M}^{\left(4\right)}\left[8\times 8;U^{\left(4\right)}\left(8;j_2;j_1;x\right);\mathsf{\alpha }\right]$ is block diagonal; therefore, solving the 4th-LASS is equivalent to solving four 1st-LASS, with different source terms. The 4th-LASS was designated as the “fourth-level” rather than “fourth-order” adjoint sensitivity system since the 3rd-LASS does not involve any explicit second-order, third-order, and/or fourth-order G derivatives of the operators underlying the original system but involves the inversion of the operators similar to those that needed to be inverted for solving the 1st-LASS.

By solving the 4th-LASS $TP\left(TP+1\right)\left(TP+2\right)/6$ times, the fourth-order mixed sensitivities ${\partial }^{4}R\left[u\left(x\right);\mathsf{\alpha }\right]/\partial {\mathsf{\alpha }}_{j_1}\partial {\mathsf{\alpha }}_{j_3}\partial {\mathsf{\alpha }}_{j_3}\partial {\mathsf{\alpha }}_{j_4}$ will be computed four times, in four different ways using distinct adjoint functions. Consequently, the multiple symmetries intrinsic to the fourth-order sensitivities provide an intrinsic numerical verification that the components of the first-, second-, third-, and fourth-level adjoint functions are computed accurately.

The structure of the 4th-LASS enables full flexibility for prioritizing the computation of the third-order sensitivities. The computation of the fourth-order sensitivities would be prioritized based on the relative magnitudes of the third-order sensitivities, so that the unimportant fourth-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them.

4. Discussion

This work has presented the “fourth-order comprehensive sensitivity analysis methodology for nonlinear systems (abbreviated as “4th-CASAM-N”), which enables the hitherto very difficult, if not intractable, exact computation of all of the first-, second-, third-, and fourth-order response sensitivities for large-scale nonlinear systems involving many parameters. The qualifier “comprehensive” indicates that the fourth-CASAM-N methodology enables the exact and efficient computation not only of response sensitivities with respect to the customary model parameters (including computational input data, correlations, initial and/or boundary conditions) but also with respect to imprecisely known material boundaries, which would be caused by manufacturing tolerances.

It has been shown that the first-order sensitivities of the system response under consideration with respect to the model parameters (including boundary and initial conditions) can be computed inexpensively, using just quadrature formulas, after having obtained the first-level adjoint state function. The first-level adjoint state function is obtained by solving once the first-level adjoint sensitivity system (1st-LASS), which is the sole large-scale computation needed for obtaining all of the first-order sensitivities. This is because the 1st-LASS is independent of parameter variations, and therefore needs to be solved only once, regardless of the number of model parameters (denoted as “$TP$”) under consideration. Furthermore, solving the 1st-LASS requires less computational effort than solving the nonlinear underlying the original models, since the 1st-LASS is linear in the first-level adjoint function.

The second-order sensitivities which correspond to each first-order sensitivity are also computed by using inexpensive quadrature formulas, after having obtained the second-level adjoint state function. The computation of the second-level adjoint state function requires solving the second-level adjoint sensitivity system (2nd-LASS). Since the 2nd-LASS is block diagonal, solving it is equivalent to solving two 1st-LASS, with two different source terms. Thus, the “solvers” and the computer program used for solving the 1st-LASS can also be used for solving the 2nd-LASS. The 2nd-LASS was designated as the “second-level” rather than the “second-order” adjoint sensitivity system, since the 2nd-LASS does not involve any second-order G derivatives of the operators underlying the original system but involves the inversion of the same operators that need to be inverted to solve the 1st-LASS. The structure of the 2nd-LASS enables full flexibility for prioritizing the computation of the second-order sensitivities. The computation of the second-order sensitivities would logically be prioritized based on the relative magnitudes of the first-order sensitivities: the largest relative first-order response sensitivity should have the highest priority for computing the corresponding second-order mixed sensitivities; then, the second largest relative first-order response sensitivity should be considered next, and so on. The unimportant second-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them. Computing second-order sensitivities that correspond to vanishing first-order sensitivities may also be of interest, since vanishing first-order sensitivities may indicate critical points of the response in the phase space of model parameters. If the 2nd-LASS is solved $TP$ times, the second-order mixed sensitivities will be computed twice, in two different ways, in terms of two distinct second-level adjoint functions. Consequently, the symmetry property enjoyed by the second-order sensitivities provides an intrinsic (numerical) verification that the components of the first- and second-level adjoint functions are computed accurately.

The exact computation of all of the partial third-order sensitivities that correspond to a second-order sensitivities is also accomplished by using quadrature formulas after having obtained the third-level adjoint function. Thus, for each second-order sensitivity, the third-level adjoint function is obtained by solving the third-level adjoint sensitivity system (3rd-LASS) once, which is equivalent to solving three 1st-LASS, with different source terms. The 3rd-LASS is designated as “the third-level” rather than the “third-order” adjoint sensitivity system since the 3rd-LASS does not involve any explicit second-order and/or third-order G derivatives of the operators underlying the original system. By solving the 3rd-LASS $TP\left(TP+1\right)/2$ times, the third-order mixed sensitivities will be computed three times, in three different ways. Consequently, the multiple symmetries intrinsic to the third-order sensitivities provide an intrinsic numerical verification that the components of the first-, second- and third-level adjoint functions are computed accurately. The structure of the 3rd-LASS enables full flexibility for prioritizing the computation of the third-order sensitivities. The computation of the third-order sensitivities would logically be prioritized based on the relative magnitudes of the second-order sensitivities, so that the unimportant third-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them.

The exact computation of all of the partial fourth-order sensitivities that correspond to a third-order sensitivity also involves the use of quadrature formulas after having obtained the fourth-level adjoint function by solving the fourth-level adjoint sensitivity system (4th-LASS). The 4th-LASS is designated as the “fourth-level” rather than “fourth-order” adjoint sensitivity system since the 3rd-LASS does not involve any explicit second-order, third-order, and/or fourth-order G derivatives of the operators underlying the original system but involves the inversion of the operators similar to those that needed to be inverted for solving the 1st-LASS. The 4th-LASS is block diagonal; solving it is equivalent to solving four 1st-LASS, with different source terms. If the 4th-LASS is solved $TP\left(TP+1\right)\left(TP+2\right)/6$ times, the fourth-order mixed-response sensitivities will be computed four times, in four different ways using distinct adjoint functions. Consequently, the multiple symmetries intrinsic to the fourth-order sensitivities provide an intrinsic numerical verification that the components of the first-, second-, third-, and fourth-level adjoint functions are computed accurately. The structure of the 4th-LASS enables full flexibility for prioritizing the computation of the third-order sensitivities. The computation of the fourth-order sensitivities would be prioritized based on the relative magnitudes of the third-order sensitivities, so that the unimportant fourth-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them.

In summary, the implementation of the 4th-CASAM-N requires very little additional effort beyond the construction of the adjoint sensitivity system needed for computing the first-order sensitivities. An illustrative application of the 4th-CASAM-N to a paradigm nonlinear heat conduction is presented in the accompanying work [17].