First-Order Comprehensive Adjoint Method for Computing Operator-Valued Response Sensitivities to Imprecisely Known Parameters, Internal Interfaces and Boundaries of Coupled Nonlinear Systems: I. Mathematical Framework

Abstract

This work presents the first-order comprehensive adjoint sensitivity analysis methodology (1st-CASAM) for computing efficiently the first-order sensitivities (i.e., functional derivatives) of operator-valued responses (i.e., model results) of general models of coupled nonlinear physical systems characterized by imprecisely known or and/or uncertain parameters, external boundaries, and internal interfaces between the coupled systems. The explicit mathematical formalism developed within the 1st-CASAM for computing the first-order sensitivities of operator-valued response to uncertain internal interfaces and external boundaries in the models’ phase–space enables this methodology to generalize all of the previously published methodologies for computing first-order response sensitivities. The computational resources needed for using forward versus adjoint operators in conjunction with spectral versus collocation methods for computing the response sensitivities are analyzed in detail. By enabling the exact computations of operator-valued response sensitivities to internal interfaces and external boundary parameters and conditions, the 1st-CASAM presented in this work makes it possible, inter alia, to quantify the effects of manufacturing tolerances on operator-valued responses of physical and engineering systems.

1. Introduction

The aim of sensitivity analysis is to compute the sensitivities (i.e., functional derivatives) of responses (i.e., results of interest) of a computational model with respect to the respective model’s parameters. Statistical and “brute force” methods can yield approximate values for such sensitivities, while forward and adjoint methods yield mathematically exact expressions for sensitivities, which can therefore be computed to machine accuracy. It is beyond the scope of this work to review these methods and their numerous applications, but the interested reader may wish to consult the books [1,2,3] and references therein. The specific aim of this work is to generalize the forward/adjoint sensitivity analysis methodology conceived by Cacuci [4,5] for operator-valued (as opposed to scalar-valued) model responses, to enable the explicit computation of operator-valued response sensitivities to uncertain phase–space locations of boundaries and interfaces in coupled nonlinear subsystems. Knowledge of such sensitivities is crucial in practice for enabling the quantification of the effects of manufacturing tolerances when actually constructing any physical system, from benchmark experiments to industrial-size installations. It will be shown that response sensitivities to the imprecisely known phase–space locations of domain boundaries and interfaces can arise both from the definition of the system’s response as well as from the equations, interfaces and boundary conditions defining the model and its imprecisely known domain of definition. This work is structured as follows: Section 2 presents the general mathematical framework for computing exactly (in parameter space) and efficiently the sensitivities of a generic operator type to the physical system’s imprecisely known parameters, internal and external boundaries. This mathematical framework is called the “first-order comprehensive adjoint sensitivity analysis methodology” (1st-CASAM), where the qualifier “comprehensive” indicates that that all possible uncertain model parameters, including those characterizing the phase–space locations of internal and external boundaries are explicitly taken into consideration. The total sensitivity of the operator-valued response is represented using its spectral expansion and, alternatively, using its collocation/pseudo-spectral expansion. The relative advantages and disadvantages are discussed, including using mixed spectral/collocation expansions of the sensitivities of the operator-valued response. Section 3 offers concluding remarks.

2. Mathematical Framework of the 1st-Casam for Operator-Valued Responses of Coupled Systems Comprising Imprecisely Known Parameters, Interfaces and Boundaries

The system considered in this work comprises two nonlinear sub-systems which are coupled to one another across a common internal interface (boundary) in phase–space, and which will be called “Subsystem I” and, respectively, “Subsystem II”. The first subsystem is represented mathematically as follows:

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

(1)

Bold letters will be used in this work to denote matrices and vectors. Unless explicitly stated otherwise, the vectors in this work are considered to be column vectors. The second subsystem is represented mathematically as follows:

$$N^{\left(II\right)}\left[v\left(y\right);\alpha \right]=Q^{\left(II\right)}\left(\alpha ;y\right)\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in {\mathsf{\Omega }}_y\left(\alpha \right)$$

(2)

If differential operators appear in Equations (1) and (2), a corresponding set of boundary and/or initial/final conditions must also be given, these conditions can be represented in operator form as follows:

$$B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]=\hspace{0.17em}0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left(\alpha \right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left(\alpha \right)$$

(3)

The quantities appearing in Equations (1)–(3) are defined as follows:

• $\alpha \triangleq {\left({\alpha }_1,\dots ,{\alpha }_{Z_{\alpha }}\right)}^{†}\in {ℝ}^{Z_{\alpha }}$ denotes a column vector having $Z_{\alpha }$ scalar-valued components representing all of the imprecisely known internal and boundary parameters of the physical systems, including imprecisely known parameters that characterize the interface and boundary conditions. Some of these parameters are common to both physical systems, e.g., the parameters that characterize common interfaces. These scalar parameters are considered to be subject to both random and systematic uncertainties, as is usually the case in practical applications. In order to use such parameters in practical computations, which is the scope of the methodology presented in this work, they are considered to be either “uncertain” or “imprecisely known”. “Uncertain” parameters are usually considered to follow a probability distribution having a known “mean value” and a known “standard deviation”. On the other hand, the actual values of “imperfectly known” parameters are unknown. To enable the use of such parameters in computations, “expert opinion” is invoked to assign each of such imprecisely known parameters a “nominal value” (which plays the role of a “mean value”) and a “range of variation” (which plays the role of a standard deviation). For practical computations, the actual origin of the parameter’s nominal (or mean) value and of its assigned standard deviation is immaterial, which is why the qualifiers “uncertain” and “imprecisely known” are often used interchangeably. In this work, the superscript “zero” will be used to denote the known nominal or mean values of various quantities. In particular, the vector of nominal and/or mean parameter values will be denoted as ${\alpha }^0\triangleq {\left({\alpha }_1^0,\dots ,{\alpha }_{Z_{\alpha }}^0\right)}^{†}$. The symbol “$\triangleq$” will be used to denote “is defined as” or “is by definition equal to” and transposition will be indicated by a dagger $(†)$ superscript.
• $x\triangleq {\left(x_1,\dots ,x_{Z_x}\right)}^{†}\in {ℝ}^{Z_x}$ denotes the phase–space position vector, of dimension $\hspace{0.17em}{Z}_x$, of independent variables for the system defined in Equation (1). The vector of independent variable $x$ is defined on a phase–space domain denoted as ${\mathsf{\Omega }}_x\left(\alpha \right)$, ${\mathsf{\Omega }}_x\left(\alpha \right)\triangleq \left\{-\infty \le a_i\left(\alpha \right)\le x_i\le b_i\left(\alpha \right)\le \infty ;\hspace{0.17em}\hspace{0.17em}i=1,\dots ,Z_x\right\}$, and is therefore considered to depend on the uncertain parameters $\alpha$. The lower-valued imprecisely known boundary-point of the independent variable is denoted as $a_i\left(\alpha \right)$, while the upper-valued imprecisely known boundary-point of the independent variable is denoted as $b_i\left(\alpha \right)$. For physical systems modeled by diffusion theory, for example, the “vacuum boundary condition” requires that the particle flux vanish at the “extrapolated boundary” of the spatial domain facing the vacuum; the “extrapolated boundary” depends on the imprecisely known geometrical dimensions of the system’s domain in space and also on the system’s microscopic transport cross sections and atomic number densities. The boundary $\partial {\mathsf{\Omega }}_x\left(\alpha \right)\triangleq \left\{a\left(\alpha \right)\cup b\left(\alpha \right)\right\}$ of the domain ${\mathsf{\Omega }}_x\left(\alpha \right)$ comprises all of the endpoints $a\left(\alpha \right)\triangleq {\left[a_1\left(\alpha \right),\dots ,a_{Z_x}\left(\alpha \right)\right]}^{†}$ and $b\left(\alpha \right)\triangleq {\left[b_1\left(\alpha \right),\dots ,b_{Z_x}\left(\alpha \right)\right]}^{†}$ of the intervals on which the respective components of $x$ are defined. It may happen that some components $a_i\left(\alpha \right)$ and/or $b_j\left(\alpha \right)$ are infinite, in which case they would not depend on any imprecisely known parameters.
• $u\left(x\right)\triangleq {\left[u_1\left(x\right),\dots ,u_{Z_u}\left(x\right)\right]}^{†}$ denotes a $Z_u$-dimensional column vector whose components represent the system’s dependent variables (also called “state functions”). The vector-valued function $u\left(x\right)$ is considered the unique nontrivial solution of the physical problem described by Equations (1) and (3).
• $N^{\left(I\right)}\left[u\left(x\right);\alpha \right]\triangleq {\left[N_1^{\left(I\right)}\left(u;\alpha \right),\dots ,N_i^{\left(I\right)}\left(u;\alpha \right),\dots ,N_{Z_u}^{\left(I\right)}\left(u;\alpha \right)\right]}^{†},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\dots ,Z_u,$ denotes a column vector of dimensions $Z_u$ whose components are operators that act nonlinearly on $u\left(x\right)$ and $\alpha$.
• $Q^{\left(I\right)}\left(\alpha ;x\right)\triangleq {\left[Q_1^{\left(I\right)}\left(\alpha ;x\right),\dots ,Q_{Z_u}^{\left(I\right)}\left(\alpha ;x\right)\right]}^{†}$ denotes a $Z_u$-dimensional column vector whose elements represent inhomogeneous source terms that depend either linearly or nonlinearly on $\alpha$. The components of $Q^{\left(I\right)}\left(\alpha ;x\right)$ may involve operators (rather than just finite-dimensional functions) and distributions acting on $\alpha$ and $x$.
• $y\triangleq {\left(y_1,\dots ,y_{Z_y}\right)}^{†}\in {ℝ}^{Z_y}$ denotes the $\hspace{0.17em}{Z}_y$-dimensional phase–space position vector of independent variables for the physical system defined in Equation (2). The vector of independent variable $y$ is defined on a phase–space domain denoted as ${\mathsf{\Omega }}_y\left(\alpha \right)$, which is defined as follows: ${\mathsf{\Omega }}_y\left(\alpha \right)\triangleq \left\{-\infty \le c_j\left(\alpha \right)\le y_j\le d_j\left(\alpha \right)\le \infty ;\hspace{0.17em}\hspace{0.17em}j=1,\dots ,Z_y\right\}$. The lower-valued imprecisely known boundary-point of the independent variable $y_i$ is denoted as $c_j\left(\alpha \right)$, while the upper-valued imprecisely known boundary-point of the independent variable $y_i$ is denoted as $d_j\left(\alpha \right)$. Some or all of the points $c_j\left(\alpha \right)$ may coincide with the points $b_j\left(\alpha \right)$. Additionally, some components of $y$ may coincide with some components of $x$, in which case the respective lower and upper boundary points for the respective coinciding independent variables would also coincide correspondingly. The boundary $\partial {\mathsf{\Omega }}_y\left(\alpha \right)\triangleq \left\{c\left(\alpha \right)\cup d\left(\alpha \right)\right\}$ of the domain ${\mathsf{\Omega }}_y\left(\alpha \right)$ comprises all of the endpoints $c\left(\alpha \right)\triangleq {\left[c_1\left(\alpha \right),\dots ,c_{Z_y}\left(\alpha \right)\right]}^{†}$ and $d\left(\alpha \right)\triangleq {\left[d_1\left(\alpha \right),\dots ,d_{Z_y}\left(\alpha \right)\right]}^{†}$ of the intervals on which the respective components of $y$ are defined.
• $v\left(y\right)\triangleq {\left[v_1\left(y\right),\dots ,v_{Z_v}\left(y\right)\right]}^{†}$ denotes a $Z_v$-dimensional column vector whose components represent the system’s dependent variables (also called “state functions”). The vector-valued function $v\left(y\right)$ is considered the unique nontrivial solution of the physical problem described by Equations (2) and (3).
• $N^{\left(II\right)}\left[u\left(x\right);\alpha \right]\triangleq {\left[N_1^{\left(II\right)}\left(u;\alpha \right),\dots ,N_i^{\left(II\right)}\left(u;\alpha \right),\dots ,N_{Z_v}^{\left(II\right)}\left(u;\alpha \right)\right]}^{†},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\dots ,Z_v,$ denotes a column vector of dimensions $Z_v$ whose components are operators acting nonlinearly on $v\left(y\right)$ and $\alpha$.
• $Q^{\left(II\right)}\left(\alpha ;y\right)\triangleq {\left[Q_1^{\left(II\right)}\left(\alpha ;y\right),\dots ,Q_{Z_v}^{\left(II\right)}\left(\alpha ;y\right)\right]}^{†}$ denotes a $Z_v$-dimensional column vector whose elements represent inhomogeneous source terms that depend either linearly or nonlinearly on $\alpha$. The components of $Q^{\left(II\right)}\left(\alpha ;y\right)$ may involve operators and distributions acting on $\alpha$ and $y$.
• The vector-valued operator $B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]$ comprises all of the boundary, interface, and initial/final conditions for the coupled physical systems. If the boundary, interface and/or initial/final conditions are inhomogeneous, which is most often the case, then $B\left[0,0;\alpha ;x,y\right]\ne \hspace{0.17em}0$.
• Since $Q^{\left(I\right)}\left(\alpha ;x\right)$ and $Q^{\left(II\right)}\left(\alpha ;y\right)$ may involve operators and distributions acting on $\alpha$ and $y$, all of the equalities in this work, including Equations (1)–(3), are considered to hold in the weak (“distributional”) sense.

The nominal (or “base-case”) solutions of Equations (1)–(3), denoted as $u^0\left(x\right)$ and $v^0\left(y\right)$, are obtained by solving these equations at the nominal parameter values ${\alpha }^0$, i.e.,

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

(4)

$$N^{\left(II\right)}\left[v^0\left(y\right);{\alpha }^0\right]=Q^{\left(II\right)}\left({\alpha }^0;y\right)\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in {\mathsf{\Omega }}_y\left({\alpha }^0\right)$$

(5)

$$B\left[u^0\left(x\right),v^0\left(y\right);{\alpha }^0;x,y\right]=\hspace{0.17em}0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left({\alpha }^0\right)$$

(6)

The response considered in this work is a generic nonlinear function-valued operator, denoted as follows:

$$R\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]$$

(7)

The nominal value of the response, denoted as $R^0\triangleq R\left[u^0\left(x\right),v^0\left(y\right);{\alpha }^0;x,y\right]$, is determined by computing the response at the nominal values ${\alpha }^0$, $u^0\left(x\right)$ and $v^0\left(y\right)$. The true values of imprecisely known model, interface and boundary parameters may differ from their nominal (average, or “base-case”) values by variations denoted as $\delta \alpha \triangleq \left(\delta {\alpha }_1,\dots ,\delta {\alpha }_{N_{\alpha }}\right)$, where $\delta {\alpha }_i\triangleq {\alpha }_i-{\alpha }_i^0$, $i=1,\dots ,N_{\alpha }$. In turn, the parameter variations $\delta \alpha$ will cause variations $\delta u\left(x\right)\triangleq {\left[\delta u_1\left(x\right),\dots ,\delta u_{Z_u}\left(x\right)\right]}^{†}$ and $\delta v\left(y\right)\triangleq {\left[\delta v_1\left(y\right),\dots ,\delta v_{Z_v}\left(y\right)\right]}^{†}$ in the state functions and all of these variations will cause variations in the response $R\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]$ around the nominal response value $R^0$. Sensitivity analysis aims at computing the functional derivatives (called “sensitivities”) of the response to the imprecisely known parameters $\alpha$. Subsequently, these sensitivities can be used for a variety of purposes, including quantifying the uncertainties induced in responses by the uncertainties in the model and boundary parameters, combining the uncertainties in computed responses with uncertainties in measured response (“data assimilation”) to obtain more accurate predictions of responses and/or parameters (“model calibration”, “predictive modeling”, etc.).

As has been shown by Cacuci [1], the most general definition of the 1st-order total sensitivity of an operator-valued model response to parameter variations is provided by the first-order “Gateaux-variation” (G-variation) of the response under consideration. To determine the first G-variation of the response $R\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]$, it is convenient to denote the functions appearing in the argument of the response as being the components of a vector $e\triangleq {\left[u\left(x\right),v\left(y\right);\alpha \right]}^{†}$, which represents an arbitrary “point” in the combined phase–space of the state functions and parameters. The point which corresponds to the nominal values of the state functions and parameters in this phase space is denoted as $e^0\triangleq {\left[u^0\left(x\right),v^0\left(y\right);{\alpha }^0\right]}^{†}$. Analogously, it is convenient to consider the variations in the model’s state functions and parameters to be the components of a “vector of variations”, $\delta e$, defined as follows: $\delta e\triangleq {\left[\delta u\left(x\right),\delta v\left(y\right);\delta \alpha \right]}^{†}$. The 1st-order Gateaux- (G-) variation of the response $R\left(e\right)$, which will be denoted as $\delta R\left(e^0;\delta e\right)$, for arbitrary variations $\delta e$ in the model parameters and state functions in a neighborhood $\left(e^0+\epsilon \delta e\right)$ around $e^0$, is obtained, by definition, as follows:

$$\delta R\left(e^0;\delta e\right)\triangleq {\left\{\frac{d}{d\epsilon }R\left[u^0\left(x\right)+\epsilon \delta u\left(x\right),v^0\left(y\right)+\epsilon \delta v\left(y\right);{\alpha }^0+\epsilon \delta \alpha ;x,y\right]\right\}}_{\epsilon =0}$$

(8)

The unknown variations $\delta u\left(x\right)$ and $\delta v\left(y\right)$ in the state functions are related to the variations $\delta \alpha$ through the equations obtained by applying the definition of the G-differential to the equations underlying the coupled nonlinear, i.e., Equations (1)–(3), to obtain the following relations:

$${\left\{\frac{d}{d\epsilon }{N}^{\left(I\right)}\left[u^0\left(x\right)+\epsilon \delta u\left(x\right);{\alpha }^0+\epsilon \delta \alpha \right]\right\}}_{\epsilon =0}=\hspace{0.17em}{\left\{\frac{d}{d\epsilon }{Q}^{\left(I\right)}\left({\alpha }^0+\epsilon \delta \alpha ;x\right)\right\}}_{\epsilon =0},x\in {\mathsf{\Omega }}_x\left({\alpha }^0\right),$$

(9)

$${\left\{\frac{d}{d\epsilon }{N}^{\left(II\right)}\left[v^0\left(y\right)+\epsilon \delta v\left(y\right);{\alpha }^0+\epsilon \delta \alpha \right]\right\}}_{\epsilon =0}=\hspace{0.17em}{\left\{\frac{d}{d\epsilon }{Q}^{\left(II\right)}\left({\alpha }^0+\epsilon \delta \alpha ;y\right)\right\}}_{\epsilon =0},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in {\mathsf{\Omega }}_y\left({\alpha }^0\right),$$

(10)

$${\left\{\frac{d}{d\epsilon }B\left[u^0\left(x\right)+\epsilon \delta u\left(x\right),v^0\left(y\right)+\epsilon \delta v\left(y\right);{\alpha }^0+\epsilon \delta \alpha ;x,y\right]\right\}}_{\epsilon =0}=\hspace{0.17em}0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left({\alpha }^0\right).$$

(11)

Performing in Equations (9)–(11) the differentiations with respect to $\epsilon$ and setting $\epsilon =0$ in the resulting expression yields the following system of equations:

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

(12)

$$\delta N^{\left(II\right)}\left[v^0\left(y\right),{\alpha }^0;\delta v\left(y\right),\delta \alpha \right]=\delta Q^{\left(II\right)}\left({\alpha }^0;\delta \alpha \right)\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in {\mathsf{\Omega }}_y\left({\alpha }^0\right),$$

(13)

$$\delta B\left(e^0;\delta e\right)=\hspace{0.17em}0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left({\alpha }^0\right)$$

(14)

The system of equations comprising Equations (12)–(14) is called the “First-Level Forward Sensitivity System” (1st-LFSS) and could be solved to obtain the variations $\delta u\left(x\right)$ and $\delta v\left(y\right)$ in the state functions in terms of the parameter variations $\delta \alpha$ which appear as sources in the 1st-LFSS equations. Subsequently, the variations $\delta u\left(x\right)$ and $\delta v\left(y\right)$ thus obtained could be used to compute the total sensitivity $\delta R\left(e^0;\delta e\right)$ defined in Equation (8).

The existence of the G-variations of the operators underlying the 1st-LFSS and total sensitivity $\delta R\left(e^0;\delta e\right)$ does not guarantee their numerical computability. Numerical methods most often require that $\delta R\left(e^0;\delta e\right)$ and the operators underlying the 1st-LFSS be linear in the variations $\delta e$ in a neighborhood $\left(e^0+\epsilon \delta e\right)$ around $e^0$. The necessary and sufficient conditions for the G-differential $\delta W\left(e^0;\delta e\right)$ of a nonlinear operator $W\left(e\right)$ to be linear in $\delta e$ in a neighborhood $\left(e^0+\epsilon \delta e\right)$ around $e^0$, and thus admit partial and total G-derivatives, are as follows [6]:

(i)

$W\left(e\right)$ satisfies a weak Lipschitz condition at $e^0$;

$$\Vert W\left(e^0+\epsilon h;x\right)-W\left(e^0;x\right)\Vert \le k\Vert \epsilon e^0\Vert ,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k<\infty$$

(15)

(ii)

for two arbitrary vectors of variations $\hspace{0.17em}\delta e_1$ and $\delta e_2$, the operator $W\left(e\right)$ satisfies the following relation:

$$W\left(e^0+\epsilon \delta e_1+\epsilon \delta e_2\right)-W\left(e^0+\epsilon \delta e_1\right)-W\left(e^0+\epsilon \delta e_2\right)+W\left(e^0\right)=o\left(\epsilon \right)\hspace{0.17em}$$

(16)

It will henceforth be assumed that the operators $N^{\left(I\right)}$, $N^{\left(II\right)}$, $B$, $Q^{\left(I\right)}$, $Q^{\left(II\right)}$ and $R$ satisfy the conditions indicated in Equations (15) and (16). Hence, Equations (12)–(14) can be written in the following form:

$${\left\{\hspace{0.17em}\frac{\partial N^{\left(I\right)}\left(u;\alpha \right)}{\partial u}\right\}}_{\left(u^0;{\alpha }^0\right)}\delta u\left(x\right)=\hspace{0.17em}{\left\{Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\hspace{0.17em}\right\}}_{\left(u^0;{\alpha }^0\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in {\mathsf{\Omega }}_x\left({\alpha }^0\right),$$

(17)

$${\left\{\frac{\partial N^{\left(II\right)}\left(v;\alpha \right)}{\partial v}\right\}}_{\left(v^0,{\alpha }^0\right)}\delta v\left(y\right)=\hspace{0.17em}{\left\{Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\hspace{0.17em}\right\}}_{\left(v^0,{\alpha }^0\right)},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in {\mathsf{\Omega }}_y\left({\alpha }^0\right),$$

(18)

$$\begin{array}{l}{\left\{\frac{\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]}{\partial u}\right\}}_{\left(e^0\right)}\delta u\left(x\right)+{\left\{\frac{\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]}{\partial v}\right\}}_{\left(e^0\right)}\delta v\left(y\right)\hspace{0.17em}\\ +{\left\{\frac{\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]}{\partial \alpha }\right\}}_{\left(e^0\right)}\delta \alpha =\hspace{0.17em}0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left({\alpha }^0\right).\end{array}$$

(19)

where

$${\left\{Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\hspace{0.17em}\right\}}_{\left(u^0;{\alpha }^0\right)}\triangleq {\left\{\frac{\partial \left[Q^{\left(I\right)}\left(\alpha ;x\right)-N^{\left(I\right)}\left[u\left(x\right);\alpha \right]\right]}{\partial \alpha }\hspace{0.17em}\right\}}_{\left(u^0;{\alpha }^0\right)}\delta \alpha ,$$

(20)

$$\hspace{0.17em}{\left\{Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\hspace{0.17em}\right\}}_{\left(v^0,{\alpha }^0\right)}\triangleq {\left\{\frac{\partial \left[Q^{\left(II\right)}\left(\alpha ;y\right)-N^{\left(II\right)}\left[v\left(y\right);\alpha \right]\right]}{\partial \alpha }\hspace{0.17em}\right\}}_{\left(v^0,{\alpha }^0\right)}\delta \alpha ,$$

(21)

The partial G-derivatives $\partial N^{\left(I\right)}\left[u\left(x\right);\alpha \right]/\partial u$, $\partial N^{\left(II\right)}\left[v\left(y\right);\alpha \right]/\partial v$, $\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]/\partial u$, $\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]/\partial v$, $\partial N^{\left(I\right)}\left[u\left(x\right);\alpha \right]/\partial \alpha$, $\partial N^{\left(II\right)}\left[v\left(y\right);\alpha \right]/\partial \alpha$, $\partial Q^{\left(I\right)}\left[u\left(x\right);\alpha \right]/\partial \alpha$, $\partial Q^{\left(II\right)}\left[u\left(x\right);\alpha \right]/\partial \alpha$ and $\partial B\left[u\left(x\right),v\left(y\right);\alpha ;x,y\right]/\partial \alpha$, which appear in Equations (17)–(21), are matrices of corresponding dimensions. When the G-variation $\delta R\left(e^0;\delta e\right)$ is linear in $\delta e$, it is called the G-differential of $R\left(e\right)$ and is usually denoted as $DR\left(e^0;\delta e\right)$. Furthermore, the result of the differentiations indicated on the right-side of the definition provided in Equation (8) can be written as follows:

$$DR\left(e^0;\delta e\right)={\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{dir}+{\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{ind},$$

(22)

where the so-called “direct-effect” term is defined as follows:

$${\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{dir}\triangleq {\left\{\frac{\partial R}{\partial \alpha }\right\}}_{\left(e^0\right)}\delta \alpha ={\displaystyle \sum _{i=1}^{Z_{\alpha }}{\left\{\frac{\partial R}{\partial {\alpha }_i}\right\}}_{\left(u^0;{\alpha }^0\right)}\delta {\alpha }_i}\hspace{0.17em},$$

(23)

while the so-called “indirect-effect” term is defined as follows:

$${\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{ind}\triangleq {\left\{\frac{\partial R}{\partial u}\right\}}_{\left(e^0\right)}\delta u\left(x\right)+{\left\{\frac{\partial R}{\partial v}\right\}}_{\left(e^0\right)}\delta v\left(y\right)\hspace{0.17em}.$$

(24)

In Equations (23) and (24), the vectors $\partial R/\partial u$, $\partial R/\partial v$ and $\partial R/\partial \alpha$ comprise, as components, the first-order partial G-derivatives computed at the phase–space point $e^0$. The G-differential $DR\left(e^0;\delta e\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-differential $DR\left(e^0;\delta e\right)$ satisfies the relation $R\left(e^0+\epsilon \delta e\right)-R\left(e^0\right)=DR\left(e^0;\delta e\right)+\Delta \left(\delta e\right)\hspace{0.17em},$ with $\underset{\epsilon \to 0}{\lim }\left[\Delta \left(\epsilon \delta e\right)\right]/\epsilon =0$.

The “direct effect” term ${\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{dir}$ depends only on the parameter variations $\delta \alpha$ so it can be computed immediately, since it does not depend on the variations $\delta u$ and $\delta v$. On the other hand, “indirect effect” term ${\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{ind}$ depends indirectly on the parameter variations $\delta \alpha$ through the yet unknown variations $\delta u\left(x\right)$ and $\delta v\left(y\right)$ in the state functions, and these functions can be determined only by solving repeatedly the 1st-LFSS for every possible parameter variation $\delta {\alpha }_i,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}i=1,\dots ,Z_{\alpha }$. The need for these prohibitively expensive computations can be circumvented by extending the concepts underlying the “Adjoint Sensitivity Analysis Methodology” (ASAM) conceived by Cacuci [1] to construct a “First-Level Adjoint Sensitivity System” (1st-LASS), the solution of which will be independent of the variations $\delta \alpha$, $\delta u\left(x\right)$ and $\delta v\left(y\right)$. Subsequently, the solution of the 1st-LASS will be used to compute the indirect-effect term ${\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{ind}$ by constructing an equivalent expression (for this indirect-effect term) which would not involve the unknown variations $\delta u\left(x\right)$ and $\delta v\left(y\right)$.

2.1. Spectral Representation of the System Response’s Indirect-Effect Term

Since the indirect-effect term ${\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{ind}$ is defined on the same domain ${\mathsf{\Omega }}_x\left({\alpha }^0\right)\cup {\mathsf{\Omega }}_y\left({\alpha }^0\right)$ as $R\left(e^0\right)$, and has the same range as $R\left(e^0\right)$, it follows that it can be represented in the following form:

$$\begin{array}{l}{\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{ind}=\\ {\displaystyle \sum _{m_1=0}^{\infty }\dots }{\displaystyle \sum _{m_{Z_x}=0}^{\infty }{\displaystyle \sum _{n_1=0}^{\infty }\dots }{\displaystyle \sum _{n_{Z_y}=0}^{\infty }{\left\{F_{m_1\dots m_{Z_x}{n}_1\dots n_{Z_y}}\left(u;v;\alpha ;\delta u\right)+G_{m_1\dots m_{Z_x}{n}_1\dots n_{Z_y}}\left(u;v;\alpha ;\delta v\right)\right\}}_{\left(e^0\right)}}}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \hspace{0.17em}{P}_{m_1}\left(x_1\right)\dots P_{m_{Z_x}}\left(x_{Z_x}\right)O_{n_1}\left(y_1\right)\dots O_{n_{Z_y}}\left(y_{Z_y}\right).\end{array}$$

(25)

where

$$\begin{array}{l}{\left\{F_{m_1\dots m_{Z_x}{n}_1\dots n_{Z_y}}\left(u;v;\alpha ;\delta u\right)\right\}}_{\left(e^0\right)}\triangleq {\displaystyle \underset{a_1\left({\alpha }^0\right)}{\overset{b_1\left({\alpha }^0\right)}{\int }}d{x}_1\hspace{0.17em}\dots {\displaystyle \underset{a_{Z_x}\left({\alpha }^0\right)}{\overset{b_{Z_x}\left({\alpha }^0\right)}{\int }}d{x}_{Z_x}{\displaystyle \underset{c_1\left({\alpha }^0\right)}{\overset{d_1\left({\alpha }^0\right)}{\int }}d{y}_1\dots {\displaystyle \underset{c_{Z_y}\left({\alpha }^0\right)}{\overset{d_{Z_y}\left({\alpha }^0\right)}{\int }}d{y}_{Z_y}{\left\{\frac{\partial R\left[u\left(x\right);v\left(y\right);\alpha ;x;y\right]}{\partial u}\delta u\right\}}_{\left(e^0\right)}\hspace{0.17em}}}}} \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \hspace{0.17em}{P}_{m_1}\left(x_1\right)\dots P_{m_{Z_x}}\left(x_{Z_x}\right)O_{n_1}\left(y_1\right)\dots O_{n_{Z_y}}\left(y_{Z_y}\right)\end{array}$$

(26)

and

$$\begin{array}{l}{\left\{G_{m_1\dots m_{Z_x}{n}_1\dots n_{Z_y}}\left(u;v;\alpha ;\delta v\right)\right\}}_{\left(e^0\right)}\triangleq {\displaystyle \underset{a_1\left({\alpha }^0\right)}{\overset{b_1\left({\alpha }^0\right)}{\int }}d{x}_1\hspace{0.17em}\dots {\displaystyle \underset{a_{Z_x}\left({\alpha }^0\right)}{\overset{b_{Z_x}\left({\alpha }^0\right)}{\int }}d{x}_{Z_x}{\displaystyle \underset{c_1\left({\alpha }^0\right)}{\overset{d_1\left({\alpha }^0\right)}{\int }}d{y}_1\dots {\displaystyle \underset{c_{Z_y}\left({\alpha }^0\right)}{\overset{d_{Z_y}\left({\alpha }^0\right)}{\int }}d{y}_{Z_y}{\left\{\frac{\partial R\left[u\left(x\right);v\left(y\right);\alpha ;x;y\right]}{\partial v}\delta v\right\}}_{\left(e^0\right)}\hspace{0.17em}}}}} \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \hspace{0.17em}{P}_{m_1}\left(x_1\right)\dots P_{m_{Z_x}}\left(x_{Z_x}\right)O_{n_1}\left(y_1\right)\dots O_{n_{Z_y}}\left(y_{Z_y}\right).\end{array}$$

(27)

The following designations have been used in Equations (26) and (27): (i) the quantities $P_{m_i}\left(x_i\right)$, $i=1,\dots ,Z_x$, denote the corresponding spectral basis functions (e.g., orthogonal polynomials, Fourier exponential/trigonometric functions) for the domain defined as the domain ${\mathsf{\Omega }}_x$; (ii) the quantities $O_{m_i}\left(y_i\right)$, $i=1,\dots ,Z_y$, denote the spectral functions corresponding to the domain ${\mathsf{\Omega }}_y$; and (iii) the quantities ${\left\{F_{m_1\dots m_{Z_x}{n}_1\dots n_{Z_y}}\left(u;v;\alpha ;\delta u\right)\right\}}_{\left(e^0\right)}$ and ${\left\{G_{m_1\dots m_{Z_x}{n}_1\dots n_{Z_y}}\left(u;v;\alpha ;\delta v\right)\right\}}_{\left(e^0\right)}$ denote the corresponding generalized spectral (Fourier) coefficients.

The appearance of the “difficult to compute” variations $\delta u$ and $\delta v$ in the functionals defined in Equations (26) and (27), respectively, can be eliminated by expressing the right-sides of Equations (26) and (27) in terms of adjoint functions that will be obtained by implementing the following sequence of steps:

• Introduce a Hilbert space pertaining to the domain ${\mathsf{\Omega }}_x\left({\alpha }^0\right)\cup {\mathsf{\Omega }}_y\left({\alpha }^0\right)$, denoted as $\mathrm{H}$, comprising square-integrable vector-valued elements of the form $f^{\left(\alpha \right)}\left(x,y\right)\triangleq {\left[g^{\left(\alpha \right)}\left(x,y\right),h^{\left(\alpha \right)}\left(x,y\right)\right]}^{†}\in \mathrm{H}$ and $f^{\left(\beta \right)}\left(x,y\right)\triangleq {\left[g^{\left(\beta \right)}\left(x,y\right),h^{\left(\beta \right)}\left(x,y\right)\right]}^{†}\in \mathrm{H}$, where $g^{\left(\alpha \right)}\left(x,y\right)\triangleq {\left[g_1^{\left(\alpha \right)}\left(x,y\right),\dots ,g_{Z_u}^{\left(\alpha \right)}\left(x,y\right)\right]}^{†}$, $g^{\left(\beta \right)}\left(x,y\right)\triangleq {\left[g_1^{\left(\beta \right)}\left(x,y\right),\dots ,g_{Z_u}^{\left(\beta \right)}\left(x,y\right)\right]}^{†}$, $h^{\left(\alpha \right)}\left(x,y\right)\triangleq {\left[h_1^{\left(\alpha \right)}\left(x,y\right),\dots ,h_{Z_v}^{\left(\alpha \right)}\left(x,y\right)\right]}^{†}$, $h^{\left(\beta \right)}\left(x,y\right)\triangleq {\left[h_1^{\left(\beta \right)}\left(x,y\right),\dots ,h_{Z_v}^{\left(\beta \right)}\left(x,y\right)\right]}^{†}$.
• Define the inner product, denoted as $\langle f^{\left(\alpha \right)}\left(x,y\right),f^{\left(\beta \right)}\left(x,y\right)\rangle$, between two elements of $\mathrm{H}$, as follows:

  $$\langle f^{\left(\alpha \right)}\left(x,y\right),f^{\left(\beta \right)}\left(x,y\right)\rangle \triangleq {\displaystyle \underset{a_1\left({\alpha }^0\right)}{\overset{b_1\left({\alpha }^0\right)}{\int }}\dots {\displaystyle \underset{a_{Z_x}\left({\alpha }^0\right)}{\overset{b_{Z_x}\left({\alpha }^0\right)}{\int }}{\displaystyle \underset{c_1\left({\alpha }^0\right)}{\overset{d_1\left({\alpha }^0\right)}{\int }}\dots {\displaystyle \underset{c_{Z_x}\left({\alpha }^0\right)}{\overset{d_{Z_x}\left({\alpha }^0\right)}{\int }}\left[g^{\left(\alpha \right)}\left(x,y\right)·g^{\left(\beta \right)}\left(x,y\right)+h^{\left(\alpha \right)}\left(x,y\right)·h^{\left(\beta \right)}\left(x,y\right)\right]}}}}dxdy$$

  (28)

  where

  $$g^{\left(\alpha \right)}\left(x,y\right)·g^{\left(\beta \right)}\left(x,y\right)\triangleq {\displaystyle \sum _{n=1}^{Z_u}{g}_n^{\left(\alpha \right)}\left(x,y\right)g_n^{\left(\beta \right)}\left(x,y\right)}$$

  (29)

  and

  $$h^{\left(\alpha \right)}\left(x,y\right)·h^{\left(\beta \right)}\left(x,y\right)\triangleq {\displaystyle \sum _{n=1}^{Z_v}{h}_n^{\left(\alpha \right)}\left(x,y\right)h_n^{\left(\beta \right)}\left(x,y\right)}$$

  (30)
• Recast Equations (17) and (18) in the following matrix from:

  $${\left\{\left[\begin{array}{cc}\frac{\partial N^{\left(I\right)}\left(u;\alpha \right)}{\partial u}& 0\\ 0& \frac{\partial N^{\left(II\right)}\left(v;\alpha \right)}{\partial v}\end{array}\right]\hspace{0.17em}\right\}}_{\left(e^0\right)}\left(\begin{array}{c}\delta u\left(x\right)\\ \delta v\left(y\right)\end{array}\right)={\left\{\left(\begin{array}{c}{Q}_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\\ Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\end{array}\hspace{0.17em}\right)\hspace{0.17em}\right\}}_{\left(e^0\right)}$$

  (31)
• Use the definition provided in Equation (28) to form the inner product of Equation (31) with a square-integrable vector ${\mathsf{\psi }}^{\left(1\right)}\left(x,y\right)\triangleq {\left[{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right),{\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)\right]}^{†}\in \mathrm{H}$ to obtain the following relation:

  $$\begin{array}{l}\langle {\left(\begin{array}{c}{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right)\\ {\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left[\begin{array}{cc}\frac{\partial N^{\left(I\right)}\left(u;\alpha \right)}{\partial u}& 0\\ 0& \frac{\partial N^{\left(II\right)}\left(v;\alpha \right)}{\partial v}\end{array}\right]\hspace{0.17em}\right\}}_{\left(e^0\right)}\left(\begin{array}{c}\delta u\left(x\right)\\ \delta v\left(y\right)\end{array}\right)\hspace{0.17em}\rangle \\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\langle {\left(\begin{array}{c}{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right)\\ {\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left(\begin{array}{c}{Q}_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\\ Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\end{array}\hspace{0.17em}\right)\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\rangle .\end{array}$$

  (32)
• Using the definition of the adjoint operator in the Hilbert space $\mathrm{H}$, recast the left-side of Equation (32) as follows:

  $$\begin{array}{l}\langle {\left(\begin{array}{c}{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right)\\ {\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left[\begin{array}{cc}\frac{\partial N^{\left(I\right)}\left(u;\alpha \right)}{\partial u}& 0\\ 0& \frac{\partial N^{\left(II\right)}\left(v;\alpha \right)}{\partial v}\end{array}\right]\hspace{0.17em}\right\}}_{\left(e^0\right)}\left(\begin{array}{c}\delta u\left(x\right)\\ \delta v\left(y\right)\end{array}\right)\hspace{0.17em}\rangle \\ =\langle {\left(\begin{array}{c}\delta u\left(x\right)\\ \delta v\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}{\left\{\left[\begin{array}{cc}{A}^{*}\left(u;\alpha \right)& 0\\ 0& B^{*}\left(v;\alpha \right)\end{array}\right]\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\left(\begin{array}{c}{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right)\\ {\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)\end{array}\right)\hspace{0.17em}\rangle \\ +{\left\{B{C}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right),{\mathsf{\psi }}^{\left(II\right)}\left(x,y\right);\delta u\left(x\right),\delta v\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)},\end{array}$$

  (33)

  where the operator $A^{*}\left(u;\alpha \right)$ denotes the formal adjoint of $\partial N^{\left(I\right)}\left(u;\alpha \right)/\partial u$, the operator $B^{*}\left(v;\alpha \right)$ denotes the formal adjoint of $\partial N^{\left(II\right)}\left(v;\alpha \right)/\partial v$ and where ${\left\{B{C}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right),{\mathsf{\psi }}^{\left(II\right)}\left(x,y\right);\delta u\left(x\right),\delta v\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)}$ denotes the bilinear concomitant evaluated on the boundary $\delta {\mathsf{\Omega }}_x\left({\alpha }^0\right)\cup \delta {\mathsf{\Omega }}_y\left({\alpha }^0\right)$. The superscript “1” which appears in the notation of the bilinear concomitant $B{C}^{\left(1\right)}$ indicates that this quantity arises in conjunction with the construction of the “First-Level Adjoint Sensitivity System (1st-LASS)”.
• Replace the left-side of Equation (33) with the right-side of Equation (32) to obtain the following relation:

  $$\begin{array}{l}\langle {\left(\begin{array}{c}\delta u\left(x\right)\\ \delta v\left(y\right)\end{array}\right)}^{†},\hspace{0.17em}{\left\{\left[\begin{array}{cc}{A}^{*}\left(u;\alpha \right)& 0\\ 0& B^{*}\left(v;\alpha \right)\end{array}\right]\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\left(\begin{array}{c}{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right)\\ {\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)\end{array}\right)\hspace{0.17em}\rangle \\ =\langle {\left(\begin{array}{c}{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right)\\ {\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)\end{array}\right)}^{†},\hspace{0.17em}\hspace{0.17em}{\left\{\left(\begin{array}{c}{Q}_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)\\ Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\end{array}\hspace{0.17em}\right)\hspace{0.17em}\right\}}_{\left(e^0\right)}\hspace{0.17em}\rangle \\ -{\left\{B{C}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right),{\mathsf{\psi }}^{\left(II\right)}\left(x,y\right);\delta u\left(x\right),\delta v\left(y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)}.\end{array}$$

  (34)
• Require the left-side of Equation (34) to represent the indirect-effect term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}_{indirect}$ defined in Equation (25), which can be fulfilled by requiring the yet undetermined (adjoint) functions ${\mathsf{\psi }}^{\left(I\right)}\left(x,y\right)$ and ${\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)$ to satisfy the following equations:

  $${\left\{A^{*}\left(u;\alpha \right)\hspace{0.17em}\right\}}_{\left(e^0\right)}{\mathsf{\psi }}^{\left(I\right)}\left(x,y\right)={\left\{\frac{\partial R\left[u\left(x\right);v\left(y\right);\alpha ;x;y\right]}{\partial u}\delta u\right\}}_{\left(e^0\right)}{P}_{m_1}\left(x_1\right)\dots P_{m_{Z_x}}\left(x_{Z_x}\right)O_{n_1}\left(y_1\right)\dots O_{n_{Z_y}}\left(y_{Z_y}\right),$$

  (35)

  $${\left\{B^{*}\left(v;\alpha \right)\hspace{0.17em}\right\}}_{\left(e^0\right)}{\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)={\left\{\frac{\partial R\left[u\left(x\right);v\left(y\right);\alpha ;x;y\right]}{\partial v}\delta v\right\}}_{\left(e^0\right)}{P}_{m_1}\left(x_1\right)\dots P_{m_{Z_x}}\left(x_{Z_x}\right)O_{n_1}\left(y_1\right)\dots O_{n_{Z_y}}\left(y_{Z_y}\right).$$

  (36)
• Since the source terms on the right-sides of Equations (35) and (36) depend on the indices of the spectral bases functions, it follows that the adjoint functions ${\mathsf{\psi }}^{\left(I\right)}\left(x,y\right)$ and ${\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)$ also depend on the respective indices, which will henceforth be explicitly displayed by writing ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}\left(x,y\right)$ and ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}\left(x,y\right)$, respectively.
• The boundary, interface, and initial/final conditions for the functions ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}\left(x,y\right)$ and ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}\left(x,y\right)$ are now determined by imposing the following requirements:

  (a)

  Implement the boundary, interface and initial/final conditions given in Equation (19) into the bilinear concomitant in Equation (34).

  (b)

  Eliminate the remaining unknown boundary, interface and initial/final conditions involving the functions $\delta u\left(x\right)$ and $\delta v\left(y\right)$ from the expression of the bilinear concomitant in Equation (34) by selecting boundary, interface and initial/final conditions for the adjoint functions ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}\left(x,y\right)$ and ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}\left(x,y\right)$ such that the selected conditions for these adjoint functions must be independent of unknown values of $\delta u\left(x\right)$, $\delta v\left(y\right)$ and $\delta \alpha$ while ensuring that Equations (35) and (36) are well posed. The boundary conditions thus chosen for the adjoint functions ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}\left(x,y\right)$ and ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}\left(x,y\right)$ can be represented in operator form as follows:

  $$\begin{array}{l}{\left\{B_A^{\left(1\right)}\left[u\left(x\right);v\left(y\right);{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}\left(x,y\right),{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}\left(x,y\right);\alpha ;x,y\right]\right\}}_{\left(e^0\right)}=0,\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left({\alpha }^0\right),\end{array}$$

  (37)

  where the subscript “A” indicates “adjoint” and the superscript “1” indicates that these boundary conditions arises in conjunction with the construction of 1st-LASS. The selection of the boundary conditions for the adjoint functions ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}\left(x,y\right)$ and ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}\left(x,y\right)$ represented by Equation (37) eliminates the appearance of any unknown values of the variations $\delta u\left(x\right)$ and $\delta v\left(y\right)$ in the bilinear concomitant in Equation (34), reducing it to a residual quantity that contains boundary terms involving only known values of $\delta \alpha$, $u\left(x\right)$, $v\left(y\right)$, ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}\left(x,y\right)$, ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}\left(x,y\right)$, and $\alpha$. This residual bilinear concomitant will be denoted as ${\left\{R{C}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}\left(x,y\right),{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}\left(x,y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\delta \hspace{0.17em}{\mathsf{\Omega }}_x\cup }}_{\delta \hspace{0.17em}{\mathsf{\Omega }}_y}}\right\}}_{\left(e^0\right)}$.

  In general, this residual bilinear concomitant does not automatically vanish, although it may do so in particular instances. In principle, this residual bilinear concomitant could be forced to vanish, if necessary, by considering extensions, in the operator sense, of the linear operators $A^{*}\left(u;\alpha \right)\hspace{0.17em}$ and/or $B^{*}\left(v;\alpha \right)$, but such extensions seldom need to be used in practice.
• Using Equations (34)–(36) in conjunction with Equations (26) and (27) in Equation (25) yields the following expression for the indirect-effect term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{ind}$:

  $$\begin{array}{l}{\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{ind}\\ ={\displaystyle \sum _{m_1=0}^{\infty }\dots }{\displaystyle \sum _{m_{Z_x}=0}^{\infty }{\displaystyle \sum _{m_1=0}^{\infty }\dots }{\displaystyle \sum _{m_{Z_y}=0}^{\infty }\{\langle {\mathsf{\psi }}^{\left(I\right)}\left(x,y\right)Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)+{\mathsf{\psi }}^{\left(II\right)}\left(x,y\right)Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\rangle }}\\ {-R{C}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}\left(x,y\right),{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}\left(x,y\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\}}_{\left(e^0\right)}\\ \times \hspace{0.17em}{P}_{m_1}\left(x_1\right)\dots P_{m_{Z_x}}\left(x_{Z_x}\right)O_{n_1}\left(y_1\right)\dots O_{n_{Z_y}}\left(y_{Z_y}\right)\\ \triangleq {\left\{DR\left(e^0;{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)},{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)};\delta \alpha \right)\right\}}^{ind}.\end{array}$$

  (38)

As the expression in Equation (38) indicates, the desired elimination from ${\left\{DR\left(e^0;\delta \alpha \right)\right\}}^{ind}$ of the unknown variations $\delta u$ and $\delta v$ has been accomplished by having replaced them by the adjoint functions ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}$ and ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}$, which do not depend on any parameter variations; this fact that has been underscored by having explicitly indicated that the indirect-effect term can now be written in the form ${\left\{DR\left(e^0;{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)},{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)};\delta \alpha \right)\right\}}^{ind}$.

When first introduced in Equation (32), it was not known that the adjoint functions would ultimately depend on the indices $m_1,\dots ,m_{Z_x}$ and $n_1,\dots ,n_{Z_y}$; this fact has become apparent only after having constructed the right-sides (i.e., sources) of Equations (35) and(36) to emphasize this fact. These equations are re-written below:

$$\begin{array}{c}{\left\{A^{*}\left(u;\alpha \right)\hspace{0.17em}\right\}}_{\left(e^0\right)}{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}\left(x,y\right)={\left\{\frac{\partial R\left[u\left(x\right);v\left(y\right);\alpha ;x;y\right]}{\partial u}\delta u\right\}}_{\left(e^0\right)}\\ \times P_{m_1}\left(x_1\right)\dots P_{m_{Z_x}}\left(x_{Z_x}\right)O_{n_1}\left(y_1\right)\dots O_{n_{Z_y}}\left(y_{Z_y}\right),\end{array}$$

(39)

$$\begin{array}{c}{\left\{B^{*}\left(v;\alpha \right)\hspace{0.17em}\right\}}_{\left(e^0\right)}{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}\left(x,y\right)={\left\{\frac{\partial R\left[u\left(x\right);v\left(y\right);\alpha ;x;y\right]}{\partial v}\delta v\right\}}_{\left(e^0\right)}\\ \times P_{m_1}\left(x_1\right)\dots P_{m_{Z_x}}\left(x_{Z_x}\right)O_{n_1}\left(y_1\right)\dots O_{n_{Z_y}}\left(y_{Z_y}\right).\end{array}$$

(40)

The system of Equations (39) and (40), together with the adjoint boundary/initial conditions represented by Equation (37) will be called the “First-Level Adjoint Sensitivity System (1st-LASS).” The 1st-LASS is independent of the parameter variations $\delta \alpha$ but depends on the indices $m_1,\dots ,m_{Z_x}$ and $n_1,\dots ,n_{Z_y}$. In principle, therefore, the 1st-LASS needs to be solved as many times as there are nonzero spectral basis functions, which act as sources on the right side of the equations underlying the 1st-LASS. It is therefore very important to represent the indirect-effect term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{ind}$ defined in Equation (25) using as few basis-functions as possible, within a criterion of accuracy that is set by the user, a priori. Once the adjoint functions ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)}$ and ${\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)}$ are available, they can be used in Equation (38) to compute the indirect-effect term ${\left\{DR\left(e^0;{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(I\right)},{\mathsf{\psi }}_{m_1..m_{Z_x}{n}_1\dots n_{Z_y}}^{\left(II\right)};\delta \alpha \right)\right\}}^{ind}$ exactly and efficiently, using quadrature formulas, which are many orders of magnitude faster to compute than solving the operator (differential, integral) equations that underlie the 1st-LFSS.

In practice, orthogonal polynomials will often be selected to serve as basis-functions for the spectral Fourier representations of the responses of interest. As is well-known, orthogonal polynomials possess many recurrence relations which can be advantageously used to reduce massively the number of computations that would actually require solving the 1st-LASS.

In the particular case when the response is a scalar-valued functional of the system’s dependent variables, the expansion in Equation (25) reduces to a single term, so that the summations in the expression of the indirect-effect term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{ind}$ in Equation (38) also reduce to a single term.

2.2. Pseudo-Spectral Representation of the System Response’s Indirect-Effect Term

Alternatively, Lagrange interpolation, see e.g., [7], can be used to express the indirect-effect term defined in Equation (24) approximately as follows:

$${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{ind}\cong {\displaystyle \sum _{i,j=0}^N{R}_{ij}\left(x_i,y_j\right)C{F}_{ij}\left(x,y\right)}.$$

(41)

where the quantities $C{F}_{ij}\left(x,y\right)$ represent the “cardinal functions”, where $x_i\hspace{0.17em}$ and $y_j$ denote the collocation (or interpolation) points, and where

$$\begin{array}{l}{R}_{ij}\left(x_i,y_j\right)={\displaystyle \underset{a_1\left({\alpha }^0\right)}{\overset{b_1\left({\alpha }^0\right)}{\int }}..{\displaystyle \underset{a_{Z_x}\left({\alpha }^0\right)}{\overset{b_{Z_x}\left({\alpha }^0\right)}{\int }}{\displaystyle \underset{c_1\left({\alpha }^0\right)}{\overset{d_1\left({\alpha }^0\right)}{\int }}..{\displaystyle \underset{c_{Z_x}\left({\alpha }^0\right)}{\overset{d_{Z_x}\left({\alpha }^0\right)}{\int }}\left[{\left(\frac{\partial R}{\partial u}\right)}_{\left(e^0\right)}\delta u\left(x\right)+{\left(\frac{\partial R}{\partial v}\right)}_{\left(e^0\right)}\delta v\left(y\right)\right]}}}}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\times \delta \left(x\hspace{0.17em}-x_i\hspace{0.17em}\right)\delta \left(y-y_j\hspace{0.17em}\right)dxdy\hspace{0.17em}.\end{array}$$

(42)

The cardinal functions $C{F}_{ij}\left(x,y\right)$ are also called [3] the “fundamental polynomials for pointwise interpolation”, the “elements of the cardinal basis”, the “Lagrange basis”, or the “shape functions”. Depending on the domains of definition $\hspace{0.17em}x\in {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}y\in {\mathsf{\Omega }}_y\left({\alpha }^0\right)$ and choices of weight functions, particularly important cardinal functions are the Chebyshev, Legendre, Gegenbauer, Hermite, Laguerre polynomials, and Whittaker’s “sinc” function. In several dimensions, it is most efficient to use a tensor product basis, i.e., use basis functions that are products of one-dimensional basis functions. Particularly efficient computational procedures can be constructed when both the basis functions and grid are tensor products of one-dimensional functions and grids, respectively. Using trigonometric functions, Chebyshev polynomials, or rational Chebyshev functions as basis functions enables the use of the Fast Fourier Transform, which further enhances computational efficiency.

Following established practice [3], “collocation points” and “interpolation points” will be used as synonyms in this work, as will be the terms “collocation” and “pseudospectral” when referring to the fact that interpolatory methods will be used to determine the yet unknown indirect-effect term ${\left\{DR\left(e^0;\delta u,\delta v\right)\right\}}^{ind}$ by expressing it in terms of adjoint functions specifically developed for each of the collocation/interpolation points. The reason that “collocation” methods are alternatively labeled “pseudospectral” is that the optimum choice of the interpolation points makes collocation methods identical with the Galerkin method if the inner products are evaluated by “Gaussian integration”. It is important to note that neither the cardinal functions $C{F}_{ij}\left(x,y\right)$ nor the collocation points $x_i\hspace{0.17em}$ and $y_j$ are subject to model parameter uncertainties.

The functionals $F_{ij}\left(x_i,y_j\right)$ defined in Equation (42) can be evaluated by using adjoint functions that are the solutions of a 1st-LASS constructed by following the same conceptual steps as those leading to Equations (39) and (40), and the adjoint boundary conditions defined by Equation (37). Omitting these intermediate steps, the final result is as follows:

$$\begin{array}{l}{F}_{ij}\left(x_i,y_j\right)\\ =\{\langle {\mathsf{\psi }}^{\left(A\right)}\left(x,y;x_i,y_j\right)Q_1^{\left(1\right)}\left(u;\alpha ;\delta \alpha \right)+{\mathsf{\psi }}^{\left(B\right)}\left(x,y;x_i,y_j\right)Q_2^{\left(1\right)}\left(v;\alpha ;\delta \alpha \right)\rangle \\ {-{\widehat{C}}^{\left(1\right)}{\left[u\left(x\right);v\left(y\right);{\mathsf{\psi }}^{\left(A\right)}\left(x,y;x_i,y_j\right),{\mathsf{\psi }}^{\left(B\right)}\left(x,y;x_i,y_j\right);\alpha ;x,y;\delta \alpha \right]}_{{{}_{\partial {\mathsf{\Omega }}_x\cup }}_{\partial {\mathsf{\Omega }}_y}}\}}_{\left(e^0\right)}\end{array}$$

(43)

where the adjoint functions ${\mathsf{\psi }}^{\left(A\right)}\left(x,y;x_i,y_j\right)$ and ${\mathsf{\psi }}^{\left(B\right)}\left(x,y;x_i,y_j\right)$ are the solutions of the following 1st-LASS:

$${\left\{A^{*}\left(u;\alpha \right)\hspace{0.17em}\right\}}_{\left(e^0\right)}{\mathsf{\psi }}^{\left(A\right)}\left(x,y;x_i,y_j\right)={\left\{\frac{\partial R\left[u\left(x\right);v\left(y\right);\alpha ;x;y\right]}{\partial u}\delta u\right\}}_{\left(e^0\right)}\delta \left(x\hspace{0.17em}-x_i\hspace{0.17em}\right)\delta \left(y-y_j\hspace{0.17em}\right)$$

(44)

$${\left\{B^{*}\left(v;\alpha \right)\hspace{0.17em}\right\}}_{\left(e^0\right)}{\mathsf{\psi }}^{\left(B\right)}\left(x,y;x_i,y_j\right)={\left\{\frac{\partial R\left[u\left(x\right);v\left(y\right);\alpha ;x;y\right]}{\partial v}\delta v\right\}}_{\left(e^0\right)}\delta \left(x\hspace{0.17em}-x_i\hspace{0.17em}\right)\delta \left(y-y_j\hspace{0.17em}\right)$$

(45)

$$\begin{array}{l}{\left\{C_A^{\left(1\right)}\left[u\left(x\right);v\left(y\right);{\mathsf{\psi }}^{\left(A\right)}\left(x,y;x_i,y_j\right),{\mathsf{\psi }}^{\left(B\right)}\left(x,y;x_i,y_j\right);\alpha ;x,y\right]\right\}}_{\left(e^0\right)}=0,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}x\in \partial {\mathsf{\Omega }}_x\left({\alpha }^0\right),\hspace{0.17em}\hspace{0.17em}y\in \partial {\mathsf{\Omega }}_y\left({\alpha }^0\right),\end{array}$$

(46)

It is evident from Equations (44)–(46) that the 1st-LASS must be solved anew for each of the collocation/interpolation points considered in the expansion of the indirect-effect term shown in Equation (41). The choice between using the spectral expansion shown in Equation (25) or using the collocation/interpolation pseudo-spectral expansion shown in Equation (41) depends on the specific problem under consideration, but for comparable accuracy in the computation of the response sensitivities, using the collocation/interpolation pseudo-spectral expansion shown in Equation (41) is often more efficient computationally than using the full spectral expansion.

The practical implementation of the mathematical methodology underlying the 1st-CASAM is illustrated in Figure 1 and Figure 2. The derivation of the 1st-LFSS is illustrated in Figure 1. The path on the left side of Figure 1 depicts the derivation of the (non-discretized) 1st-LFSS starting from the differential equations underlying the original nonlinear system. On the other hand, the path on the right-side of Figure 1 depicts the derivation of the discretized 1st-LFSS starting from the discretized form of the original nonlinear equations. If this path is followed, it must be ensured that the discretized 1st-LFSS is consistent with the differential form of the 1st-LFSS in the limit of vanishing size of the discretization interval considered for the independent variables.

The derivation of the 1st-LASS is illustrated in Figure 2. The path on the left side of Figure 2 depicts the derivation of the (non-discretized) 1st-LASS starting from the differential form of the 1st-LFSS. On the other hand, the path on the right side of Figure 2 depicts the derivation of the discretized 1st-LASS starting from the discretized 1st-LFSS. If this path is chosen, the consistency of the discretized 1st-LASS with the differential form of the 1st-LFSS must again be ensured.

3. Concluding Remarks

This work has presented the First-Order Comprehensive Adjoint Sensitivity Analysis Methodology (1st-CASAM) for computing efficiently the exact first-order sensitivities (i.e., functional derivatives) of operator-valued responses (i.e., model results) of general models of coupled nonlinear physical systems characterized by imprecisely known parameters, internal interfaces between the coupled systems and external boundaries. When the model response is a (scalar-valued) functional of the system’s dependent variables (i.e., state functions), the total sensitivity of a scalar-valued functional response to all of the model’s state functions is (also) a functional of the variations in the model’s state variables. By being a functional of the variations in the model’s state variables, the total response sensitivity naturally defines an inner product in terms of which it can be expressed uniquely by virtue of the well-known Riesz Representation Theorem (which ensures that every functional defined in a Hilbert space can be expressed uniquely as an inner product). The existence of such a natural inner-product induced by a functional response enables the construction of an appropriate adjoint sensitivity system, the solution of which (i.e., the respective adjoint sensitivity functions) can always be used to compute, exactly and most efficiently, the sensitivities of a functional response to the model’s scalar parameters. When the response is a functional of the state variables, a single adjoint computation (i.e., solution of the adjoint sensitivity system) suffices for subsequently computing exactly all of the model’s response sensitivities to all of the model’s scalar parameters. The adjoint sensitivity system has the same dimensions as the original system, but it is always linear in the adjoint state functions. This is in contradistinction to the original system, which is usually nonlinear in its state functions. Solving the original forward system and the adjoint sensitivity system involve large-scale computations, since these systems invariably involve inversion of large matrices stemming from differential, difference, integral, and/or algebraic equations. Since the adjoint sensitivity analysis methodology requires solving just once the adjoint sensitivity system, this methodology is the most advantageous to use computationally in practice for large-scale systems involving many parameters.

On the other hand, the total sensitivity (to model parameters and state functions) of a model response which is a function-valued (as opposed to a scalar-valued) operator of the model’s state functions does not provide a natural inner product for the model/system under consideration. Without an inner product, it is not possible to construct an adjoint sensitivity system, the solution of which would subsequently be used for computing the response sensitivities to the model’s parameters. Therefore, an inner product must first be constructed to enable expressing the operator-valued total response sensitivity to the variations in the state functions in terms of functionals of the system’s dependent variables (state functions). The requisite inner product can be constructed by representing the total sensitivity of the operator-valued response to the system’s state functions in terms of scalar-valued (functionals) response using: (i) spectral expansions; or (ii) collocation/pseudo-spectral expansions, or (iii) combined spectral/collocation expansions. The coefficients in any of these expansions are functionals that can be represented in terms of an inner product. In turn, this inner product enables the construction of an adjoint sensitivity system, the solution of which can subsequently be used to compute exactly and efficiently the sensitivities of these coefficients to the model’s parameters. A different source for the adjoint sensitivity system is developed for each spectral coefficient or for each collocation point. Altogether, therefore, as many adjoint computations would be needed as there are spectral coefficients and/or collocation points in the phase–space of independent variables. Thus, for operator-valued responses, the fundamental issue is to establish the number of collocation points in the phase–space of independent variables and/or the number of Fourier coefficients which would be needed for representing the response within an a priori established accuracy in the phase–space of independent variables. Subsequently, for each Fourier coefficient and/or at each collocation point, the 1st-CASAM provides the exact sensitivities in the parameter space, in the computationally most efficient manner. By enabling the exact computations of operator-valued response sensitivities to internal interfaces and external boundary parameters and conditions, the 1st-CASAM presented in this work makes it possible, inter alia, to quantify the effects of manufacturing tolerances on the responses of physical and engineering systems.

An accompanying work [7] will present the application of the 1st-CASAM developed in this work to a benchmark problem [8] that models coupled heat conduction and convection in a physical system comprising an electrically heated rod surrounded by a coolant which simulates the geometry of a nuclear reactor. In particular, this benchmark [8] was used to verify [8,9] the numerical results produced by the FLUENT Adjoint Solver [10].