The nth-order features adjoint sensitivity analysis methodology for response-coupled forward/adjoint linear systems (nth-FASAM-L): I. mathematical framework

This work presents the mathematical/theoretical framework of the “nth-Order Feature Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems” (abbreviated as “nth-FASAM-L”), which enables the most efficient computation of exactly obtained mathematical expressions of arbitrarily-high-order (nth-order) sensitivities of a generic system response with respect to all of the parameters (including boundary and initial conditions) underlying the respective forward/adjoint systems. Responses of linear models can depend simultaneously on both the forward and the adjoint state functions. This is in contradistinction to responses for nonlinear systems, which can depend only on the forward state functions since nonlinear operators do not admit bona-fide adjoint operators. Among the best-known model responses that depend simultaneously on both the forward and adjoint state functions are Lagrangians used for system optimization, the Schwinger and Roussopoulos functionals for analyzing reaction rates and ratios thereof, and the Rayleigh quotient for analyzing eigenvalues and/or separation constants. The sensitivity analysis of such responses makes it necessary to treat linear models/systems in their own right, rather than treating them just as particular cases of nonlinear systems. The unparalleled efficiency and accuracy of the nth-FASAM-L methodology stems from the maximal reduction of the number of adjoint computations (which are “large-scale” computations) for computing high-order sensitivities, since the number of large-scale computations when applying the nth-FASAM-N methodology is proportional to the number of model features as opposed to the number of model parameters (which are considerably more than the number of features). The mathematical framework underlying the nth-FASAM-L is developed in linearly increasing higher-dimensional Hilbert spaces, as opposed to the exponentially increasing “parameter-dimensional” spaces in which response sensitivities are computed by other methods (statistical, finite differences, etc.), thus providing the basis for overcoming the curse of dimensionality in sensitivity analysis and all other fields (uncertainty quantification, predictive modeling, etc.) which need such sensitivities.

1 Introduction

The analysis of computational models fundamentally relies on the use of functional derivatives (called “sensitivities”) of the results (called “model responses”) with respect to the imprecisely known parameters underlying the computational model. Sensitivities are used for many purposes, including: (a) ranking the importance of the various parameters and performing “reduced-order modeling” by eliminating unimportant parameters and/or processes; (b) quantifying the uncertainties induced in a model response due to uncertainties in the model’s parameters; (c) performing “model validation” by comparing computational and experimental results to address the question “does the model represent reality?”; (d) performing data assimilation and model calibration as part of forward and inverse “predictive modeling” to obtain best-estimate predicted results with reduced predicted uncertainties; (e) prioritizing improvements while optimizing the model.

Response sensitivities are computed by using either deterministic or statistical methods. The simplest deterministic method for computing response sensitivities is to use finite-difference schemes in conjunction with re-computations using the model with “judiciously chosen” altered parameter values. Evidently, such methods can at best compute approximate values of a very limited number of sensitivities. Deterministic methods that can compute more exactly the values of first-order sensitivities include the “Green’s function method” (Kramer et al., 1981), the “forward sensitivity analysis methodology” (Cacuci, 1981), and the “direct method” (Dunker, 1984), which rely on analytical or numerical differentiation of the computational model under investigation to compute local response sensitivities exactly. On the other hand, “statistical methods” construct an approximate response distribution (often called “response surface”) in the parameters space, and subsequently use scatter plots, regression, rank transformation, correlations, and so-called “partial correlation analysis,” in order to identify approximate expectation values, variances and covariances for the responses. These statistical quantities are subsequently used to construct quantities that play the role of approximate first-order response sensitivities. Thus, statistical methods commence with “uncertainty analysis” and subsequently attempt an approximate “sensitivity analysis” of the approximately computed model response (called a “response surface”) in the phase-space of the parameters under consideration. The currently popular statistical methods for uncertainty and sensitivity analysis are broadly categorized as sampling-based methods (Iman et al., 1981a; Iman et al., 1981b), variance-based methods (Cukier et al., 1978; Hora and Iman, 1986), and Bayesian methods (Rios Insua, 1990). Various variants of the statistical methods for uncertainty and sensitivity analysis are reviewed in the book edited by Saltarelli et al. (2000).

For a computational model comprising many parameters, the conventional deterministic and statistical methods become impractical for computing sensitivities higher than first-order because they are subject to the “curse of dimensionality,” a term coined by Bellman (1957) to describe phenomena in which the number of computations increases exponentially in the respective phase-space. It is known that the “adjoint method of sensitivity analysis” has been the most efficient method for computing exactly first-order sensitivities, since it requires a single large-scale (adjoint) computation for computing all of the first-order sensitivities, regardless of the number of model parameters. The idea underlying the computation of response sensitivities with respect to model parameters using adjoint operators was first used by Wigner (1945) to analyze first-order perturbations in nuclear reactor physics and shielding models based on the linear neutron transport (or diffusion) equation, as subsequently described in textbooks on these subjects (Weiberg and Wigner, 1958; Weisbin et al., 1978; Williams, 1986; Shultis and Faw, 2000; Stacey, 2001). Cacuci (1981) is credited (see, e.g., Práger and Kelemen, 2014; Luo et al., 2020) for having conceived the rigorous “1st-order adjoint sensitivity analysis methodology” for generic large-scale nonlinear (as opposed to linearized) systems involving generic operator responses and having introduced these principles to the earth, atmospheric and other sciences.

Cacuci (2015), Cacuci (2016) has extended his 1st-order adjoint sensitivity analysis methodology to enable the comprehensive and exact computation of 2nd-order sensitivities of model responses to model parameters (including imprecisely known domain boundaries and interfaces) for large-scale linear and nonlinear systems. The unparalleled efficiency of the 2nd-order adjoint sensitivity analysis methodology for linear systems (Cacuci, 2015) was demonstrated (see Cacuci and Fang, 2023, and references therein) by applying this methodology to compute exactly the 21,976 first-order sensitivities and 482,944,576 second-order sensitivities (of which 241,483,276 are distinct from each other) for an OECD/NEA reactor physics benchmark (Valentine, 2006). This benchmark is modeled by the neutron transport equation involving 21,976 uncertain parameters, the solving of which is representative of “large-scale computations.” The neutron transport equation was solved using the software package PARTISN (Alcouffe et al., 2008) in conjunction with the MENDF71X cross section library (Conlin et al., 2013), which comprises 618-group cross sections based on ENDF/B-VII.1 nuclear data (Chadwick et al., 2011). The spontaneous fission source was computed using the code SOURCES4C (Wilson et al., 2002). Contrary to the widely held belief that second- and higher-order sensitivities are negligeable for reactor physics systems, it was found (see Cacuci and Fang, 2023, and references therein) that many 2nd-order sensitivities of this OECD/NEA benchmark’s leakage response to the benchmark’s uncertain parameters were much larger than the largest 1st-order ones, which motivated the investigation of the largest 3rd-order sensitivities, many of which were found to be even larger than the 2nd-order ones. This finding has motivated the development of the mathematical framework for determining and computing the 4th-order sensitivities, many of which were found to be larger than the 3rd-order ones. This sequence of findings has motivated the development by Cacuci (2022) of the “nth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems” (which is abbreviated as “nth-CASAM-L”). The “nth-CASAM-L” mathematical framework was developed specifically for linear systems because the most important model responses produced by such systems can depend simultaneously on both the forward and adjoint state functions governing the respective linear system. Among the most important responses of linear systems that involve both the forward and adjoint functions are various Lagrangian functionals, the Raleigh quotient for computing eigenvalues and/or separation constants when solving partial differential equations, and the Schwinger functional for first-order “normalization-free” solutions (see, e.g., Lewins, 1965; Williams and Engle, 1977; Stacey, 2001). These functionals play fundamental roles in optimization and control procedures, derivation of numerical methods for solving equations (differential, integral, integro-differential), etc. Nonlinear operators do not admit adjoint operators, so responses in nonlinear systems can only depend on the system’s forward state functions. Therefore, the sensitivity analysis of responses that simultaneously involve both forward and adjoint state functions makes it necessary to treat linear models/systems in their own right, rather than treating them as particular cases of nonlinear systems.

The traditional methods of sensitivity analysis aim at computing sensitivities of responses directly to the primary parameters (i.e., microscopic cross sections, isotopic number densities, etc.) involved in the computational model of the physical system under consideration. Although the sensitivities to the primary model parameters are ultimately of interest for subsequent use in predictive modeling activities (which includes the quantification of the uncertainties induced in responses by uncertainties in the primary model parameters, assimilation of experimental data for calibrating the model’s parameters and improving the model’s predictions), the primary parameters seldom appear explicitly in the equations underlying the model. For example, the primary model parameters (e.g., microscopic cross sections, atomic number densities) do not appear explicitly in the forward and adjoint transport equations modeling (Cacuci and Fang, 2023) the above-mentioned OECD/NEA reactor physics benchmark. What appear explicitly in these equations are the macroscopic cross sections, which are functions of the primary model parameters, and which can be considered to be features of the transport equation. This fact has motivated the development by Cacuci (2024a), Cacuci (2024b) of the “nth-Order Features Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (nth-FASAM-N),” which significantly reduces the computational effort computing efficiently and exactly sensitivities of model responses to model features (i.e., functions of the primary model responses), and subsequently compute the sensitivities to responses to the primary model parameters by using the sensitivities to the model features.

Paralleling the mathematical framework of the nth-FASAM-N, it is the purpose of this work to develop a methodology which will enable the efficient and exact computation of sensitivities of model responses to model features for response-coupled forward and adjoint linear systems; this new methodology will be abbreviated as the “nth-FASAM-L” methodology. The mathematical framework of this methodology is established in Section 2 of this work by using the proof by “mathematical induction” as follows: (i) establish the mathematical framework underlying the nth-CASAM-L for n = 1; (ii) assume that the mathematical framework is valid for an arbitrarily high-order, n; (iii) prove that the mathematical framework proposed for n is also valid for n+1. Section 3 presents a concluding discussion that prepares the ground for an illustrative application of the nth-FASAM-L methodology to a representative energy-dependent neutron slowing down model of fundamental importance to reactor physics, which will be presented in an accompanying manuscript (because of word limitations per article), designated as “Part II (Cacuci, 2024c).”

2 The Nth-order function/feature adjoint sensitivity analysis methodology for response-coupled forward and adjoint linear systems (Nth-FASAM-L)

The mathematical framework of the nth-FASAM-L methodology, to be presented in this Section, was established while striving to maximize the computational efficiency of the mathematical framework of the “nth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Coupled Forward/Adjoint Linear Systems” (abbreviated as: nth-CASAM-L)” conceived by Cacuci (2022). The starting point for both the nth-CASAM-L and the nth-FASAM-L is the generic mathematical modeling of a response-coupled forward/adjoint linear system, which is presented in Section 2.1, for convenient referencing.

The validity of mathematical framework underlying the nth-FASAM-L methodology will be established in this Section by using the “proof by mathematical induction” comprising the usual steps, as follows:

1. Conjecture the pattern underlying the nth-FASAM-L, for arbitrary n, based on prior experience.

2. Prove that the conjectured pattern for arbitrary n, is valid for the lowest value of $n$, i.e., for $n=1$.

3. Assuming that that the pattern underlying the nth-FASAM-L is valid for an arbitrarily high-order $n$, prove that this pattern is also valid for $n\to n+1$, i.e., for the (n + 1)^th-FASAM-L.

2.1 Mathematical modeling of response-coupled linear forward and adjoint systems establishing the mathematical framework of the nth-FASAM-L methodology

The mathematical model of a process and/or state of a physical system comprises equations that relate the system’s independent variables and parameters to the system’s state/dependent variables. A linear physical system can generally be modeled by a system of coupled equations written generically in operator form as follows:

$\mathbf{L}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]\mathbf{\phi }\left(\mathbf{x}\right)=\mathbf{Q}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right],\mathbf{x}\in \mathrm{\Omega }\left(\mathbf{\alpha }\right).(1)$

The quantities that appear in Eq. 1 are defined as follows:

1. The vector $\mathbf{\phi }\left(\mathbf{x}\right)\triangleq {\left[{\phi }_1\left(\mathbf{x}\right),\dots ,{\phi }_{TD}\left(\mathbf{x}\right)\right]}^{†}$ is a $TD$-dimensional column vector of dependent variables and where the sub/superscript “$TD$” denotes the “Total (number of) Dependent variables.” The functions ${\phi }_i\left(\mathbf{x}\right),i=1,...,TD$, denote the system’s “dependent variables” (also called “state functions”). The symbol “$\triangleq$” denotes “is defined as” or “is by definition equal to.” Transposition is indicated by a dagger $\left(†\right)$ superscript.

2. The components of the vector $\mathbf{\alpha }\triangleq {\left({\alpha }_1,\dots ,{\alpha }_{TP}\right)}^{†}\in {\mathbb{R}}^{TP}$ denote the primary model parameters, where the subscript/superscript “TP” indicates “Total number of Primary Parameters” and where ${\mathbb{R}}^{TP}$ denotes the TP-dimensional subset of the set of real scalars. Without loss of generality, the model parameters can be considered to be real scalar quantities, having known nominal (or mean) values and, possibly, known 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. These imprecisely known model parameters 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. The nominal parameter values will be denoted as ${\mathbf{\alpha }}^0\triangleq {\left[{\alpha }_1^0,...,{\alpha }_i^0,..,{\alpha }_{TP}^0\right]}^{†}$; the superscript “0” will be used throughout this work to denote “nominal” or “mean” values.

3. The components of the $TI$-dimensional column vector $\mathbf{x}\triangleq {\left(x_1,\dots ,x_{TI}\right)}^{†}\in {\mathbb{R}}^{TI}$ denote the model’s independent variables $x_i,i=1,...,TI$, where the sub/superscript “$TI$” denotes the “Total number of Independent variables.” The vector $\mathbf{x}\in {\mathbb{R}}^{TI}$ of independent variables is considered to be defined on a phase-space domain $\mathrm{\Omega }\left(\mathbf{\alpha }\right)\triangleq \left\{-\infty \le {\lambda }_i\left(\mathbf{\alpha }\right)\le x_i\le {\omega }_i\left(\mathbf{\alpha }\right)\le \infty ;i=1,...,TI\right\}$, the boundaries of which may depend on some of the model parameters $\mathbf{\alpha }$. The lower boundary-point of an independent variable is denoted as ${\lambda }_i\left(\mathbf{\alpha }\right)$, while the corresponding upper boundary-point is denoted as ${\omega }_i\left(\mathbf{\alpha }\right)$. The boundary of $\mathrm{\Omega }\left(\mathbf{\alpha }\right)$, which will be denoted as $\partial \mathrm{\Omega }\left[\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right]$, comprises the set of all of the endpoints ${\lambda }_i\left(\mathbf{\alpha }\right),{\omega }_i\left(\mathbf{\alpha }\right),i=1,...,TI,$ of the respective intervals on which the components of $\mathbf{x}$ are defined, i.e., $\partial \mathrm{\Omega }\left[\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right]\triangleq \left\{{\lambda }_i\left(\mathbf{\alpha }\right)\cup {\omega }_i\left(\mathbf{\alpha }\right),i=1,...,TI\right\}$.

4. The components ${\mathbf{L}}_{ij}\left(\mathbf{x};\mathbf{\alpha }\right)$ of the $TD\times TD$-dimensional matrix $\mathbf{L}\left(\mathbf{x};\mathbf{\alpha }\right)\triangleq \left[L_{ij}\left(\mathbf{x};\mathbf{\alpha }\right)\right],$ $i,j=1,...,TD$, are operators that act linearly on the dependent variables ${\phi }_j\left(\mathbf{x}\right)$ and also depend on the uncertain model parameters $\mathbf{\alpha }$.

5. The vector $\mathbf{g}\left(\mathbf{\alpha }\right)\triangleq \left[g_1\left(\mathbf{\alpha }\right),\dots ,g_{TG}\left(\mathbf{\alpha }\right)\right]$ is a $TG$-dimensional vector having components $g_i\left(\mathbf{\alpha }\right),i=1,...,TG$, which are real-valued functions of (some of) the primary model parameters $\mathbf{\alpha }\in {\mathbb{R}}^{TP}$. The quantity $TG$ denotes the total number of such functions which appear exclusively in the definition of the model’s underlying equations. Such functions customarily appear in models in the form of correlations that describe “features” of the system under consideration, such as material properties, flow regimes. etc. Usually, the number of functions $g_i\left(\mathbf{\alpha }\right)$ is considerably smaller than the total number of model parameters, i.e., $TG\ll TP$. For example, the numerical model (Cacuci and Fang, 2023) of the OECD/NEA “Polyethylene-Reflected Plutonium” reactor physics benchmark (Valentine, 2006) comprises 21,976 uncertain primary model parameters (including microscopic cross sections and isotopic number densities) but the neutron transport equation, which is solved numerically to determine the neutron flux distribution within the benchmark, does not use these primary parameters directly but instead uses just several hundreds of “group-averaged macroscopic cross sections” which are functions/features of the microscopic cross sections and isotopic number densities (which in turn are uncertain quantities that would be components of the vector of primary model parameters). In particular, a component $g_j\left(\mathbf{\alpha }\right)$ may simply be one of the primary model parameters ${\alpha }_j$, i.e., $g_j\left(\mathbf{\alpha }\right)\equiv {\alpha }_j$.

6. The $TD$-dimensional column vector $\mathbf{Q}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]\triangleq {\left(q_1,....,q_{TD}\right)}^{†}$, having components $q_i\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right],i=1,...,TD$, denotes inhomogeneous source terms, which usually depend nonlinearly on the uncertain parameters $\mathbf{\alpha }$. Since the right-side of Eq. 1 may contain distributions, the equality in this equation is considered to hold in the weak (i.e., “distributional”) sense. Similarly, all of the equalities that involve differential equations in this work will be considered to hold in the distributional sense.

7. When $\mathbf{L}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]$ contains differential operators, corresponding boundary and initial conditions which define the domain of $\mathbf{L}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]$ must also be given. Since the complete mathematical model is considered to be linear in $\mathbf{\phi }\left(\mathbf{x}\right)$, the boundary and/or initial conditions needed to define the domain of $\mathbf{L}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]$ must also be linear in $\mathbf{\phi }\left(\mathbf{x}\right)$. Such linear boundary and initial conditions are represented in the following operator form:

$\mathbf{B}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right]\mathbf{\phi }\left(\mathbf{x}\right)=\mathbf{C}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right],\mathbf{x}\in \partial \mathrm{\Omega }\left[\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right](2)$

In Eq. 2, the quantity $\mathbf{B}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right]$ denotes a matrix of dimensions $N_B\times TD$ having components denoted as $B_{ij}\left(\mathbf{x};\mathbf{\alpha }\right);i=1,...,N_B;j=1,...,TD$, which are operators that act linearly on $\mathbf{\phi }\left(\mathbf{x}\right)$ and nonlinearly on the components of $\mathbf{g}\left(\mathbf{\alpha }\right)$; the quantity $N_B$ denotes the total number of boundary and initial conditions. The $N_B$-dimensional column vector $\mathbf{C}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right]$ comprises components that are operators which, in general, act nonlinearly on the components of $\mathbf{g}\left(\mathbf{\alpha }\right)$.

Physical problems modeled by linear systems and/or operators are naturally defined in Hilbert spaces. The dependent variables ${\phi }_i\left(\mathbf{x}\right),i=1,...,TD$, for the physical system represented by Eqs 1, 2 are considered to be square-integrable functions of the independent variables and are considered to belong to a Hilbert space which will be denoted as ${\mathrm{H}}_0\left(\mathrm{\Omega }\right)$, where the subscript “zero” denotes “zeroth-level” or “original.” Higher-level Hilbert spaces, which will be denoted as ${\mathrm{H}}_1\left(\mathrm{\Omega }\right),{\mathrm{H}}_2\left(\mathrm{\Omega }\right),$ etc., will also be introduced and used in this work. The Hilbert space ${\mathrm{H}}_0\left(\mathrm{\Omega }\right)$ is considered to be endowed with the following inner product, denoted as ${〈\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right)〉}_0$, between two elements $\mathbf{\phi }\left(\mathbf{x}\right)\in {\mathrm{H}}_0\left(\mathrm{\Omega }\right)$ and $\mathbf{\psi }\left(\mathbf{x}\right)\in {\mathrm{H}}_0\left(\mathrm{\Omega }\right)$:

${〈\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right)〉}_0\triangleq {\displaystyle \prod _{i=1}^{TI}}{\displaystyle \underset{{\lambda }_i\left(\mathbf{\alpha }\right)}{\overset{{\omega }_i\left(\mathbf{\alpha }\right)}{\int }}}\mathbf{\phi }\left(\mathbf{x}\right)·\mathbf{\psi }\left(\mathbf{x}\right)d\mathbf{x}={\displaystyle \sum _{j=1}^{TD}}{\displaystyle \underset{{\lambda }_1\left(\mathbf{\alpha }\right)}{\overset{{\omega }_1\left(\mathbf{\alpha }\right)}{\int }}}...{\displaystyle \underset{{\lambda }_i\left(\mathbf{\alpha }\right)}{\overset{{\omega }_i\left(\mathbf{\alpha }\right)}{\int }}}...{\displaystyle \underset{{\lambda }_{TI}\left(\mathbf{\alpha }\right)}{\overset{{\omega }_{TI}\left(\mathbf{\alpha }\right)}{\int }}}{\phi }_j\left(\mathbf{x}\right){\psi }_j\left(\mathbf{x}\right)d{x}_1...d{x}_i...d{x}_{TI}.(3)$

The “dot” in Eq. 3 indicates the “scalar product of two vectors,” which is defined in Eq. 4, below, as follows:

$\mathbf{\phi }\left(\mathbf{x}\right)·\mathbf{\psi }\left(\mathbf{x}\right)\triangleq {\displaystyle \sum _{i=1}^{TD}}{\phi }_i\left(\mathbf{x}\right){\psi }_i\left(\mathbf{x}\right).(4)$

The product-notation ${\displaystyle \prod _{i=1}^{TI}}{\displaystyle \underset{{\lambda }_i\left(\mathbf{\alpha }\right)}{\overset{{\omega }_i\left(\mathbf{\alpha }\right)}{\int }}}\left[\right]d{x}_i$ in Eq. 3 denotes the respective multiple integrals.

The linear operator $\mathbf{L}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]$ admits an adjoint operator, which will be denoted as ${\mathbf{L}}^{*}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]$ and which is defined through the following relation for a vector $\mathbf{\psi }\left(\mathbf{x}\right)\in {\mathrm{H}}_0$:

$\begin{array}{l}{〈\mathbf{\psi }\left(\mathbf{x}\right),\mathbf{L}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]\mathbf{\phi }\left(\mathbf{x}\right)〉}_0={〈{\mathbf{L}}^{*}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]\mathbf{\psi }\left(\mathbf{x}\right),\mathbf{\phi }\left(\mathbf{x}\right)〉}_0\\ ={〈\mathbf{\psi }\left(\mathbf{x}\right),\mathbf{Q}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]〉}_0={〈{\mathbf{Q}}^{*}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right],\mathbf{\phi }\left(\mathbf{x}\right)〉}_0.\end{array}(5)$

In Eq. 5, the formal adjoint operator ${\mathbf{L}}^{*}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]$ is the $TD\times TD$ matrix comprising elements $L_{ji}^{*}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]$ which are obtained by transposing the formal adjoints of the forward operators $L_{ij}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]$. Hence, the system adjoint to the linear system represented by Eqs 1, 2 can generally be represented as follows:

${\mathbf{L}}^{*}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]\mathbf{\psi }\left(\mathbf{x}\right)={\mathbf{Q}}^{*}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right],\mathbf{x}\in \mathrm{\Omega }\left(\mathbf{\alpha }\right),(6)$

${\mathbf{B}}^{*}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right]\mathbf{\psi }\left(\mathbf{x}\right)={\mathbf{C}}^{*}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right],\mathbf{x}\in \partial \mathrm{\Omega }\left[\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right].(7)$

When the forward operator $\mathbf{L}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]$ comprises differential operators, the operations (e.g., integration by parts) that implement the transition from the left-side to the right side of Eq. 5 give rise to boundary terms which are collectively called the “bilinear concomitant.” The domain of ${\mathbf{L}}^{*}\left[\mathbf{x};\mathbf{g}\left(\mathbf{\alpha }\right)\right]$ is determined by selecting adjoint boundary and/or initial conditions so as to ensure that the bilinear concomitant vanishes when the selected adjoint boundary conditions are implemented together with the forward boundary conditions given in Eq. 2. The adjoint boundary conditions thus selected are represented in operator form by Eq. 7.

The results computed using a mathematical model are customarily called “model responses” (or “system responses” or “objective functions” or “indices of performance”). For linear physical systems, the system’s response may depend not only on the model’s state-functions and on the system parameters but may simultaneously also depend on the adjoint state function. As has been discussed by Cacuci (2022, 2023a), Cacuci D. G. (2023), any response of a linear system can be formally represented (using expansions or interpolation, if necessary) and fundamentally analyzed in terms of the following generic integral representation:

$R\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]\triangleq {\displaystyle \underset{{\lambda }_1\left(\mathbf{\alpha }\right)}{\overset{{\omega }_1\left(\mathbf{\alpha }\right)}{\int }}}...{\displaystyle \underset{{\lambda }_{TI}\left(\mathbf{\alpha }\right)}{\overset{{\omega }_{TI}\left(\mathbf{\alpha }\right)}{\int }}}S\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{h}\left(\mathbf{\alpha }\right);\mathbf{x}\right]d{x}_1...d{x}_{TI},(8)$

where $S\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{h}\left(\mathbf{\alpha }\right);\mathbf{x}\right]$ is a suitably differentiable nonlinear function of $\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right)$, and $\mathbf{\alpha }$. The integral representation of the response provided in Eq. 8 can represent “averaged” and/or “point-valued” quantities in the phase-space of independent variables. For example, if $R\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ represents the computation or the measurement (which would be a “detector-response”) of a quantity of interest at a point ${\mathbf{x}}_d$ in the phase-space of independent variables, then $S\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{h}\left(\mathbf{\alpha }\right);\mathbf{x}\right]$ would contain a Dirac-delta functional of the form $\delta \left(\mathbf{x}-{\mathbf{x}}_d\right)$. Responses that represent “differentials/derivatives of quantities” would contain derivatives of Dirac-delta functionals in the definition of $S\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{h}\left(\mathbf{\alpha }\right);\mathbf{x}\right]$. The vector $\mathbf{h}\left(\mathbf{\alpha }\right)\triangleq \left[h_1\left(\mathbf{\alpha }\right),\dots ,h_{TH}\left(\mathbf{\alpha }\right)\right]$, having components $h_i\left(\mathbf{\alpha }\right),i=1,...,TH$, which appears among the arguments of the function $S\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{h}\left(\mathbf{\alpha }\right);\mathbf{x}\right]$, represents functions of primary parameters that often appear solely in the definition of the response but do not appear in the mathematical definition of the model, i.e., in Eqs 1, 2, 6, 7. The quantity $TH$ denotes the total number of such functions which appear exclusively in the definition of the model’s response. Evidently, the response will depend directly and/or indirectly (through the “feature”-functions) on all of the primary model parameters. This fact has been indicated in Eq. 8 by using the vector-valued function $\mathbf{f}\left(\mathbf{\alpha }\right)$ as an argument in the definition of the response $R\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ to represent the concatenation of all of the “features” of the model and response under consideration. The vector $\mathbf{f}\left(\mathbf{\alpha }\right)$ of “model features” is thus defined as follows:

$\mathbf{f}\left(\mathbf{\alpha }\right)\triangleq {\left[\mathbf{g}\left(\mathbf{\alpha }\right);\mathbf{h}\left(\mathbf{\alpha }\right);\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right]}^{†}\triangleq {\left[f_1\left(\mathbf{\alpha }\right),...,f_{TF}\left(\mathbf{\alpha }\right)\right]}^{†};TF\triangleq TG+TH+2TI.(9)$

As defined in Eq. 9, the quantity $TF$ denotes the total number of “feature functions of the model’s parameters” which appear in the definition of the nonlinear model’s underlying equations and response.

Solving Eqs 1, 2, at the nominal (or mean) values, denoted as ${\mathbf{\alpha }}^0\triangleq {\left[{\alpha }_1^0,...,{\alpha }_i^0,..,{\alpha }_{TP}^0\right]}^{†}$, of the model parameters, yields the nominal forward solution, which will be denoted as ${\mathbf{\phi }}^0\left(\mathbf{x}\right)$. Solving Eqs 6, 7 at the nominal values, ${\mathbf{\alpha }}^0$, of the model parameters yields the nominal adjoint solution, which will be denoted as ${\mathbf{\psi }}^0\left(\mathbf{x}\right)$. The nominal value of the response, $R\left[{\mathbf{\phi }}^0\left(\mathbf{x}\right),{\mathbf{\psi }}^0\left(\mathbf{x}\right);\mathbf{f}\left({\mathbf{\alpha }}^0\right)\right]$, is determined by using the nominal parameter values ${\mathbf{\alpha }}^0$, the nominal value ${\mathbf{\phi }}^0\left(\mathbf{x}\right)$ of the forward state function, and the nominal value ${\mathbf{\psi }}^0\left(\mathbf{x}\right)$ of the adjoint state function.

The definition provided by Eq. 8 implies that the model response $R\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ depends on the components of the feature function $\mathbf{f}\left(\mathbf{\alpha }\right)$, and would therefore admit a Taylor-series expansion around the nominal value ${\mathbf{f}}^0\triangleq \mathbf{f}\left({\mathbf{\alpha }}^0\right)$, having the following form:

$R\left[\mathbf{f}\left(\mathbf{\alpha }\right)\right]=R\left({\mathbf{f}}^0\right)+{\displaystyle \sum _{j_1=1}^{TF}}{\left\{\frac{\partial R\left(\mathbf{f}\right)}{\partial f_{j_1}}\right\}}_{{\mathbf{f}}^0}\delta f_{j_1}+\frac{1}{2}{\displaystyle \sum _{j_1=1}^{TF}}{\displaystyle \sum _{j_2=1}^{TF}}{\left\{\frac{{\partial }^{2}R\left(\mathbf{f}\right)}{\partial f_{j_1}\partial f_{j_2}}\right\}}_{{\mathbf{f}}^0}\delta f_{j_1}\delta f_{j_2}+...(10)$

where $\delta f_j\triangleq \left[f_j\left(\mathbf{\alpha }\right)-f_j^0\right];f_j^0\triangleq f_j\left({\mathbf{\alpha }}^0\right);j=1,...,TF$. The “sensitivities of the model response with respect to the (feature) functions” are naturally defined as being the functional derivatives of $R\left[\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ with respect to the components (“features”) $f_j\left(\mathbf{\alpha }\right)$ of $\mathbf{f}\left(\mathbf{\alpha }\right)$. The notation ${\left\{·\right\}}_{{\mathbf{f}}^0}$ indicates that the quantity enclosed within the braces is to be evaluated at the nominal values ${\mathbf{f}}^0\triangleq \mathbf{f}\left({\mathbf{\alpha }}^0\right)$. Since $TF\ll TP$, there will be fewer derivatives of the response with respect to the feature functions than there are response derivatives with respect to the primary model parameters. Hence, the computations of the functional derivatives of $R\left[\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ with respect to the functions $f_j\left(\mathbf{\alpha }\right)$, which appear in Eq. 10, will be considerably less expensive computationally than the computation of the functional derivatives involved in the Taylor-series of the response with respect to the model parameters. The functional derivatives of the response with respect to the primary parameters can be obtained from the functional derivatives of the response with respect to the “feature” functions $f_j\left(\mathbf{\alpha }\right)$ by simply using the chain rule, i.e.,:

${\left\{\frac{\partial R\left(\mathbf{\alpha }\right)}{\partial {\alpha }_{j_1}}\right\}}_{{\mathbf{\alpha }}^0}={\displaystyle \sum _{i_1=1}^{TF}}{\left\{\frac{\partial R\left(\mathbf{f}\right)}{\partial f_{i_1}}\frac{\partial f_{i_1}\left(\mathbf{\alpha }\right)}{\partial {\alpha }_{j_1}}\right\}}_{{\mathbf{\alpha }}^0};{\left\{\frac{{\partial }^{2}R\left(\mathbf{\alpha }\right)}{\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}}\right\}}_{{\mathbf{\alpha }}^0}=\frac{\partial }{\partial {\alpha }_{j_2}}{\displaystyle \sum _{i_1=1}^{TF}}{\left\{\frac{\partial R\left(\mathbf{f}\right)}{\partial f_{i_1}}\frac{\partial f_{i_1}\left(\mathbf{\alpha }\right)}{\partial {\alpha }_{j_1}}\right\}}_{{\mathbf{\alpha }}^0};(11)$

and so on. The evaluation/computation of the functional derivatives $\partial f_{i_1}\left(\mathbf{\alpha }\right)/\partial {\alpha }_{j_1}$, ${\partial }^2{f}_{i_1}\left(\mathbf{\alpha }\right)/\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}$, etc., does not require computations involving the model, and is therefore trivial (computationally) by comparison to the evaluation of the functional derivatives (“sensitivities”) of the response with respect to either the functions (“features”) $f_j\left(\mathbf{\alpha }\right)$ or the model parameters ${\alpha }_i,i=1,...,TP$.

The range of validity of the Taylor-series shown in Eq. 10 is defined by its radius of convergence. The accuracy −as opposed to the “validity”− of the Taylor-series in predicting the value of the response at an arbitrary point in the phase-space of model parameters depends on the order of sensitivities retained in the Taylor-expansion: the higher the respective order, the more accurate the respective response value predicted by the Taylor-series. In the particular cases when the response happens to be a polynomial function of the “feature” functions $f_j\left(\mathbf{\alpha }\right)$, the Taylor series is actually exact.

In turn, the functions $f_i\left(\mathbf{\alpha }\right)$ can also be formally expanded in a multivariate Taylor-series around the nominal (mean) parameter values ${\mathbf{\alpha }}^0$, namely,:

$\begin{array}{l}{f}_i\left(\mathbf{\alpha }\right)=f_i\left({\mathbf{\alpha }}^0\right)+{\displaystyle \sum _{j_1=1}^{TP}}{\left\{\frac{\partial f_i\left(\mathbf{\alpha }\right)}{\partial {\alpha }_{j_1}}\right\}}_{{\mathbf{\alpha }}^0}\delta {\alpha }_{j_1}+\frac{1}{2}{\displaystyle \sum _{j_1=1}^{TP}}{\displaystyle \sum _{j_2=1}^{TP}}{\left\{\frac{{\partial }^2{f}_i\left(\mathbf{\alpha }\right)}{\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}}\right\}}_{{\mathbf{\alpha }}^0}\delta {\alpha }_{j_1}\delta {\alpha }_{j_2}\\ +\frac{1}{3!}{\displaystyle \sum _{j_1=1}^{TP}}{\displaystyle \sum _{j_2=1}^{TP}}{\displaystyle \sum _{j_3=1}^{TP}}{\left\{\frac{{\partial }^3{f}_i\left(\mathbf{\alpha }\right)}{\partial {\alpha }_{j_1}\partial {\alpha }_{j_2}\partial {\alpha }_{j_3}}\right\}}_{{\mathbf{\alpha }}^0}\delta {\alpha }_{j_1}\delta {\alpha }_{j_2}\delta {\alpha }_{j_3}+...,\end{array}(12)$

The domain of validity of the Taylor-series in Eq. 12 is defined by its own radius of convergence.

The choice of feature functions $f_i\left(\mathbf{\alpha }\right)$ is not unique but can be tailored by the user to the problem at hand. The two most important guiding principles for constructing the feature functions $f_i\left(\mathbf{\alpha }\right)$ based on the primary parameters are as follows:

(i) As shown in Section 2.2 while establishing the mathematical framework underlying the nth-FASAM-L, the number of large-scale computations needed to determine the numerical value of the second- and higher-order sensitivities is proportional to the number of first-order sensitivities of the model’s response with respect to the feature functions $f_i\left(\mathbf{\alpha }\right)$. Consequently, it is important to minimize the number of feature functions $f_i\left(\mathbf{\alpha }\right)$, while ensuring that all of the primary model parameters are encompassed within the expressions constructed for the feature functions $f_i\left(\mathbf{\alpha }\right)$. In the extreme case when some primary parameters, ${\alpha }_j$, cannot be grouped into the expressions of the feature functions $f_i\left(\mathbf{\alpha }\right)$, then each of the respective primary model parameters ${\alpha }_j$ becomes a feature function $f_j\left(\mathbf{\alpha }\right)$.

(ii) The expressions of the features functions $f_i\left(\mathbf{\alpha }\right)$ must be independent of the model’s state functions; they must be exact, closed-form, scalar-valued functions of the primary model parameters ${\alpha }_j$, so the exact expressions of the derivatives of $f_i\left(\mathbf{\alpha }\right)$ with respect to the primary model parameters ${\alpha }_j$ can be obtained analytically (with “pencil and paper”). The motivation for this requirement is to ensure that the numerical determination of the subsequent derivatives of the features functions $f_i\left(\mathbf{\alpha }\right)$ with respect to the primary model parameters ${\alpha }_j$ becomes trivial computationally. In the extreme case when no feature function can be constructed, the feature functions are the primary parameters themselves, in which case the nth-FASAM-L methodology becomes identical to the previously established nth-CASAM-L methodology (Cacuci, 2022)

2.2 Establishing the mathematical framework of the nth-FASAM-L methodology

Cacuci D. G. (2023), Cacuci (2024a), Cacuci (2024b) has recently developed the “nth-Order Features Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (nth-FASAM-N)” which enables the computation of arbitrarily-high-order sensitivities of responses to features/functions of parameters for nonlinear models/systems. Together, the nth-CASAM-L and the nth-FASAM-N provide the basis for the development of the “nth-Order Features Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward and Adjoint Linear Systems (Nth-FASAM-L)” to be presented in this Section. In particular, comparing the mathematical framework of the 1st-FASAM-L to the framework of the 1st-CASAM-L (Cacuci and Fang, 2023) suggests that the components $f_i\left(\mathbf{\alpha }\right),i=1,...,TF$, of the “feature function” $\mathbf{f}\left(\mathbf{\alpha }\right)\triangleq {\left[f_1\left(\mathbf{\alpha }\right),...,f_{TF}\left(\mathbf{\alpha }\right)\right]}^{†}$ will play within the 1st-FASAM-L the same role as played by the components ${\alpha }_j,j=1,...,TP$, of the “vector of primary model parameters” $\mathbf{\alpha }\triangleq {\left({\alpha }_1,\dots ,{\alpha }_{TP}\right)}^{†}$ within the framework of the 1st-CASAM-L. It can therefore be conjectured that the same correspondence would be expected to hold in general, between the general frameworks of the nth-FASAM-L and the nth-CASAM-L methodologies. As will be demonstrated in this Section, this conjecture is indeed correct.

Considering the analogy to the framework of the nth-CASAM-L methodology (Cacuci, 2022), it is conjectured that that the G-differential of the (n-1)^th-order sensitivity of the model’s response $R\left[\mathbf{u}\left(\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ with respect to the components $f_1,...,f_{TF}$ of the “feature” function $\mathbf{f}\left(\mathbf{\alpha }\right)\triangleq {\left[f_1\left(\mathbf{\alpha }\right),...,f_{TF}\left(\mathbf{\alpha }\right)\right]}^{†}$ will have the following form:

$\begin{array}{l}{\left\{\delta R^{\left(n-1\right)}\left[j_{n-1};...;j_1;{\mathbf{u}}^{\left(n\right)};\mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}={\displaystyle \sum _{j_n=1}^{TF}}{\left\{\frac{\partial }{\partial f_{j_n}}{\displaystyle \underset{{\lambda }_1\left(\mathbf{\alpha }\right)}{\overset{{\omega }_1\left(\mathbf{\alpha }\right)}{\int }}}d{x}_1...{\displaystyle \underset{{\lambda }_{TI}\left(\mathbf{\alpha }\right)}{\overset{{\omega }_{TI}\left(\mathbf{\alpha }\right)}{\int }}}d{x}_{TI}{S}^{\left(n-1\right)}\left(j_{n-1},...,j_1;{\mathbf{u}}^{\left(n\right)};\mathbf{\alpha }\right)\delta f_{j_n}\right\}}_{{\mathbf{\alpha }}^0}\\ +{\left\{{〈{\mathbf{a}}^{\left(n\right)}\left(j_{n-1},...,j_1;\mathbf{x}\right),{\mathbf{q}}_V^{\left(n\right)}\left[2^n;{\mathbf{u}}^{\left(n\right)}\left(2^n;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]〉}_n\right\}}_{{\mathbf{\alpha }}^0}\\ -{\left\{{\widehat{P}}^{\left(n\right)}\left({\mathbf{a}}^{\left(n\right)};{\mathbf{u}}^{\left(n\right)};\mathbf{f};\delta \mathbf{f}\right)\right\}}_{{\mathbf{\alpha }}^0}\equiv {\displaystyle \sum _{j_n=1}^{TF}}{\left\{R^{\left(n\right)}\left(j_n,...,j_1;{\mathbf{u}}^{\left(n\right)};{\mathbf{a}}^{\left(n\right)};\mathbf{f}\right)\right\}}_{\left({\mathbf{\alpha }}^0\right)}\delta f_{j_n},\end{array}(13)$

such that the nth-order sensitivity of the model’s response $R\left[\mathbf{u}\left(\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ with respect to the components $f_{j_1},...,f_{j_n}$ of the “feature” function $\mathbf{f}\left(\mathbf{\alpha }\right)\triangleq {\left[f_{j_1}\left(\mathbf{\alpha }\right),...,f_{j_n}\left(\mathbf{\alpha }\right)\right]}^{†}$ is expected to have the following functional form:

$\begin{array}{l}{R}^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n\right)}\left(2^n;j_{n-2},...,j_1;\mathbf{x}\right);{\mathbf{a}}^{\left(n\right)}\left(2^n;j_{n-1},...,j_1;\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]\\ \triangleq {\displaystyle \underset{{\lambda }_1\left(\mathbf{\alpha }\right)}{\overset{{\omega }_1\left(\mathbf{\alpha }\right)}{\int }}}d{x}_1...{\displaystyle \underset{{\lambda }_{TI}\left(\mathbf{\alpha }\right)}{\overset{{\omega }_{TI}\left(\mathbf{\alpha }\right)}{\int }}}d{x}_{TI}{S}^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n\right)};{\mathbf{a}}^{\left(n\right)};\mathbf{f}\left(\mathbf{\alpha }\right)\right]\\ \triangleq {\partial }^{n}R\left[\mathbf{u}\left(\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]/\partial f_{j_n}...\partial f_{j_1};j_1=1,...,TF;j_n=1,...,j_{n-1};n=2,3,...\end{array}(14)$

where $TF$ denotes the “total number of features,” i.e., functions of the primary model parameters.

It is also conjectured that the nth-level adjoint functions ${\mathbf{u}}^{\left(n\right)}\left(2^n;j_{n-2},...,j_1;\mathbf{x}\right)$ and ${\mathbf{a}}^{\left(n\right)}\left(j_{n-1},...,j_1;2^n;\mathbf{x}\right)$, which are needed to compute the nth-order sensitivities shown in Eq. 14, are obtained as follows:

(i) ${\mathbf{u}}^{\left(n\right)}\left(2^n;j_{n-2},...,j_1;\mathbf{x}\right)\triangleq {\left[{\mathbf{u}}^{\left(n-1\right)}\left(2^{n-1};j_{n-3},...,j_1;\mathbf{x}\right),{\mathbf{a}}^{\left(n-1\right)}\left(2^{n-1};j_{n-2},...,j_1;\mathbf{x}\right)\right]}^{†}$ are the solutions of the following nth-Level Forward/Adjoint System (nth-LFAS) for $j_n=1,...,j_{n-1};n\ge 3$:

${\mathbf{F}}^{\left(n\right)}\left[2^n\times 2^n;\mathbf{f}\left(\mathbf{\alpha }\right)\right]{\mathbf{u}}^{\left(n\right)}\left(2^n;j_{n-2},...,j_1;\mathbf{x}\right)={\mathbf{q}}_F^{\left(n\right)}\left[2^n;{\mathbf{u}}^{\left(n-1\right)};\mathbf{f}\left(\mathbf{\alpha }\right)\right];\mathbf{x}\in \mathrm{\Omega };(15)$

${\mathbf{b}}_F^{\left(n\right)}\left(2^n;{\mathbf{u}}^{\left(n\right)};\mathbf{f}\right)=\mathbf{0};\mathbf{x}\in \partial \mathrm{\Omega };(16)$

(ii) ${\mathbf{a}}^{\left(n\right)}\left(j_{n-1},...,j_1;2^n;\mathbf{x}\right)$ are the solutions of the following nth-Level Adjoint Sensitivity System (nth-LASS) for $j_n=1,...,j_{n-1};n\ge 3$:

${\mathbf{A}}^{\left(n\right)}\left[2^n\times 2^n;\mathbf{f}\left(\mathbf{\alpha }\right)\right]{\mathbf{a}}^{\left(n\right)}\left(2^n;j_{n-1},...,j_1;\mathbf{x}\right)={\mathbf{s}}^{\left(n\right)}\left(2^n;{\mathbf{u}}^{\left(n\right)};\mathbf{f}\right);\mathbf{x}\in \mathrm{\Omega };(17)$

${\left\{{\mathbf{b}}_A^{\left(n\right)}\left[{\mathbf{u}}^{\left(n\right)}\left(2^n;j_{n-1},...,j_1;\mathbf{x}\right);{\mathbf{a}}^{\left(n\right)}\left(2^n;j_{n-1},...,j_1;\mathbf{x}\right);\mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}=\mathbf{0},\mathbf{x}\in \partial \mathrm{\Omega }.(18)$

Through their implicit dependence on lower-level forward and adjoint functions, the block-matrix valued operators ${\mathbf{F}}^{\left(n\right)}\left[2^n\times 2^n;\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ and ${\mathbf{A}}^{\left(n\right)}\left[2^n\times 2^n;\mathbf{f}\left(\mathbf{\alpha }\right)\right]$, as well as the source terms ${\mathbf{q}}_F^{\left(n\right)}\left[2^n;{\mathbf{u}}^{\left(n-1\right)};\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ and ${\mathbf{s}}^{\left(n\right)}\left(2^n;{\mathbf{u}}^{\left(n\right)};\mathbf{f}\right)$, also depend on lower-level indices $j_k,k<n$, but this dependence is not material to establishing the general framework of the nth-FASAM-L and has therefore been omitted, to keep the notation as simple as possible.

2.3 Proving that the conjectured mathematical framework of the nth-FASAM-L methodology is correct for $n=1$

The proof that the framework conjectured in Section 2.2 for the nth-FASAM-L methodology is indeed correct/valid when $n=1$ (for the 1st-FASAM-L methodology) parallels the proof used in (Cacuci, 2022) to show that the framework of the nth-CASAM-L methodology reduces to the corresponding 1st-CASAM-L methodology when $n=1$. In preparation for subsequent generalizations towards establishing the generic pattern for computing sensitivities of arbitrarily high-order, the function ${\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right)\right]}^{†}$ will be called the “1st-level forward/adjoint function” and the system of equations satisfied by this function (which is obtained by concatenating the original forward and adjoint equations together with their respective boundary/initial conditions) will be called “the 1st-Level Forward/Adjoint System (1st-LFAS)” and will be re-written in the following concatenated matrix-form:

${\mathbf{F}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right)={\mathbf{q}}_F^{\left(1\right)}\left(2;\mathbf{x};\mathbf{f}\right);\mathbf{x}\in \mathrm{\Omega }\left(\mathbf{\alpha }\right);(19)$

${\mathbf{b}}_F^{\left(1\right)}\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\right]=\mathbf{0};\mathbf{x}\in \partial \mathrm{\Omega }\left[\mathbf{\lambda }\left(\mathbf{\alpha }\right);\mathbf{\omega }\left(\mathbf{\alpha }\right)\right];(20)$

where the following definitions were used:

${\mathbf{F}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]\triangleq \left(\begin{array}{cc}\mathbf{L}\left(\mathbf{x};\mathbf{f}\right)& \mathbf{0}\\ \mathbf{0}& {\mathbf{L}}^{*}\left(\mathbf{x};\mathbf{f}\right)\end{array}\right);{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right)\right]}^{†};(21)$

${\mathbf{q}}_F^{\left(1\right)}\left(2;\mathbf{x};\mathbf{f}\right)\triangleq \left(\begin{array}{c}\mathbf{Q}\left(\mathbf{x};\mathbf{g}\right)\\ {\mathbf{Q}}^{*}\left(\mathbf{x};\mathbf{g}\right)\end{array}\right);{\mathbf{b}}_F^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\right]\triangleq \left(\begin{array}{c}\mathbf{B}\left(\mathbf{x};\mathbf{f}\right)\mathbf{\phi }\left(\mathbf{x}\right)-\mathbf{C}\left(\mathbf{f}\right)\\ {\mathbf{B}}^{*}\left(\mathbf{x};\mathbf{f}\right)\mathbf{\psi }\left(\mathbf{x}\right)-{\mathbf{C}}^{*}\left(\mathbf{f}\right)\end{array}\right).(22)$

In the list of arguments of the matrix ${\mathbf{F}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]$, the argument “$2\times 2$” indicates that this square matrix comprises four component sub-matrices, as indicated in Eq. 21. Similarly, the argument “2” that appears in the block-vectors ${\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right)$, ${\mathbf{q}}_F^{\left(1\right)}\left(2;\mathbf{x};\mathbf{f}\right)$, and ${\mathbf{b}}_F^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\right]$ defined in Eq. 22 indicates that each of these column block-vectors comprises two sub-vectors as components. Also, throughout this work, the quantity “$\mathbf{0}$” will be used to denote either as a vector with zero-valued components or a matrix zero-valued components, depending on the context. For example, the vector “$\mathbf{0}$” in Eq. 20 is considered to have as many components as the vector ${\mathbf{b}}_F^{\left(1\right)}\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\right]$. On the other hand, the quantity “$\mathbf{0}$” which appears in Eq. 21 may represent either a (sub) matrix or a vector of the requisite dimensions.

The primary parameters $\mathbf{\alpha }$ are subject to uncertainties; their nominal (or mean) values, denoted as ${\mathbf{\alpha }}^0$, are considered to be known, but these values will differ from the true values $\mathbf{\alpha }$, which are unknown, by variations $\delta \mathbf{\alpha }\triangleq {\left(\delta {\alpha }_1,\dots ,\delta {\alpha }_{TP}\right)}^{†}$, where $\delta {\alpha }_i\triangleq {\alpha }_i-{\alpha }_i^0$. The parameter variations $\delta \mathbf{\alpha }$ will induce variations $\delta \mathbf{f}\left(\mathbf{\alpha }\right)\triangleq {\left[\delta f_1\left(\mathbf{\alpha }\right),\dots ,\delta f_{TF}\left(\mathbf{\alpha }\right)\right]}^{†}$ in the vector-valued “feature” function $\mathbf{f}\left(\mathbf{\alpha }\right)$, around the nominal value ${\mathbf{f}}^0\triangleq \mathbf{f}\left({\mathbf{\alpha }}^0\right)$, and will also induce variations $\delta \mathbf{\phi }\left(\mathbf{x}\right)$ and $\delta \mathbf{\psi }\left(\mathbf{x}\right)$, respectively, around the nominal solutions ${\mathbf{\phi }}^0,{\mathbf{\psi }}^0$, through the equations underlying the model. All of these variations will induce variations in the model response $R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\right]\equiv R\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$.

Formally, the first-order sensitivities of the response $R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\right]$ with respect to the components of the feature function $\mathbf{f}\left(\mathbf{\alpha }\right)$ are provided by the first-order Gateaux (G-)variation of $R\left(\mathbf{\phi },\mathbf{\psi },\mathbf{f}\right)$ at the phase-space point $\left({\mathbf{\phi }}^0,{\mathbf{\psi }}^0,{\mathbf{f}}^0\right)$, which is defined as follows:

$\begin{array}{l}\delta R\left({\mathbf{\phi }}^0,{\mathbf{\psi }}^0,{\mathbf{f}}^0;\delta \mathbf{\phi },\delta \mathbf{\psi },\delta \mathbf{f}\right)\triangleq {\left\{\frac{d}{d\epsilon }R\left[{\mathbf{\phi }}^0\left(\mathbf{x}\right)+\epsilon \delta \mathbf{\phi }\left(\mathbf{x}\right),{\mathbf{\psi }}^0\left(\mathbf{x}\right)+\epsilon \delta \mathbf{\psi }\left(\mathbf{x}\right);{\mathbf{f}}^0+\epsilon \delta \mathbf{f}\right]\right\}}_{\epsilon =0}\\ \equiv {\left\{\frac{d}{d\epsilon }R\left[{\mathbf{u}}^{\left(1,0\right)}\left(2;\mathbf{x}\right)+\epsilon {\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{f}}^0+\epsilon \delta \mathbf{f}\right]\right\}}_{\epsilon =0}\equiv \delta R\left[{\mathbf{u}}^{\left(1,0\right)}\left(2;\mathbf{x}\right);{\mathbf{f}}^0;{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right),\delta \mathbf{f}\right].\end{array}(23)$

The definitions provided in Eq. 24, below, were used in Eq. 23:

${\mathbf{u}}^{\left(1,0\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{\phi }}^0\left(\mathbf{x}\right),{\mathbf{\psi }}^0\left(\mathbf{x}\right)\right]}^{†};{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[\delta \mathbf{\phi }\left(\mathbf{x}\right),\delta \mathbf{\psi }\right]}^{†}.(24)$

The numerical methods (e.g., Newton’s method and variants thereof) for solving large-scale systems require the existence of the first-order G-derivatives of the original model equations and of the model’s response; these will be assumed to exist. When the 1st-order G-derivatives exists, the variation $\delta R\left[{\mathbf{u}}^{\left(1,0\right)}\left(2;\mathbf{x}\right);{\mathbf{f}}^0;{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right),\delta \mathbf{f}\right]$ can be written as follows:

$\begin{array}{l}\delta R\left[{\mathbf{u}}^{\left(1,0\right)}\left(2;\mathbf{x}\right);{\mathbf{f}}^0;{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right),\delta \mathbf{f}\right]={\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{dir}\\ +{\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\right]\right\}}_{ind}.\end{array}(25)$

In Eq. 25, the “direct-effect” term ${\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{dir}$ comprises only dependencies on $\delta \mathbf{f}\left(\mathbf{\alpha }\right)$ and is defined as follows:

${\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{dir}\triangleq {\left\{\frac{\partial R\left({\mathbf{u}}^{\left(1\right)};\mathbf{f}\right)}{\partial \mathbf{f}}\delta \mathbf{f}\right\}}_{{\mathbf{\alpha }}^0}.(26)$

The following convention/definition was used in Eq. 26:

$\frac{\partial \left[\right]}{\partial \mathbf{f}}\delta \mathbf{f}\triangleq {\displaystyle \sum _{i=1}^{TF}}\frac{\partial \left[\right]}{\partial f_i}\delta f_i={\displaystyle \sum _{i=1}^{TG}}\frac{\partial \left[\right]}{\partial g_i}\delta g_i+{\displaystyle \sum _{i=1}^{TH}}\frac{\partial \left[\right]}{\partial h_i}\delta h_i+{\displaystyle \sum _{i=1}^{TI}}\frac{\partial \left[\right]}{\partial {\omega }_i}\delta {\omega }_i+{\displaystyle \sum _{i=1}^{TI}}\frac{\partial \left[\right]}{\partial {\lambda }_i}\delta {\lambda }_i.(27)$

The notation on the left-side of Eq. 27 represents the inner product between two vectors, but the “dagger” symbol “($†$)” which indicates “transposition” has been omitted in order to keep the notation as simple as possible. “Daggers” indicating transposition will also be omitted in other inner products, whenever possible, while avoiding ambiguities.

In Eq. 25, the “indirect-effect” term ${\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\right]\right\}}_{ind}$ depends only on the variations ${\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[\delta \mathbf{\phi }\left(\mathbf{x}\right),\delta \mathbf{\psi }\right]}^{†}$ in the state functions, and is defined as follows:

$\begin{array}{l}{\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\right]\right\}}_{ind}\triangleq {\left\{{\displaystyle \underset{{\lambda }_1\left(\mathbf{\alpha }\right)}{\overset{{\omega }_1\left(\mathbf{\alpha }\right)}{\int }}}d{x}_1...{\displaystyle \underset{{\lambda }_{TI}\left(\mathbf{\alpha }\right)}{\overset{{\omega }_{TI}\left(\mathbf{\alpha }\right)}{\int }}}d{x}_{TI}\frac{\partial S\left(\mathbf{\phi },\mathbf{\psi };\mathbf{g};\mathbf{h}\right)}{\partial {\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right)}{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\right\}}_{{\mathbf{\alpha }}^0}\\ \triangleq {\left\{{\displaystyle \underset{{\lambda }_1\left(\mathbf{\alpha }\right)}{\overset{{\omega }_1\left(\mathbf{\alpha }\right)}{\int }}}d{x}_1...{\displaystyle \underset{{\lambda }_{TI}\left(\mathbf{\alpha }\right)}{\overset{{\omega }_{TI}\left(\mathbf{\alpha }\right)}{\int }}}d{x}_{TI}\frac{\partial S\left(\mathbf{\phi },\mathbf{\psi };\mathbf{g};\mathbf{h}\right)}{\partial \mathbf{\phi }}\delta \mathbf{\phi }\right\}}_{{\mathbf{\alpha }}^0}+{\left\{{\displaystyle \underset{{\lambda }_1\left(\mathbf{\alpha }\right)}{\overset{{\omega }_1\left(\mathbf{\alpha }\right)}{\int }}}d{x}_1...{\displaystyle \underset{{\lambda }_{TI}\left(\mathbf{\alpha }\right)}{\overset{{\omega }_{TI}\left(\mathbf{\alpha }\right)}{\int }}}d{x}_{TI}\frac{\partial S\left(\mathbf{\phi },\mathbf{\psi };\mathbf{g};\mathbf{h}\right)}{\partial \mathbf{\psi }}\delta \mathbf{\psi }\right\}}_{{\mathbf{\alpha }}^0}.\end{array}(28)$

In Eqs 26, 28, the notation ${\left\{\hspace{1em}\right\}}_{{\mathbf{\alpha }}^0}$ has been used to indicate that the quantity within the brackets is to be evaluated at the nominal values of the parameters and state functions. This simplified notation is justified by the fact that when the parameters take on their nominal values, it implicitly means that the corresponding state functions also take on their corresponding nominal values. This simplified notation will be used throughout this work.

The direct-effect term can be computed after having solved the forward system modeled by Eqs 1, 2, as well as the adjoint system modeled by Eqs 6, 7, using the nominal parameter values to obtain the nominal values ${\mathbf{\phi }}^0,{\mathbf{\psi }}^0$ of the forward and adjoint dependent variables.

On the other hand, the indirect-effect term ${\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\right]\right\}}_{ind}$ defined in Eq. 28 can be quantified only after having determined the variations ${\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[\delta \mathbf{\phi }\left(\mathbf{x}\right),\delta \mathbf{\psi }\right]}^{†}$ in the state functions of the 1st-Level Forward/Adjoint System (1st-LFAS). The variations ${\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)$ are obtained as the solutions of the system of equations obtained by taking the first-order G-differentials of the 1st-LFAS defined by Eqs 19, 20, which are obtained by definition as follows:

${\left\{\frac{d}{d\epsilon }{\mathbf{F}}^{\left(1\right)}\left[2\times 2;\mathbf{x};{\mathbf{f}}^0+\epsilon \delta \mathbf{f}\right]\left[{\mathbf{u}}^{\left(1,0\right)}\left(2;\mathbf{x}\right)+\epsilon {\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\right]\right\}}_{\epsilon =0}={\left\{\frac{d}{d\epsilon }{\mathbf{q}}_F^{\left(1\right)}\left[2;\mathbf{x};{\mathbf{f}}^0+\epsilon \delta \mathbf{f}\right]\right\}}_{\epsilon =0},(29)$

${\left\{\frac{d}{d\epsilon }{\mathbf{b}}_F^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1,0\right)}\left(2;\mathbf{x}\right)+\epsilon {\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{f}}^0+\epsilon \delta \mathbf{f}\right]\right\}}_{\epsilon =0}=\mathbf{0}\left[2\right].(30)$

Carrying out the differentiations with respect to $\epsilon$ in Eqs 29, 30 and setting $\epsilon =0$ in the resulting expressions yields the following matrix-vector equations:

${\left\{{\mathbf{V}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\right\}}_{{\mathbf{\alpha }}^0}={\left\{{\mathbf{q}}_V^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0};\mathbf{x}\in \mathrm{\Omega }\left({\mathbf{\alpha }}^0\right);(31)$

${\left\{{\mathbf{b}}_v^{\left(1\right)}\left({\mathbf{u}}^{\left(1\right)};{\mathbf{v}}^{\left(1\right)};\mathbf{f};\delta \mathbf{f}\right)\right\}}_{{\mathbf{\alpha }}^0}=\mathbf{0};\mathbf{x}\in \partial \mathrm{\Omega }\left[\mathbf{\lambda }\left({\mathbf{\alpha }}^0\right);\mathbf{\omega }\left({\mathbf{\alpha }}^0\right)\right];(32)$

where:

${\mathbf{V}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]\triangleq \left(\begin{array}{cc}\mathbf{L}\left(\mathbf{x};\mathbf{f}\right)& \mathbf{0}\\ \mathbf{0}& {\mathbf{L}}^{*}\left(\mathbf{x};\mathbf{f}\right)\end{array}\right)={\mathbf{F}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right];(33)$

${\mathbf{q}}_V^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1\right)};\mathbf{f};\delta \mathbf{f}\right]\triangleq \left(\begin{array}{c}{\mathbf{q}}_1^{\left(1\right)}\left(\mathbf{\phi };\mathbf{f};\delta \mathbf{f}\right)\\ {\mathbf{q}}_2^{\left(1\right)}\left(\mathbf{\psi };\mathbf{f};\delta \mathbf{f}\right)\end{array}\right);{\mathbf{b}}_v^{\left(1\right)}\left({\mathbf{u}}^{\left(1\right)};{\mathbf{v}}^{\left(1\right)};\mathbf{f};\delta \mathbf{f}\right)\triangleq \left(\begin{array}{c}{\mathbf{b}}_1^{\left(1\right)}\left(\mathbf{\phi };\delta \mathbf{\phi };\mathbf{f};\delta \mathbf{f}\right)\\ {\mathbf{b}}_2^{\left(1\right)}\left(\mathbf{\psi };\delta \mathbf{\psi };\mathbf{f};\delta \mathbf{f}\right)\end{array}\right);(34)$

${\mathbf{q}}_1^{\left(1\right)}\left(\mathbf{\phi };\mathbf{f};\delta \mathbf{f}\right)\triangleq \frac{\partial \left[\mathbf{Q}-\mathbf{L}\mathbf{\phi }\left(\mathbf{x}\right)\right]}{\partial \mathbf{f}}\delta \mathbf{f}\triangleq {\displaystyle \sum _{j_1=1}^{TF}}{\mathbf{s}}_1^{\left(1\right)}\left(j_1;\mathbf{\phi };\mathbf{f}\right)\delta f_{j_1}(35)$

${\mathbf{q}}_2^{\left(1\right)}\left(\mathbf{\psi },\mathbf{f};\delta \mathbf{f}\right)\triangleq \frac{\partial \left[{\mathbf{Q}}^{*}-{\mathbf{L}}^{*}\mathbf{\psi }\left(\mathbf{x}\right)\right]}{\partial \mathbf{f}}\delta \mathbf{f}\triangleq {\displaystyle \sum _{j_1=1}^{TF}}{\mathbf{s}}_2^{\left(1\right)}\left(j_1;\mathbf{\psi };\mathbf{f}\right)\delta f_{j_1}(36)$

${\mathbf{b}}_1^{\left(1\right)}\left(\mathbf{\phi };\delta \mathbf{\phi };\mathbf{f};\delta \mathbf{f}\right)\triangleq \mathbf{B}\delta \mathbf{\phi }+\frac{\partial \left(\mathbf{B}\mathbf{\phi }-\mathbf{C}\right)}{\partial \mathbf{f}}\delta \mathbf{f};(37)$

${\mathbf{b}}_2^{\left(1\right)}\left(\mathbf{\psi };\delta \mathbf{\psi };\mathbf{f};\delta \mathbf{f}\right)\triangleq {\mathbf{B}}^{*}\delta \mathbf{\psi }+\frac{\partial \left({\mathbf{B}}^{*}\mathbf{\psi }-{\mathbf{C}}^{*}\right)}{\partial \mathbf{f}}\delta \mathbf{f}.(38)$

In order to keep the notation as simple as possible in Eqs 31‒38, the differentials with respect to the various components of the feature function $\mathbf{f}\left(\mathbf{\alpha }\right)$ have all been written in the form $\left(\partial \left[\right]/\partial \mathbf{f}\right)\delta \mathbf{f}$, keeping in mind the convention/notation introduced in Eq. 27. The system of equations comprising Eqs 31, 32 will be called the “1st-Level Variational Sensitivity System (1st-LVSS)” and its solution, ${\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)$, will be called the “1st-level variational sensitivity function,” which is indicated by the superscript “(1)”. The solution, ${\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)$, of the 1st-LVSS will be a function of the components of the vector of variations $\delta \mathbf{f}$. In principle, therefore, if the response sensitivities with respect to the components of the feature function $\mathbf{f}\left(\mathbf{\alpha }\right)$ are of interest, then the 1st-LVSS would need to be solved as many times as there are components in the variational features-function $\delta \mathbf{f}$. On the other hand, if the response sensitivities with respect to the primary parameters are of interest, then the 1st-LVSS would need to be solved as many times as there are primary parameters. Solving the 1st-LVSS involves “large-scale computations.”

Solving the 1st-LVSS can be avoided altogether by using the ideas underlying the “adjoint sensitivity analysis methodology” originally conceived by Cacuci (1981), and subsequently generalized by Cacuci (2022), Cacuci D. G. (2023) to enable the computation of arbitrarily high-order response sensitivities to primary model parameters for both linear and nonlinear models. Thus, the need for solving repeatedly the 1st-LVSS for every variation in the components of the feature function (or for every variation in the model’s parameters) is eliminated by expressing the indirect-effect term ${\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\right]\right\}}_{ind}$ defined in Eq. 28 in terms of the solutions of the “1st-Level Adjoint Sensitivity System” (1st-LASS), which will be constructed by implementing the following sequence of steps:

1. Introduce a Hilbert space, denoted as ${\mathrm{H}}_1$, comprising vector-valued elements of the form ${\mathbf{\chi }}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{\chi }}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{\chi }}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$, where the components ${\mathbf{\chi }}_i^{\left(1\right)}\left(\mathbf{x}\right)\triangleq {\left[{\chi }_{i,1}^{\left(1\right)}\left(\mathbf{x}\right),...,{\chi }_{i,j}^{\left(1\right)}\left(\mathbf{x}\right),...,{\chi }_{i,TD}^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$, $i=1,2$, are square-integrable functions. Consider further that this Hilbert space is endowed with an inner product denoted as ${〈{\mathbf{\chi }}^{\left(1\right)}\left(2;\mathbf{x}\right),{\mathbf{\theta }}^{\left(1\right)}\left(2;\mathbf{x}\right)〉}_1$ between two elements, ${\mathbf{\chi }}^{\left(1\right)}\left(2;\mathbf{x}\right)\in {\mathrm{H}}_1$, ${\mathbf{\theta }}^{\left(1\right)}\left(2;\mathbf{x}\right)\in {\mathrm{H}}_1$, which is defined as follows:

${〈{\mathbf{\chi }}^{\left(1\right)}\left(2;\mathbf{x}\right),{\mathbf{\theta }}^{\left(1\right)}\left(2;\mathbf{x}\right)〉}_1\triangleq {\displaystyle \sum _{i=1}^2}{〈{\mathbf{\chi }}_i^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{\theta }}_i^{\left(1\right)}\left(\mathbf{x}\right)〉}_0.(39)$

2. In the Hilbert ${\mathrm{H}}_1$, use Eq. 39 to form the inner product of Eq. 31 with a yet undefined vector-valued function ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}\in {\mathrm{H}}_1$ to obtain the following relation:

${\left\{{〈{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right),{\mathbf{V}}^{\left(1\right)}\left[2\times 2;\mathbf{x};{\mathbf{f}}^0\right]{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)〉}_1\right\}}_{{\mathbf{\alpha }}^0}={\left\{{〈{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right),{\mathbf{q}}_V^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]〉}_1\right\}}_{{\mathbf{\alpha }}^0}.(40)$

3. Using the definition of the adjoint operator in the Hilbert space ${\mathrm{H}}_1$, recast the left-side of Eq. 40 as follows:

$\begin{array}{l}{\left\{{〈{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right),{\mathbf{V}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)〉}_1\right\}}_{{\mathbf{\alpha }}^0}={\left\{{〈{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right),{\mathbf{A}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)〉}_1\right\}}_{{\mathbf{\alpha }}^0}\\ +{\left\{P^{\left(1\right)}\left[{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0},\end{array}(41)$

where ${\left\{P^{\left(1\right)}\left[{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}$ denotes the bilinear concomitant defined on the phase-space boundary $\mathbf{x}\in \partial \mathrm{\Omega }\left({\mathbf{\alpha }}^0\right)$, and where ${\mathbf{A}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]$ is the operator formally adjoint to ${\mathbf{V}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]$, as defined in Eq. 42 below:

${\mathbf{A}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]\triangleq {\left\{{\mathbf{V}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]\right\}}^{*}=\left(\begin{array}{cc}{\mathbf{L}}^{*}\left(\mathbf{x};\mathbf{f}\right)& \mathbf{0}\\ \mathbf{0}& \mathbf{L}\left(\mathbf{x};\mathbf{f}\right)\end{array}\right).(42)$

4. Require the first term on right-side of Eq. 41 to represent the indirect-effect term defined in Eq. 28, to obtain the following relation:

${\mathbf{A}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)={\mathbf{q}}_A^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\right],\mathbf{x}\in \mathrm{\Omega }\left({\mathbf{\alpha }}^0\right);(43)$

where the source term on the right-side of Eq. 43 is defined in Eq. 44, below:

${\mathbf{q}}_A^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\right]\triangleq {\left[\frac{\partial S\left({\mathbf{u}}^{\left(1\right)};\mathbf{f}\right)}{\partial {\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right)}\right]}^{†}\triangleq \left(\begin{array}{c}{\left[\partial S\left({\mathbf{u}}^{\left(1\right)};\mathbf{f}\right)/\partial \mathbf{\phi }\right]}^{†}\\ {\left[\partial S\left({\mathbf{u}}^{\left(1\right)};\mathbf{f}\right)/\partial \mathbf{\psi }\right]}^{†}\end{array}\right).(44)$

5. Implement the boundary conditions represented by Eq. 32 into Eq. 41 and eliminate the remaining unknown boundary-values of the function ${\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)$ from the expression of the bilinear concomitant ${\left\{P^{\left(1\right)}\left[{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}$ by selecting appropriate boundary conditions for the function ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$, to ensure that Eq. 43 is well-posed while being independent of unknown values of ${\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)$ and of $\delta \mathbf{f}$. The boundary conditions thus chosen for the function ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$ can be represented in operator form as follows:

${\left\{{\mathbf{b}}_A^{\left(1\right)}\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}=\mathbf{0},\mathbf{x}\in \partial \mathrm{\Omega }\left[\mathbf{\lambda }\left({\mathbf{\alpha }}^0\right);\mathbf{\omega }\left({\mathbf{\alpha }}^0\right)\right].(45)$

The selection of the boundary conditions for ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$ represented by Eq. 45 eliminates the appearance of the unknown values of ${\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)$ in ${\left\{P^{\left(1\right)}\left[{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}$ and reduces this bilinear concomitant to a residual quantity containing boundary terms which involve only known values of ${\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right)$, ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)$, $\mathbf{f}$, and $\delta \mathbf{f}$. This residual quantity will be denoted as ${\left\{{\widehat{P}}^{\left(1\right)}\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}$. In general, this residual quantity does not automatically vanish, although it may do so occasionally.

6. The system of equations comprising Eq. 43 together with the boundary conditions represented Eq. 45 will be called the 1st-Level Adjoint Sensitivity System (1st-LASS). The solution ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$ of the 1st-LASS will be called the 1st-level adjoint sensitivity function. The 1st-LASS is called “first-level” (as opposed to “first-order”) because it does not contain any differential or functional-derivatives, but its solution, ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)$, will be used below to compute the first-order sensitivities of the response with respect to the components of the feature function $\mathbf{f}\left(\mathbf{\alpha }\right)$.

7. Using Eq. 40 together with the forward and adjoint boundary conditions represented by Eqs 32, 45 in Eq. 41 reduces the latter to the following relation:

$\begin{array}{l}{\left\{{〈{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right),{\mathbf{q}}_V^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]〉}_1\right\}}_{{\mathbf{\alpha }}^0}={\left\{{〈{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right),{\mathbf{A}}^{\left(1\right)}\left[2\times 2;\mathbf{x};\mathbf{f}\right]{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)〉}_1\right\}}_{{\mathbf{\alpha }}^0}\\ +{\left\{{\widehat{P}}^{\left(1\right)}\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}.\end{array}(46)$

8. In view of Eqs 28, 43, the first term on the right-side of Eq. 46 represents the indirect-effect term ${\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};{\mathbf{v}}^{\left(1\right)}\right]\right\}}_{ind}$. It therefore follows from Eq. 46 that the indirect-effect term can be expressed in terms of the 1st-level adjoint sensitivity function ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$ as follows:

$\begin{array}{l}{\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};{\mathbf{v}}^{\left(1\right)}\left(2;\mathbf{x}\right)\right]\right\}}_{ind}={\left\{{〈{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right),{\mathbf{q}}_V^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]〉}_1\right\}}_{{\mathbf{\alpha }}^0}\\ -{\left\{{\widehat{P}}^{\left(1\right)}\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}\equiv {\left\{\delta R\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{ind}.\end{array}(47)$

As indicated by the identity shown in Eq. 47, the variations $\delta \mathbf{\phi }$ and $\delta \mathbf{\psi }$ have been eliminated from the original expression of the indirect-effect term, which now depends on the 1st-level adjoint sensitivity function ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$. Adding the expression obtained in Eq. 47 with the expression for the direct-effect term defined in Eq. 26 yields, according to Eq. 25 the following expression for the total 1st-order sensitivity ${\left\{\delta R\left(\mathbf{\phi },\mathbf{\psi },\mathbf{f};\delta \mathbf{\phi },\delta \mathbf{\psi },\delta \mathbf{f}\right)\right\}}_{{\mathbf{\alpha }}^0}$ of the response $R\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{f}\right]$ with respect to the components of the feature function $\mathbf{f}\left(\mathbf{\alpha }\right)$:

$\begin{array}{l}{\left\{\delta R\left(\mathbf{\phi },\mathbf{\psi },\mathbf{f};\delta \mathbf{\phi },\delta \mathbf{\psi },\delta \mathbf{f}\right)\right\}}_{{\mathbf{\alpha }}^0}={\left\{\frac{\partial \mathbf{R}\left({\mathbf{u}}^{\left(1\right)};\mathbf{f}\right)}{\partial \mathbf{f}}\delta \mathbf{f}\right\}}_{{\mathbf{\alpha }}^0}+{\left\{{〈{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right),{\mathbf{q}}_V^{\left(1\right)}\left[2;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]〉}_1\right\}}_{{\mathbf{\alpha }}^0}\\ -{\left\{{\widehat{P}}^{\left(1\right)}\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}\equiv {\displaystyle \sum _{j_1=1}^{TF}}{\left\{R^{\left(1\right)}\left[j_1;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]\delta f_{j_1}\right\}}_{{\mathbf{\alpha }}^0}.\end{array}(48)$

The identity which appears in Eq. 48 emphasizes the fact that the variations $\delta \mathbf{\phi }$ and $\delta \mathbf{\psi }$, which are expensive to compute, have been eliminated from the final expressions of the 1st-order sensitivities $R^{\left(1\right)}\left[j_1;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ of the response with respect to the components $f_{j_1}\left(\mathbf{\alpha }\right),j_1=1,...,TF,$ of the “features function” $\mathbf{f}\left(\mathbf{\alpha }\right)$. The dependence on the variations $\delta \mathbf{\phi }$ and $\delta \mathbf{\psi }$ has been replaced in the expression of $R^{\left(1\right)}\left[j_1;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ by the dependence on the 1st-level adjoint sensitivity function ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$. It is very important to note that the 1st-LASS is independent of variations $\delta \mathbf{f}\left(\mathbf{\alpha }\right)$ in the components of the feature function and is consequently also independent of any variations $\delta \mathbf{\alpha }$ in the primary model parameters. Hence, the 1st-LASS needs to be solved only once to obtain the 1st-level adjoint sensitivity function ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$. Subsequently, the “indirect-effect term” is computed efficiently and exactly by simply performing the integrations required to compute the inner product over the adjoint function ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$, as indicated on the right-side of Eq. 48. Solving the 1st-Level Adjoint Sensitivity System (1st-LASS) requires the same computational effort as solving the original coupled linear system, entailing the following operations: (i) inverting (i.e., solving): the left-side of the original adjoint equation with the source ${\left[\partial S\left({\mathbf{u}}^{\left(1\right)};\mathbf{\alpha }\right)/\partial \mathbf{\phi }\right]}^{†}$ to obtain the 1st-level adjoint sensitivity function ${\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right)$; and (ii) inverting the left-side of the original forward equation with the source ${\left[\partial S\left({\mathbf{u}}^{\left(1\right)};\mathbf{\alpha }\right)/\partial \mathbf{\psi }\right]}^{†}$ to obtain the 1st-level adjoint sensitivity function ${\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)$.

The 1st-order sensitivities $R^{\left(1\right)}\left[j_1;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$, $j_1=1,...,TF$, can be expressed as an integral over the independent variables as follows:

$R^{\left(1\right)}\left[j_1;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]\triangleq {\displaystyle \underset{{\lambda }_1\left(\mathbf{\alpha }\right)}{\overset{{\omega }_1\left(\mathbf{\alpha }\right)}{\int }}}d{x}_1...{\displaystyle \underset{{\lambda }_{TI}\left(\mathbf{\alpha }\right)}{\overset{{\omega }_{TI}\left(\mathbf{\alpha }\right)}{\int }}}d{x}_{TI}{S}^{\left(1\right)}\left[j_1;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right].(49)$

In particular, if the residual bilinear concomitant is non-zero, the functions $S^{\left(1\right)}\left[j_1;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ would contain suitably defined Dirac delta-functionals for expressing the respective non-zero boundary terms as volume-integrals over the phase-space of the independent variables. Dirac-delta functionals would also be used in the expression of $S^{\left(1\right)}\left[j_1;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ to represent terms containing the derivatives of the boundary end-points with respect to the model and/or response parameters.

The response sensitivities with respect to the primary model parameters would be obtained by using the expression obtained in Eq. 49 in conjunction with the “chain rule” of differentiation provided in Eq. 11.

It is important to compare the results produced by the 1st-FASAM-L (for obtaining the sensitivities of the model response with respect to the model’s features) with the results produced by the 1st-CASAM methodology (the 1st-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems), which provides the expressions of the responses sensitivities directly with respect to the model’s primary parameters. Recall that the 1st-CASAM-L (Cacuci, 2022) yields the following expression for the 1st-order sensitivities of the response with respect to the primary model parameters:

$\begin{array}{l}{\left\{\frac{\partial R\left[j_1;{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{\alpha }\right]}{\partial {\alpha }_{j_1}}\right\}}_{{\mathbf{\alpha }}^0}={\left\{{\displaystyle \underset{{\lambda }_1\left(\mathbf{\alpha }\right)}{\overset{{\omega }_1\left(\mathbf{\alpha }\right)}{\int }}}d{x}_1...{\displaystyle \underset{{\lambda }_{TI}\left(\mathbf{\alpha }\right)}{\overset{{\omega }_{TI}\left(\mathbf{\alpha }\right)}{\int }}}d{x}_{TI}\frac{\partial S\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{\alpha }\right]}{\partial {\alpha }_{j_1}}\right\}}_{{\mathbf{\alpha }}^0}\\ +{\displaystyle \sum _{k=1}^{TI}}{\displaystyle \prod _{m=1,k\ne j}^{TI}}{\left\{{\displaystyle \underset{{\lambda }_m\left(\mathbf{\alpha }\right)}{\overset{{\omega }_m\left(\mathbf{\alpha }\right)}{\int }}}d{x}_{m}S\left[{\mathbf{u}}^{\left(1\right)}\left(2;...,{\omega }_k,...\right);\mathbf{\alpha }\right]\frac{\partial {\omega }_k\left(\mathbf{\alpha }\right)}{\partial {\alpha }_{j_1}}-S\left[{\mathbf{u}}^{\left(1\right)}\left(2;...,{\lambda }_k,...\right);\mathbf{\alpha }\right]\frac{\partial {\lambda }_k\left(\mathbf{\alpha }\right)}{\partial {\alpha }_{j_1}}\right\}}_{{\mathbf{\alpha }}^0}\\ +{\left\{{〈{\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right),\frac{\partial }{\partial {\alpha }_{j_1}}{\mathbf{q}}^{\left(1\right)}\left[{\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right);\mathbf{\alpha }\right]〉}_1\right\}}_{{\mathbf{\alpha }}^0}-{\left\{\frac{\partial }{\partial {\alpha }_{j_1}}{\widehat{P}}^{\left(1\right)}\left[{\mathbf{u}}^{\left(1\right)};{\mathbf{a}}^{\left(1\right)};\mathbf{\alpha }\right]\right\}}_{{\mathbf{\alpha }}^0};j_1=1,...,TP.\end{array}(50)$

The same 1st-level adjoint function ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)$ appears in Eq. 50 as well as in Eq. 49. Therefore, a single “large-scale computation” (needed to solve the 1st-LASS to determine the 1st-level adjoint function) is required for obtaining either the response sensitivities with respect to the components, $f_j\left(\mathbf{\alpha }\right)$, $j=1,...,TF$, of the feature function $\mathbf{f}\left(\mathbf{\alpha }\right)$ using the 1st-FASAM-L, or for obtaining the response sensitivities directly with respect to the primary model parameters ${\alpha }_j$, $j=1,...,TP$, using the 1st-CASAM-L. On the other hand, the use of the 1st-CASAM-L would require performing a number of $TP$ integrations to compute all of the response sensitivities with respect to the primary parameters, but the 1st-FASAM-L would require only $TF$ integrations ($TF\ll TP$) to compute all of the response sensitivities with respect to the components $f_j\left(\mathbf{\alpha }\right)$ of the feature function. Hence, the 1st-FASAM-L is more efficient than the 1st-CASAM-L, so the 1st-FASAM-L is the most efficient method for computing the exact expressions of the first-order sensitivities of a generic model response of the form $R\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right);\mathbf{\alpha }\right]$ with respect to the components of the “features” function $\mathbf{f}\left(\mathbf{\alpha }\right)$, and subsequently with respect to the primary model parameters. As will be shown in the sequel, the computational savings provided by the nth-FASAM-L increase massively by comparison to the nth-CASAM-L (or any other method) as the order “n” of the computed sensitivities increases.

The expression obtained in Eq. 48 is the same as the particular form taken on by general expression provided in Eq. 13 for $n=1$, where:

(i) the 1st-level forward/adjoint function ${\mathbf{u}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[\mathbf{\phi }\left(\mathbf{x}\right),\mathbf{\psi }\left(\mathbf{x}\right)\right]}^{†}$ is the solution of the 1st-LFAS defined by Eqs 19, 20, which has the same expression as the particular form taken on by the nth-LFAS, cf. Eqs 15, 16, for $n=1$;

(ii) the 1st-level adjoint sensitivity function ${\mathbf{a}}^{\left(1\right)}\left(2;\mathbf{x}\right)\triangleq {\left[{\mathbf{a}}_1^{\left(1\right)}\left(\mathbf{x}\right),{\mathbf{a}}_2^{\left(1\right)}\left(\mathbf{x}\right)\right]}^{†}$ is the solution of the 1st-LASS defined by Eqs 43, 45, which has the same expression as the particular form taken on by the nth-LFAS, cf. Eqs 17, 18, for $n=1$.

Thus, the first step in the “proof by mathematical induction” of the pattern underlying the nth-FASAM-L has been completed, having shown that this pattern holds for $n=1$.

2.4 Proving that the conjectured mathematical framework of the nth-FASAM-L methodology also holds for $n\to n+1$, i.e., for the (n + 1)^th-FASAM-L framework

The last step of the “proof by mathematical induction” to establish the validity the nth-FASAM-L framework is to show that the formalism assumed to be correct for the computation of the nth-order sensitivities also holds true for the computation of the (n + 1)^th-order sensitivities. This proof entails showing that the formulas obtained by computing the (n + 1)^th-order sensitivities using Eqs 14‒18 as the starting point will be the same as would be obtained by replacing “n” with “(n + 1)” in Eqs 14‒18.

The nth-order response sensitivity defined in Eq. 14 can be considered to be a function of the (n + 1)^th-level function ${\mathbf{u}}^{\left(n+1\right)}\left(2^{\left(n+1\right)};\mathbf{x}\right)\triangleq {\left[{\mathbf{u}}^{\left(n\right)}\left(2^n;\mathbf{x}\right),{\mathbf{a}}^{\left(n\right)}\left(2^n;\mathbf{x}\right)\right]}^{†}$, which is the solution of the (n + 1)^th-Level Forward/Adjoint System, abbreviated as “(n + 1)^th-LFAS”, which is obtained by concatenating Eqs 15‒18 and is written in the following form:

${\mathbf{F}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};\mathbf{f}\left(\mathbf{\alpha }\right)\right]{\mathbf{u}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right)={\mathbf{q}}_F^{\left(n+1\right)}\left[2^{n+1};{\mathbf{u}}^{\left(n\right)}\left(2^n;\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right];\mathbf{x}\in \mathrm{\Omega };(51)$

${\mathbf{b}}_F^{\left(n+1\right)}\left(2^{\left(n+1\right)};{\mathbf{u}}^{\left(n+1\right)};\mathbf{f}\right)\triangleq {\left({\mathbf{b}}_F^{\left(n\right)},{\mathbf{b}}_A^{\left(n\right)}\right)}^{†}=\mathbf{0};\mathbf{x}\in \partial \mathrm{\Omega }.(52)$

The following definitions were used in Eqs 51, 52, where the explicit dependence on the indices $j_k,k=1,...,n$, has been omitted, for simplicity:

${\mathbf{F}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};\mathbf{f}\left(\mathbf{\alpha }\right)\right]\triangleq diag\left({\mathbf{F}}^{\left(n\right)},{\mathbf{A}}^{\left(n\right)}\right);{\mathbf{u}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right)\triangleq {\left[{\mathbf{u}}^{\left(n\right)}\left(2^n;\mathbf{x}\right),{\mathbf{a}}^{\left(n\right)}\left(2^n;\mathbf{x}\right)\right]}^{†};(53)$

${\mathbf{q}}_F^{\left(n+1\right)}\left[2^{n+1};{\mathbf{u}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]\triangleq {\left[{\mathbf{q}}_F^{\left(n\right)}\left(2^n;\mathbf{x};\mathbf{f}\right),{\mathbf{q}}_A^{\left(n\right)}\left(2^n;{\mathbf{u}}^{\left(n\right)};\mathbf{f}\right)\right]}^{†}.(54)$

Next, it will be assumed that, for each index $j_1,...,j_n$, the 1st-order total G-differential of the nth-order sensitivities $R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right);\mathbf{f}\left(\mathbf{\alpha }\right)\right]$ exists and is linear in the variational functions ${\mathbf{v}}^{\left(n+1\right)}\left(2^{n+1};j_{n-1},...,j_1;\mathbf{x}\right)\triangleq {\left[{\mathbf{v}}^{\left(n\right)}\left(2^n;\mathbf{x}\right),\delta {\mathbf{a}}^{\left(n\right)}\left(2^n;\mathbf{x}\right)\right]}^{†}$ and $\delta \mathbf{f}$ in a neighborhood around the nominal values of the respective state functions and components of the feature function. In this case, the 1st-order total G-differential of $R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)};\mathbf{f}\right]$ is by definition obtained as follows:

$\begin{array}{l}{\left\{\delta R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)};\mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}\triangleq {\left\{\frac{d}{d\epsilon }{R}^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)}+\epsilon {\mathbf{v}}^{\left(n+1\right)};\mathbf{f}+\epsilon \delta \mathbf{f}\right]\right\}}_{\epsilon =0}\\ \triangleq {\displaystyle \sum _{j_{n+1}=1}^{TF}}{\left\{\frac{\partial R^{\left(n\right)}\left[...;{\mathbf{u}}^{\left(n+1\right)};\mathbf{f}\right]}{\partial f_{j_{n+1}}}\right\}}_{{\mathbf{\alpha }}^0}\delta f_{j_{n+1}}+{\left\{\delta R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)};{\mathbf{v}}^{\left(n+1\right)};\mathbf{f}\right]\right\}}_{ind},\end{array}(55)$

where the quantity ${\left\{\delta R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)};{\mathbf{v}}^{\left(n+1\right)};\mathbf{f}\right]\right\}}_{ind}$ denotes the “indirect-effect term” and is defined as follows:

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

The vector ${\mathbf{v}}^{\left(n+1\right)}\left(2^{n+1};j_{n-1},...,j_1;\mathbf{x}\right)$, which is needed to evaluate the indirect-effect term ${\left\{\delta R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)};{\mathbf{v}}^{\left(n+1\right)};\mathbf{f}\right]\right\}}_{ind}$, is the solution of the (n + 1)^th-Level Variational Sensitivity System, abbreviated as (n + 1)^th-LVSS, which is obtained by taking the (first-order) G-differential of the (n + 1)^th-LFAS defined by Eqs 53, 54. Performing this G-differentiation yields the following relations which define the (n + 1)^th-LVSS:

$\begin{array}{l}{\left\{\frac{d}{d\epsilon }{\mathbf{F}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};{\mathbf{f}}^0+\epsilon \delta \mathbf{f}\right]\left[{\mathbf{u}}^{\left(n+1,0\right)}\left(2^{n+1};\mathbf{x}\right)+\epsilon {\mathbf{v}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right)\right]\right\}}_{\epsilon =0}\\ ={\left\{\frac{d}{d\epsilon }{\mathbf{q}}_F^{\left(n+1\right)}\left[2^{n+1};{\mathbf{u}}^{\left(n,0\right)}\left(2^n;\mathbf{x}\right)+\epsilon {\mathbf{v}}^{\left(n\right)}\left(2^n;\mathbf{x}\right);{\mathbf{f}}^0+\epsilon \delta \mathbf{f}\right]\right\}}_{\epsilon =0};\mathbf{x}\in \mathrm{\Omega };\end{array}(57)$

${\left\{\frac{d}{d\epsilon }{\mathbf{b}}_F^{\left(n+1\right)}\left[2^{n+1};{\mathbf{u}}^{\left(n+1,0\right)}\left(2^{n+1};\mathbf{x}\right)+\epsilon {\mathbf{v}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right);{\mathbf{f}}^0+\epsilon \delta \mathbf{f}\right]\right\}}_{\epsilon =0}=\mathbf{0};\mathbf{x}\in \partial \mathrm{\Omega };(58)$

Carrying out the differentiation with respect to $\epsilon$ in Eqs 57, 58, and setting $\epsilon =0$ in the resulting expressions yields the following (n+1)^th-LVSS for the (n + 1)^th-level variational function ${\mathbf{v}}^{\left(n+1\right)}\left(2^{n+1};j_{n-1},...,j_1;\mathbf{x}\right)$:

${\left\{{\mathbf{V}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};\mathbf{x};\mathbf{f}\right]{\mathbf{v}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right)\right\}}_{{\mathbf{\alpha }}^0}={\left\{{\mathbf{q}}_V^{\left(n+1\right)}\left[2^{n+1};{\mathbf{u}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0};\mathbf{x}\in \mathrm{\Omega };(59)$

${\left\{{\mathbf{b}}_v^{\left(n+1\right)}\left({\mathbf{u}}^{\left(n+1\right)};{\mathbf{v}}^{\left(n+1\right)};\mathbf{f};\delta \mathbf{f}\right)\right\}}_{{\mathbf{\alpha }}^0}=\mathbf{0};\mathbf{x}\in \partial \mathrm{\Omega }.(60)$

Solving the (n + 1)^th-LVSS is prohibitive computationally. Therefore, the need for solving the (n + 1)^th-LVSS will be avoided by expressing the indirect-effect term ${\left\{\delta R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)};{\mathbf{v}}^{\left(n+1\right)};\mathbf{f}\right]\right\}}_{ind}$ in an alternative way, which eliminates the appearance of the variational function ${\mathbf{v}}^{\left(n+1\right)}\left(2^{n+1};j_{n-1},...,j_1;\mathbf{x}\right)$ by replacing it with the solution of the “(n + 1)^th-Level Adjoint Sensitivity System,” abbreviated as (n + 1)^th-LASS). This (n + 1)^th-LASS will be constructed below by implementing the same sequence of logical steps as were followed when constructing the first- (and lower-) level adjoint sensitivity systems, namely,:

(i) The (n + 1)^th-LASS is constructed in a Hilbert space, denoted as ${\mathrm{H}}_{n+1}$, comprising block-vectors of the form ${\mathbf{\chi }}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right)\in {\mathrm{H}}_{n+1}$, ${\mathbf{\chi }}^{\left(n+1\right)}\left(\mathbf{x}\right)\triangleq {\left[...,{\mathbf{\chi }}_k^{\left(n+1\right)}\left(\mathbf{x}\right),...\right]}^{†}$, for $k=1,...,2^{n+1}$, each of these comprising elements having the following structure: ${\mathbf{\chi }}_k^{\left(n+1\right)}\left(\mathbf{x}\right)\triangleq {\left[{\chi }_{k,1}^{\left(n+1\right)}\left(\mathbf{x}\right),...,{\chi }_{k,TD}^{\left(n+1\right)}\left(\mathbf{x}\right)\right]}^{†}$. The inner product between two elements, ${\mathbf{\chi }}^{\left(n+1\right)}\left(\mathbf{x}\right)\in {\mathrm{H}}_{n+1}$ and ${\mathbf{\theta }}^{\left(n+1\right)}\left(\mathbf{x}\right)\in {\mathrm{H}}_{n+1}$, of the Hilbert space ${\mathrm{H}}_{n+1}$, will be denoted as ${〈{\mathbf{\chi }}^{\left(n+1\right)}\left(\mathbf{x}\right),{\mathbf{\theta }}^{\left(n+1\right)}\left(\mathbf{x}\right)〉}_{n+1}$ and is defined as follows:

${〈{\mathbf{\chi }}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right),{\mathbf{\theta }}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right)〉}_{\left(n+1\right)}\triangleq {\displaystyle \sum _{i=1}^{2^{n+1}}}{〈{\mathbf{\chi }}_i^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right),{\mathbf{\theta }}_i^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right)〉}_0.(61)$

(ii) Using the definition provided in Eq. 61, form the inner product in ${\mathrm{H}}_{n+1}$ of Eq. 59 with a yet undefined vector-valued function ${\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right)\triangleq {\left[...,{\mathbf{a}}_k^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right),...\right]}^{†}\in {\mathrm{H}}_{n+1};k=1,...,2^{n+1}$, $j_1=1,...,TF$, $j_2=1,...,j_1$; $j_{n+1}=1,...,j_n$, to obtain the following relation:

$\begin{array}{l}{\left\{{〈{\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right),{\mathbf{V}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};\mathbf{x};\mathbf{f}\right]{\mathbf{v}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right)〉}_{n+1}\right\}}_{{\mathbf{\alpha }}^0}\\ =\left\{\langle {\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right)\right.{\left.{,{\mathbf{q}}_V^{\left(n+1\right)}\left[2^{n+1};{\mathbf{u}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]〉}_{n+1}\right\}}_{{\mathbf{\alpha }}^0}\\ ={\left\{{〈{\mathbf{v}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right),{\mathbf{A}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};\mathbf{x};\mathbf{f}\right]{\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right)〉}_{n+1}\right\}}_{{\mathbf{\alpha }}^0}\\ +{\left\{P^{\left(n+1\right)}\left[{\mathbf{v}}^{\left(n+1\right)};{\mathbf{a}}^{\left(n+1\right)};\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0},\end{array}(62)$

where ${\left\{P^{\left(n+1\right)}\left[{\mathbf{v}}^{\left(n+1\right)};{\mathbf{a}}^{\left(n+1\right)};\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}$ denotes the bilinear concomitant defined on the phase-space boundary $\mathbf{x}\in \partial \mathrm{\Omega }$, evaluated at the nominal values of the model parameter and respective functions, and where ${\mathbf{A}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};\mathbf{x};\mathbf{f}\right]$ is the formal adjoint of the matrix-valued operator ${\mathbf{V}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};\mathbf{x};\mathbf{f}\right]$ as defined in Eq. 63, below:

${\mathbf{A}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};\mathbf{x};\mathbf{f}\right]\triangleq {\left\{{\mathbf{V}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};\mathbf{x};\mathbf{f}\right]\right\}}^{*}.(63)$

(iii) The first term on right-side of the second equality in Eq. 62 is now required to represent the indirect-effect term ${\left\{\delta R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)};{\mathbf{v}}^{\left(n+1\right)};\mathbf{f}\right]\right\}}_{ind}$. This is achieved by requiring that the (n+1)^th-level adjoint sensitivity function ${\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right)\triangleq {\left[...,{\mathbf{a}}_k^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right),...\right]}^{†}\in {\mathrm{H}}_{n+1};k=1,...,2^{n+1}$, be the solution of the following (n+1)^th-LASS defined in Eqs. 64, 65, below, for $j_1=1,...,TP$; $j_2=1,...,j_1$;… $j_n=1,...,j_{n-1}$:

${\mathbf{A}}^{\left(n+1\right)}\left[2^{n+1}\times 2^{n+1};\mathbf{x};\mathbf{f}\right]{\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right)={\mathbf{s}}_A^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{f}\right),(64)$

${\left\{{\mathbf{b}}_A^{\left(n+1\right)}\left[{\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right);{\mathbf{u}}^{\left(n+1\right)}\left(j_{n-1},...,j_1;\mathbf{x}\right);\mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}=\mathbf{0},\mathbf{x}\in \partial \mathrm{\Omega },(65)$

where the vector ${\mathbf{s}}_A^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{f}\right)\triangleq {\left[...,{\mathbf{s}}_k^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{f}\right),...\right]}^{†}$, $k=1,...,2^{n+1}$, comprises $2^{n+1}$ components defined in Eq. 66, below, for each $j_1=1,...,TP$; $j_2=1,...,j_1$;…; $j_n=1,...,j_{n-1}$:

${\mathbf{s}}_A^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{f}\right)\triangleq \partial S^{\left(n\right)}/\partial {\mathbf{u}}^{\left(n+1\right)}\left(\mathbf{x}\right).(66)$

(iv) The (n + 1)^th-level adjoint boundary conditions represented by Eq. 65 are selected so as to eliminate, in conjunction with the boundary conditions represented by Eq. 60, all of the unknown values of the functions ${\mathbf{v}}^{\left(n+1\right)}\left(2^{n+1};j_{n-1},...,j_1;\mathbf{x}\right)$ in the expression of the bilinear concomitant ${\left\{P^{\left(n+1\right)}\left[{\mathbf{v}}^{\left(n+1\right)};{\mathbf{a}}^{\left(n+1\right)};\mathbf{f};\delta \mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}$. This bilinear concomitant may vanish after implementing the boundary conditions represented by Eqs 60, 65. However, if it does not vanish, this bilinear concomitant will be reduced to a residual quantity which will comprise only known values of ${\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right)$, ${\mathbf{u}}^{\left(n+1\right)}\left(j_{n-1}...j_1;\mathbf{x}\right)$, $\mathbf{f}\left(\mathbf{\alpha }\right)$ and $\delta \mathbf{f}\left(\mathbf{\alpha }\right)$, and which will be denoted as ${\left\{{\widehat{P}}^{\left(n+1\right)}\left({\mathbf{a}}^{\left(n+1\right)};{\mathbf{u}}^{\left(n+1\right)};\mathbf{f};\delta \mathbf{f}\right)\right\}}_{{\mathbf{\alpha }}^0}$.

(v) Using in Eq. 56 the equations underlying the (n + 1)^th-LASS together with the relation provided in Eq. 62 yields the following expression for the indirect-effect term ${\left\{\delta R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)};{\mathbf{v}}^{\left(n+1\right)};\mathbf{f}\right]\right\}}_{ind}$ in terms of the (n + 1)^th-level adjoint sensitivity functions ${\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right)$, for each $j_1=1,...,TP$; $j_2=1,...,j_1$;…; $j_n=1,...,j_{n-1}$:

$\begin{array}{l}{\left\{\delta R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)};{\mathbf{v}}^{\left(n+1\right)};\mathbf{f}\right]\right\}}_{ind}=-{\left\{{\widehat{P}}^{\left(n+1\right)}\left({\mathbf{a}}^{\left(n+1\right)};{\mathbf{u}}^{\left(n+1\right)};\mathbf{f};\delta \mathbf{f}\right)\right\}}_{{\mathbf{\alpha }}^0}\\ +{\left\{{〈{\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right),{\mathbf{q}}_V^{\left(n+1\right)}\left[2^{n+1};{\mathbf{u}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]〉}_{n+1}\right\}}_{{\mathbf{\alpha }}^0}.\end{array}(67)$

Adding the result obtained in Eq. 67 for the indirect effect term to the result provided in Eq. 55 for the direct effect term yields the following expression for the total nth-order G-variation of the response:

$\begin{array}{l}{\left\{\delta R^{\left(n\right)}\left[j_n;...;j_1;{\mathbf{u}}^{\left(n+1\right)};\mathbf{f}\right]\right\}}_{{\mathbf{\alpha }}^0}={\displaystyle \sum _{j_{n+1}=1}^{TF}}{\left\{\frac{\partial }{\partial f_{j_{n+1}}}{\displaystyle \underset{{\lambda }_1\left(\mathbf{\alpha }\right)}{\overset{{\omega }_1\left(\mathbf{\alpha }\right)}{\int }}}d{x}_1...{\displaystyle \underset{{\lambda }_{TI}\left(\mathbf{\alpha }\right)}{\overset{{\omega }_{TI}\left(\mathbf{\alpha }\right)}{\int }}}d{x}_{TI}{S}^{\left(n\right)}\left(j_n,...,j_1;{\mathbf{u}}^{\left(n+1\right)};\mathbf{\alpha }\right)\delta f_{j_{n+1}}\right\}}_{{\mathbf{\alpha }}^0}\\ +{\left\{{〈{\mathbf{a}}^{\left(n+1\right)}\left(j_n,...,j_1;\mathbf{x}\right),{\mathbf{q}}_V^{\left(n+1\right)}\left[2^{n+1};{\mathbf{u}}^{\left(n+1\right)}\left(2^{n+1};\mathbf{x}\right);\mathbf{f};\delta \mathbf{f}\right]〉}_{n+1}\right\}}_{{\mathbf{\alpha }}^0}\\ -{\left\{{\widehat{P}}^{\left(n+1\right)}\left({\mathbf{a}}^{\left(n+1\right)};{\mathbf{u}}^{\left(n+1\right)};\mathbf{f};\delta \mathbf{f}\right)\right\}}_{{\mathbf{\alpha }}^0}\equiv {\displaystyle \sum _{j_{n+1}=1}^{TF}}{\left\{R^{\left(n+1\right)}\left(j_{n+1},...,j_1;{\mathbf{u}}^{\left(n+1\right)};{\mathbf{a}}^{\left(n+1\right)};\mathbf{f}\right)\right\}}_{\left({\mathbf{\alpha }}^0\right)}\delta f_{j_{n+1}},\end{array}(68)$

where $R^{\left(n+1\right)}\left(j_{n+1},...,j_1;{\mathbf{u}}^{\left(n+1\right)};{\mathbf{a}}^{\left(n+1\right)};\mathbf{f}\right)$ denotes the (n + 1)^th-order partial sensitivity of the response $R\left[{\mathbf{u}}^{\left(1\right)}\left(\mathbf{x}\right);\mathbf{\alpha }\right]$ with respect to the components of the feature function $\mathbf{f}\left(\mathbf{\alpha }\right)$, evaluated at the nominal parameter values ${\mathbf{\alpha }}^0$.

The result obtained in Eq. 68 for the expression of the (n + 1)^th-order sensitivity, which was obtained by determining the first-order differential of the nth-order sensitivity, is identical to the expression that would be obtained by advancing the index, from n to (n + 1), in the expression of the nth-order sensitivity that was conjectured in Eq. 13. Thus, the proof by mathematical induction of the general mathematical framework underlying the nth-CASAM-L is thereby completed.

The essential characteristics of the nth-FASAM-L methodology are tabularized in Tables 1‒4, below, to underscore the conceptual parallelism between the nth-FASAM-L and the nth-CASAM-L (Cacuci, 2022) methodologies.

Table 1

Table 1. 1st-FASAM-L: 1st-order (n = 1) sensitivities of response to model features.

Table 2

Table 2. 2nd-FASAM-L: 2nd-order (n = 2) sensitivities of response to model features.

Table 3

Table 3. nth-FASAM-L: nth-order sensitivities of response to model features.

Table 4

Table 4. (n + 1)^th-FASAM-L: (n + 1)^th-order sensitivities of response to model features.

An overview, in tabular form, of the computational frameworks of the nth-CASAM-L, nth-CASAM-N, nth-FASAM-L, and nth-FASAM-N methodologies, highlighting their objectives, characteristics, and interrelationships is presented in Table 5, below.

Table 5

Table 5. The nth-CASAM-L, nth-CASAM-N, nth-FASAM-L, nth-FASAM-N methodologies: main features.

Formally, the results produced by the nth-FASAM-L can be written in the same mathematical forms as those produced by the nth-CASAM-L, with the fundamental difference that the number of large-scale computations needed within the nth-FASAM-L is dictated by the number $TF$ of “feature function components” whereas the number of large-scale computations needed within the nth-CASAM-L is dictated by the number $TP$ of primary model parameters. In particular, a single large-scale adjoint computation is needed to solve the 1st-LASS (which is the same for both the 1st-FASAM-L and the 1st-CASAM-L) to obtain the first-order sensitivities with respect to the model parameters. Obtaining the second-order sensitivities of the response with respect to the primary model parameters requires at most $TP\left(TP+1\right)/2$ large-scale computations (to solve the 2nd-LASS) within the 2nd-CASAM-L. Obtaining the same second-order sensitivities using the 2nd-FASAM-L requires at most $TF\left(TF+1\right)/2$ large-scale computations (to solve the 2nd-LASS) followed by analytical derivations to obtain the second-order sensitivities with respect to the model parameters from the second-order sensitivities with respect to the components of the feature function produced by the 2nd-FASAM-L. The same parallel holds for the computation of all of the higher-order sensitivities: the computation of the 3rd-order sensitivities with respect to the primary model parameters requires at most $TP\left(TP+1\right)\left(TP+2\right)/3!$ computations if using the 3rd-CASAM-L, as opposed to at most $TF\left(TF+1\right)\left(TF+2\right)/3!$ large-scale computations plus analytical derivations if using the 3rd-CASAM-L. The computation of the 4th-order sensitivities with respect to the primary model parameters requires at most $TP\left(TP+1\right)\left(TP+2\right)\left(TP+3\right)/4!$ computations if using the 4th-CASAM-L, as opposed to at most $TF\left(TF+1\right)\left(TF+2\right)\left(TF+3\right)/4!$ large-scale computations plus analytical derivations if using the 4th-CASAM-L; and so on. Since $TF\ll TP$, it is evident that the nth-FASAM-L methodology becomes increasingly more efficient than the nth-CASAM-L as the order of computed sensitivities increases.

3 Concluding discussion

This work has presented the “nth-Order Feature Adjoint Sensitivity Analysis Methodology for Response-Coupled Forward/Adjoint Linear Systems” (abbreviated as “nth-FASAM-L”), which is the most efficient methodology for computing exact expressions of sensitivities of model responses to features of model parameters and, subsequently, to the model parameters themselves for such linear systems. This efficiency stems from the maximal reduction of the number of adjoint computations (which are “large-scale” computations), by comparison to the extant high-order adjoint sensitivity analysis methodology nth-CASAM-N (the “nth-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems”). Specific details are as follows:

(i) Comparing the mathematical framework of the nth-FASAM-N methodology to the framework of the nth-CASAM-N methodology indicates that the components $f_i\left(\mathbf{\alpha }\right),i=1,...,TF$, of the “feature function” $\mathbf{f}\left(\mathbf{\alpha }\right)\triangleq {\left[f_1\left(\mathbf{\alpha }\right),...,f_{TF}\left(\mathbf{\alpha }\right)\right]}^{†}$ play within the nth-FASAM-N the same role as played by the components ${\alpha }_j,j=1,...,TP$, of the “vector of primary model parameters” $\mathbf{\alpha }\triangleq {\left({\alpha }_1,\dots ,{\alpha }_{TP}\right)}^{†}$ within the framework of the nt^h-CASAM-N. It is important to note that the total number of model parameters is always larger (usually by wide margin) than the total number of components of the feature function $\mathbf{f}\left(\mathbf{\alpha }\right)$, i.e., $TP\gg TF$.

(ii) The 1st-FASAM-N and the 1st-CASAM-N methodologies require a single large-scale “adjoint” computations for solving the 1st-LASS (1st-Level Adjoint Sensitivity System), so they are comparably efficient for computing the exact expressions of the first-order sensitivities of a model response to the model’s uncertain parameters, boundaries, and internal interfaces.

(iii) For computing the exact expressions of the second-order response sensitivities with respect to the primary model’s parameters, the 2nd-FASAM-N methodology requires, at most, as many large-scale “adjoint” computations as there are “feature functions of parameters” $f_i\left(\mathbf{\alpha }\right),i=1,...,TF$ (where $TF$ denotes the total number of feature functions) for solving the left-side of the 2nd-LASS with $TF$ distinct sources on its right-side. By comparison, the 2nd-CASAM-N methodology requires at most $TP$ (where $TP$ denotes the total number of model parameters) large-scale computations for solving the same left-side of the 2nd-LASS but with $TP$ distinct sources. Since $TF\ll TP$, the 2nd-FASAM-N methodology is considerably more efficient than the 2nd-CASAM-N methodology for computing the exact expressions of the second-order sensitivities of a model response to the model’s uncertain parameters, boundaries, and internal interfaces.

(iv) For computing the exact expressions of the third-order response sensitivities with respect to the primary model’s parameters, the 3rd-FASAM-N requires at most $TF\left(TF+1\right)/2$ large-scale “adjoint” computations for solving the 3rd-LASS with $TF\left(TF+1\right)/2$ distinct sources, while the 3rd-CASAM-N methodology requires at most $TP\left(TP+1\right)/2$ large-scale computations for solving the 3rd-LASS with $TP\left(TP+1\right)/2$ distinct sources. The same computational-count of ”large-scale computations” caries over when computing the higher-order sensitivities, i.e., the formula for calculating the “number of large-scale adjoint computations” is formally the same for both the nth-FASAM-N and the nth-CASAM-N methodologies, but the “variable” in the formula for determining the number of adjoint computations for the nth-FASAM-N methodology is $TF$ (i.e., total number of feature functions) while the counterpart for the formula for determining the number of adjoint computations for the nth-CASAM-N is methodology is $TP$ (i.e., total number of model parameters). Since $TF\ll TP$, it follows that the higher the order of computed sensitivities, the mode efficient the nth-FASAM-N methodology becomes by comparison to the nth-CASAM-N methodology.

(v) When a model has no “feature” functions of parameters, but only comprises primary parameters, the nth-FASAM-N methodology becomes identical to the nth-CASAM-N methodology.

(vi) Both the nth-FASAM-N and the nth-CASAM-N methodologies are formulated in linearly increasing higher-dimensional Hilbert spaces −as opposed to exponentially increasing parameter-dimensional spaces− thus overcoming the curse of dimensionality in sensitivity analysis of nonlinear systems. Both the nth-FASAM-N and the nth-CASAM-N methodologies are incomparably more efficient and more accurate than any other methods (statistical, finite differences, etc.) for computing exact expressions of response sensitivities (of any order) with respect to the model’s uncertain parameters, boundaries, and internal interfaces.

The question of “when to stop computing progressively higher-order sensitivities ?” has been addressed by Cacuci (2022), Cacuci D. G. (2023) in conjunction with the question of convergence of the Taylor-series expansion of the response in terms of the uncertain model parameters, cf; Eqs 10, 12. These Taylor-series expansions provide the fundamental premise, even if not explicitly recognized, for obtaining the expressions provided by the “propagation of errors” methodology (as originally proposed by Tukey, 1957; and generalized by Cacuci, 2022) for the cumulants of the model response distribution in the phase-space of model parameters. The convergence of these Taylor-series, which depend on both the response sensitivities with respect to parameters and the uncertainties associated with the parameter distribution, must be ensured. This can be done by ensuring that the combination of parameter uncertainties and response sensitivities are sufficiently small to fall inside the respective radius of convergence of each of these Taylor-series expansions. The application of the nth-FASAM-N to a representative response-coupled forward/adjoint linear model stemming from the field of energy-dependent particle transport in a mixture of materials will be presented in the accompanying work designated as “Part II” (Cacuci, 2024c).

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

DC: Conceptualization, Methodology, Project administration, Validation, Writing–original draft, Writing–review and editing.

Funding

The author(s) declare that no financial support was received for the research, authorship, and/or publication of this article.

Conflict of interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Alcouffe, R. E., Baker, R. S., Dahl, J. A., Turner, S. A., and Ward, R. (2008). PARTISN: a time-dependent, parallel neutral particle transport code system. Los Alamos, NM, USA: Los Alamos National Laboratory. LA-UR-08-07258.

Google Scholar

Bellman, R. E. (1957). Dynamic programming. USA: Rand Corporation, Princeton University Press. Republished: Bellman, RE (2003) Dynamic Programming. Courier Dover Publications, ISBN 978-0-486-42809-3, USA.

Google Scholar

Cacuci, D. G. (1981). Sensitivity theory for nonlinear systems: I. Nonlinear functional analysis approach. J. Math. Phys. 22, 2794–2802. doi:10.1063/1.525186

CrossRef Full Text | Google Scholar

Cacuci, D. G. (2015). Second-order adjoint sensitivity analysis methodology (2nd-ASAM) for computing exactly and efficiently first- and second-order sensitivities in large-scale linear systems: I. Computational methodology. J. Comp. Phys. 284, 687–699. doi:10.1016/j.jcp.2014.12.042

CrossRef Full Text | Google Scholar

Cacuci, D. G. (2016). The second-order adjoint sensitivity analysis methodology for nonlinear systems—I: theory. Nucl. Sci. Eng. 184, 16–30. doi:10.13182/nse16-16

CrossRef Full Text | Google Scholar

Cacuci, D. G. (2022). The nth-order comprehensive adjoint sensitivity analysis methodology (nth-CASAM): overcoming the curse of dimensionality in sensitivity and uncertainty analysis, Volume I: linear systems. Cham, Switzerland: Springer Nature, 362. doi:10.1007/978-3-030-96364-4

CrossRef Full Text | Google Scholar

Cacuci, D. G. (2023a). The nth-order comprehensive adjoint sensitivity analysis methodology (nth-CASAM): overcoming the curse of dimensionality in sensitivity and uncertainty analysis, volume III: nonlinear systems. Cham, Switzerland: Springer Nature, 369. doi:10.1007/978-3-031-22757-8

CrossRef Full Text | Google Scholar

Cacuci, D. G. (2023b). Computation of high-order sensitivities of model responses to model parameters. II: introducing the second-order adjoint sensitivity analysis methodology for computing response sensitivities to functions/features of parameters. Energies 16, 6356. doi:10.3390/en16176356

CrossRef Full Text | Google Scholar

Cacuci, D. G. (2024a). Introducing the nth-order features adjoint sensitivity analysis methodology for nonlinear systems (nth-FASAM-N): I. Mathematical framework. Am. J. Comput. Math. 14, 11–42. doi:10.4236/ajcm.2024.141002

CrossRef Full Text | Google Scholar

Cacuci, D. G. (2024b). Introducing the nth-order features adjoint sensitivity analysis methodology for nonlinear systems (nth-FASAM-N): II. Illustrative example. Am. J. Comput. Math. 14, 43–95. doi:10.4236/ajcm.2024.141003

CrossRef Full Text | Google Scholar

Cacuci, D. G. (2024c). The nth-order features adjoint sensitivity analysis methodology for response-coupled forward/adjoint linear systems (nth-FASAM-L): I. mathematical framework. Front. Energy Res 12, 1417594. doi:10.3389/fenrg.2024.1417594

CrossRef Full Text | Google Scholar

Cacuci, D. G., and Fang, R. (2023). The nth-order comprehensive adjoint sensitivity analysis methodology (nth-CASAM): overcoming the curse of dimensionality in sensitivity and uncertainty analysis, volume II: application to a large-scale system. Nature Switzerland, Cham: Springer, 463. doi:10.1007/978-3-031-19635-5

CrossRef Full Text | Google Scholar

Chadwick, M. B., Herman, M., Obložinský, P., Dunn, M. E., Danon, Y., Kahler, A. C., et al. (2011). ENDF/B-VII.1: nuclear data for science and technology: cross sections, covariances, fission product yields and decay data. Nucl. Data Sheets 112, 2887–2996. doi:10.1016/j.nds.2011.11.002

CrossRef Full Text | Google Scholar

Conlin, J. L., Parsons, D. K., Gardiner, S. J., Gray, M., Lee, M. B., and White, M. C. (2013). MENDF71X: multigroup neutron cross-section data tables based upon ENDF/B-VII.1X; los alamos national laboratory report LA-UR-15-29571. Los Alamos, NM, USA: Los Alamos National Laboratory, 2013.

Google Scholar

Cukier, R. I., Levine, H. B., and Shuler, K. E. (1978). Nonlinear sensitivity analysis of multiparameter model systems. J. Comput. Phys. 26, 1–42. doi:10.1016/0021-9991(78)90097-9

CrossRef Full Text | Google Scholar

Dunker, A. M. (1984). The decoupled direct method for calculating sensitivity coefficients in chemical kinetics. J. Chem. Phys. 81, 2385–2393. doi:10.1063/1.447938

CrossRef Full Text | Google Scholar

Hora, S. C., and Iman, R. L. (1986). A Comparison of maximum/bounding and Bayesian/Monte Carlo for fault tree uncertainty analysis. Albuquerque, NM, USA: Sandia National Laboratories. Technical Report SAND85-2839.

Google Scholar

Iman, R. L., Helton, J. C., and Campbell, J. E. (1981a). An approach to sensitivity analysis of computer models: Part I—introduction, input variable selection and preliminary variable assessment. J. Qual. Technol. 13, 174–183. doi:10.1080/00224065.1981.11978748

CrossRef Full Text | Google Scholar

Iman, R. L., Helton, J. C., and Campbell, J. E. (1981b). An approach to sensitivity analysis of computer models: Part II—ranking of input variables, response surface validation, distribution effect and technique synopsis. J. Qual. Technol. 13, 232–240. doi:10.1080/00224065.1981.11978763

CrossRef Full Text | Google Scholar

Kramer, M. A., Calo, J. M., and Rabitz, H. (1981). An improved computational method for sensitivity analysis: green’s Function Method with “AIM”. Appl. Math. Model. 5, 432–441. doi:10.1016/s0307-904x(81)80027-3

CrossRef Full Text | Google Scholar

Lewins, J. (1965). IMPORTANCE: the adjoint function. Oxford, UK: Pergamon Press Ltd.

Google Scholar

Luo, Z., Wang, X., and Liu, D. (2020). Prediction on the static response of structures with large-scale uncertain-but-bounded parameters based on the adjoint sensitivity analysis. Struct. Multidiscip. Optim. 61, 123–139. doi:10.1007/s00158-019-02349-w

CrossRef Full Text | Google Scholar

Práger, T., and Kelemen, F. D. (2014). “Adjoint methods and their application in earth sciences,” in Advanced numerical methods for complex environmental models: needs and availability. Editors I. Faragó, Á. Havasi, and Z. Zlatev (Oak Park, IL, USA: Bentham Science Publishers), 203–275. Chapter 4A.

Google Scholar

Rios Insua, D. (1990). Sensitivity analysis in multiobjective decision making. New York, USA: Springer Verlag.

Google Scholar

Saltarelli, A., Chan, K., and Scott, E. M. (2000). Sensitivity analysis (Chichester, UK: J. Wiley and Sons Ltd).

Google Scholar

Shultis, J. K., and Faw, R. E. (2000). Radiation shielding. La Grange Park, Illinois, USA: American Nuclear Society.

Google Scholar

Stacey, W. M. (2001). Nuclear reactor physics. New York. USA: John Wiley and Sons.

Google Scholar

Tukey, J. W. (1957). The propagation of errors, fluctuations and tolerances. Princeton, NJ, USA: Princeton University. Technical Reports No. 10-12.

Google Scholar

Valentine, T. E. (2006). Polyethylene-reflected plutonium metal sphere subcritical noise measurements, SUB-PU-METMIXED-001. International handbook of evaluated criticality safety benchmark experiments, NEA/NSC/DOC(95)03/I-IX, organization for economic Co-operation and development (OECD). Paris, France: Nuclear Energy Agency.

Google Scholar

Weiberg, A. M., and Wigner, E. P. (1958). The physical theory of neutron chain reactors. Chicago, Illinois, USA: University of Chicago Press.

Google Scholar

Weisbin, C. R., Oblow, E. M., Marable, J. H., Peelle, R. W., and Lucius, J. L. (1978). Application of sensitivity and uncertainty methodology to fast reactor integral experiment analysis. Nucl. Sci. Eng. 66, 307–333. doi:10.13182/nse78-3

CrossRef Full Text | Google Scholar

Wigner, E. P. (1945). Effect of small perturbations on pile period. Chicago, IL, USA: Scientific Research Publishing. Chicago Report CP-G-3048.

Google Scholar

Williams, M. L. (1986). “Perturbation theory for nuclear reactor analysis,” in Handbook of nuclear reactor calculations. Editor Y. Ronen (Boca Raton, Florida, USA: CRC Press), 63–188. Volume 3.

Google Scholar

Williams, M. L., and Engle, W. W. (1977). The concept of spatial channel theory applied to reactor shielding analysis. Nucl. Sci. Eng. 62, 92–104. doi:10.13182/nse77-a26941

CrossRef Full Text | Google Scholar

Wilson, W. B., Perry, R. T., Shores, E. F., Charlton, W. S., Parish, T. A., Estes, G. P., et al. (2002). “SOURCES4C: a code for calculating (α,n), spontaneous fission, and delayed neutron sources and spectra,” in Proceedings of the American Nuclear Society/Radiation Protection and Shielding Division 12th Biennial Topical Meeting, Santa Fe, NM, USA, 14–18 April 2002.

Google Scholar