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

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


Introduction
The computational model of a physical system comprises the following conceptual components: (a) a well-posed system of that relate the system's independent variables and parameters to the system's state (i.e., dependent) variables; (b) probability distributions, moments thereof, inequality and/or equality constraints that define the range of variations of the system's parameters; and (c) one or several quantities, customarily referred to as system responses (or objective functions, or indices of performance), which are computed using the mathematical model. This works presents a new, general-purpose methodology for computing exactly and efficiently functional derivatives (called "sensitivities") of results ("system responses"), predicted by nonlinear mathematical models of systems (physical, engineering, biological) involving imprecisely known (i.e., uncertain) parameters, including input data, correlations, initial and/or boundary conditions, as well as manufacturing tolerances that affect domain of the model's definition in phase space. This new method is called the "fourth-order comprehensive sensitivity analysis methodology for nonlinear systems" (abbreviated as "4th-CASAM-N") since it enables the hitherto very difficult, if not intractable, exact computation of all of the first-, second-, third-, and fourth-order response sensitivities for large-scale systems involving many parameters, as is usually encountered in practice. The foundation for the material presented in this work is provided by the first-order adjoint sensitivity analysis procedure for nonlinear systems that was originally formulated in a general, functional analytic framework by Cacuci [1,2].
The aims and means of sensitivity theory/analysis are occasionally confused with the aims and means of optimization theory. The algorithms underlying optimization theory compute the values at which the first-order derivatives of a response with respect to the state functions and/or parameters vanish. In contradistinction, sensitivity theory/analysis aims at computing the response sensitivities to parameters at the nominal values for the model's parameters and state functions. Therefore, although both sensitivity analysis and optimization algorithms evaluate first-order (and occasionally higher-order) response derivatives with respect to model parameters, these algorithms serve different purposes and conceptually differ from each other. Response sensitivities to model parameters, which are computed using the methods of sensitivity analysis, are needed in many activities, including: (i) determining the effects of parameter variations on the system's behavior; (ii) ranking the importance of model parameters in influencing the system response under consideration; (iii) quantifying uncertainties induced in responses by parameter uncertainties (e.g., by using the method of "propagation of uncertainties"); (iv) prioritizing possible improvements for the system under consideration and possibly reducing conservatism and redundancy; (v) validating the model under consideration by comparison to experiments, while taking into account both experimental and model uncertainties; and (vi) performing "predictive modeling" (which includes data assimilation and model calibration) for the purpose of obtaining best-estimate predicted results with reduced predicted uncertainties, using, e.g., the methodologies presented in [3][4][5] and in references therein.
First-order response sensitivities can be computed by using either statistical or deterministic methods; a comparative review of the most popular of these methods was presented in [6,7]. It is known that sensitivities cannot be computed exactly using statistical methods; this can only be performed with deterministic methods. Furthermore, for a system comprising TP parameters, the computation of the first-order sensitivities by statistical methods requires O(TP) large-scale computations while the adjoint sensitivity analysis originally developed by Cacuci [1,2] requires a single large-scale computation per response.
Since nonlinear operators do not admit bona-fide adjoint operators (only a linearized form of a nonlinear operator admits an adjoint operator), responses of nonlinear models can depend only on the forward functions. In contradistinction, model responses for linear systems may involve the solutions of both the forward and the adjoint linear models that correspond to the respective physical system. Hence, responses for linear systems cannot always be treated as particular cases of nonlinear systems but there is a need to develop a dedicated sensitivity analysis methodology for response-coupled forward and adjoint linear systems. The general methodology for computing arbitrarily high-order sensitivities for response-coupled linear forward/adjoint systems was developed by Cacuci [8][9][10][11][12]. The overwhelming impact of the higher-order (i.e., second-, third-and fourth-order) sensitivities on the model response was illustrated by Fang and Cacuci, in [13] and references therein, by means of an OECD/NEA reactor physics benchmark modeled by the linear neutron transport equation and comprising 21,976 uncertain model parameters.
The general mathematical framework for adjoint sensitivity analysis of nonlinear systems was extended by Cacuci to second-order [14] and third-order [15] sensitivities. The extension of the adjoint sensitivity analysis methodology to include the computation of sensitivities of model responses with respect to imprecisely known domain boundaries was presented in [16]. By extending all of the previous theoretical developments of the adjoint sensitivity analysis methodology, this work presents the theoretical framework for the fourth-order comprehensive sensitivity analysis methodology for nonlinear systems (abbreviated as "4th-CASAM-N"). This work is structured as follows. Section 2 presents the mathematical formulation of the that would define the computational model of a generic nonlinear physical system, including the definition of a generic response which depends on the model's state variables and parameters, which are considered to be uncertain (i.e., not known precisely). Besides initial conditions and correlations, the model parameters are also considered to include geometrical parameters that describe the system's boundaries and internal interfaces. Section 3 presents the 4th-CASAM-N methodology while Section 4 summarizes the salient features of this novel methodology. The application of the principles underlying the 4th-CASAM-N is illustrated in an accompanying work [17] by means of a paradigm nonlinear heat conduction model.

Generic Mathematical Modeling of a Nonlinear Physical System
As already mentioned, the computational model of a physical system comprises that relate the system's independent variables and parameters to the system's state variables. The model parameters usually stem from processes that are external to the system under consideration and are seldom, if ever, known precisely. The known characteristics of the model parameters may include their nominal (expected/mean) values and, possibly, higher-order moments or cumulants (i.e., variance/covariances, skewness, kurtosis), which are usually determined from experimental data and/or processes external to the physical system under consideration. Occasionally, only inequality and/or equality constraints that delimit the ranges of the system's parameters are known. Without loss of generality, the imprecisely known model parameters can be considered to be real-valued scalar quantities. These model parameters will be denoted as α 1 , . . . ,α TP , where TP denotes the "total number of imprecisely known parameters" underlying the model under consideration. For subsequent developments, it is convenient to consider that these parameters are components of a "vector of parameters" denoted as α (α 1 , . . . , α TP ) † ∈ E α ∈ R TP , where E α is also a normed linear space and where R TP denotes the TP-dimensional subset of the set of real scalars. The components of the TP-dimensional column vector α ∈ R TP are considered to include imprecisely known geometrical parameters that characterize the physical system's boundaries in the phase space of the model's independent variables. Matrices will be denoted using capital bold letters while vectors will be denoted using either capital or lower-case bold letters. The symbol " " will be used to denote "is defined as" or "is by definition equal to". Transposition will be indicated by a dagger ( †) superscript. The generic nonlinear model is considered to comprise TI independent variables which will be denoted as x i , i = 1, . . . , TI, and which are considered to be components of a TI-dimensional column vector denoted as x (x 1 , . . . , x TI ) † ∈ R TI , where the sub/superscript "TI" denotes the "total number of independent variables". The vector x ∈ R TI of independent variables is considered to be defined on a phase-space domain which will be denoted as Ω(α) and which is defined as follows: . . , TI}. The lower boundary point of an independent variable is denoted as λ i (α) and the corresponding upper boundary point is denoted as ω i (α). A typical example of boundaries that depend on both geometrical parameters and material properties are the "boundaries facing vacuum" in models based on diffusion theory, where conditions are imposed on the "extrapolated boundary" of the respective spatial domain. The "extrapolated boundary" depends both on the imprecisely known physical dimensions of the problem's domain and also on the medium's properties, such as atomic number densities and microscopic transport cross sections. The boundary of Ω(α), which will be denoted as ∂Ω(α), comprises the set of all of the endpoints λ i (α), ω i (α), i = 1, . . . , TI, of the respective intervals on which the components of x are defined, i.e., A nonlinear physical system can be generally modeled by means of coupled which can be represented in operator form as follows: The quantities which appear in Equation (1) are defined as follows: . . , u TD (x)] † is a TD-dimensional column vector of dependent variables; the abbreviation "TD" denotes "total number of dependent variables". The functions u i (x), i = 1, . . . , TD denote the system's "dependent variables" (also called "state functions"); u(x) ∈ E u , where E u is a normed linear space over the scalar field F of real numbers.

N[u(x)
; α] [N 1 (u; α), . . . , N TD (u; α)] † denotes a TD-dimensional column vector. The components N i (u; α), i = 1, . . . , TD are operators (including differential, difference, integral, distributions, and/or finite or infinite matrices) acting (usually) nonlinearly on the dependent variables u(x), the independent variables x, and the model parameters α. The mapping N(u; α) is defined on the combined domains of the model's parameters and state functions, i.e., N : . . ., q TD (x; α)] † is a TD-dimensional column vector which represents inhomogeneous source terms, which usually depend nonlinearly on the uncertain parameters α. The vector Q(x, α) is defined on a normed linear space denoted as E Q , i.e., The equalities in this work are considered to hold in the weak ("distributional") sense. The right sides of Equation (1) and of other various to be derived in this work may contain "generalized functions/functionals", particularly Dirac distributions and derivatives thereof.
Boundary and/or initial conditions must also be provided if differential operators appear in Equation (1). In operator form, these boundaries and/or initial conditions are represented as follows: where the column vector 0 has TD components, all of which are identically zero, i.e., In Equation (2), . . , B TD (u; α)] † are nonlinear operators in u(x) and α, which are defined on the boundary ∂Ω x (α) of the model's domain Ω x (α). The components C i (x; α), i = 1, . . . , TD of C(x, α) [C 1 (x; α), . . . , C TD (x; α)] † comprise inhomogeneous boundary sources which are nonlinear functions of α.
Solving Equations (1) and (2) at the nominal parameter values, denoted as α 0 α 0 1 , . . . , α 0 i , .., α 0 TP † , provides the "nominal solution" u 0 (x), i.e., the vectors u 0 (x) and α 0 satisfy the following: The superscript "0" will be used throughout this work to denote "nominal values". The results computed using a mathematical model are customarily called "model responses" (or "system responses" or "objective functions" or "indices of performance"). In general, a function-valued (i.e., operator-type) response R[u(x); α] can be represented by a spectral expansion in multidimensional orthogonal polynomials or Fourier series of the form: where the quantities P m i (x i ), i = 1, . . . , TI denote the corresponding spectral functions (e.g., orthogonal polynomials or Fourier exponential/trigonometric functions) and where the spectral Fourier coefficients c m 1 ...m TI [u(x); α] are defined as follows: . . .
The coefficients c m 1 ...m TI [u(x); α] can themselves be considered as "model responses" since the spectral polynomials P m i (x i ) are perfectly well known while the expansion coefficients will contain all of the dependencies (directly or indirectly-through the state functions) of the respective response on the imprecisely known model parameters. This way, the sensitivity analysis of an operator-valued response R[u(x); α] can be reduced to the sensitivity analysis of the scalar-valued responses c m 1 ...m TI [u(x); α].
A measurement of a physical quantity that depends on the model's state functions and parameters can be considered to be a response denoted as R p [u(x); α], which is to denotes the location in the phase space of the specific "measurement point". Such a measurement (or measurement-like) response can be represented mathematically as follows: where the function F[u(x); α; x] denotes the mathematical dependence of the measurement device on the model's dependent variable(s), and where the quantity δ x i − x p i (α) denotes the Dirac delta functional. The measurement's location in phase space, x p (α), may itself be afflicted by measurement (experimental) uncertainties. Hence, it is convenient to consider the components of x p (α) to be included among the components of the vector α of model parameters, even though x p (α) appears only in the definition of the response but does not appear in Equations (1) and (2), which mathematically define the physical model. Thus, the physical "system" is defined to comprise both the system's computational model and the system's response. In most cases, the coordinates x p k (α), k = 1, . . . , TI will be independent (albeit uncertain) model parameters, in which case ∂x The representations shown in Equations (6)- (8) indicate that model responses can be fundamentally analyzed by considering the following generic integral representation: where S[u(x); α] is suitably differentiable nonlinear function of u(x) and of α. It is important to note that the components of α not only include parameters that appear in the defining the computational model per se, i.e., in Equations (1) and (2), but also include parameters that specifically occur only in the definition of the response under consideration. It is also important to note that the system's definition domain, Ω(α), in phase space is considered to be imprecisely known, subject to uncertainties in the components of the vector of model parameters α. Therefore, the system domain's boundary, ∂Ω(α), as well as the model response R[u(x); α], will be affected by the boundary uncertainties that affect the endpoints λ i (α), ω i (α), i = 1, . . . , TI. Such boundary uncertainties stem most often from manufacturing uncertainties.

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

The First-Order Comprehensive Adjoint Sensitivity Analysis Methodology for Nonlinear Systems (1st-CASAM-N)
The model and boundary parameters α are considered to be uncertain quantities, having unknown true values. The nominal (or mean) parameter vales α 0 are considered to be known, and these will differ from the true values by quantities denoted as δα (δα 1 , . . . , δα TP ), where δα i α i − α 0 i . Since the forward state functions u(x) are related to the model and boundary parameters α through Equations (1) and (2), it follows that the variations δα in the model and boundary parameters will cause corresponding variations v (1) (x) [δu 1 (x), . . . , δu TD (x)] † around the nominal solution u 0 (x) in the forward state functions. In turn, the variations δα and v (1) (x) will induce variations in the system's response. Cacuci [1,2] has shown that the most general definition of the sensitivity of an operator-valued model response R(e), where e (α, u) ∈ E, to variations h δα, v (1) in the model parameters and state functions in a neighborhood around the nominal functions and parameter values e 0 α 0 , u 0 ∈ E, is given by the first-order Gateaux (G) variation, which will be denoted as δR e 0 ; h and is defined as follows: for a scalar ε ∈ F and for all (i.e., arbitrary) vectors h ∈ E = E α × E u in a neighborhood e 0 + εh around e 0 = α 0 , u 0 ∈ E. The G variation δR e 0 ; h is an operator defined on the same domain as R(e) and has the same range as R(e). The G variation δR e 0 ; h satisfies the relation: The existence of the G variation δR e 0 ; h does not guarantee its numerical computability.
Numerical methods most often require that δR e 0 ; h is linear in h δα; v (1) in a neighborhood e 0 + εh around e 0 = u 0 , α 0 ∈ E. Formally, the necessary and sufficient conditions for the G variation δR e 0 ; h of a nonlinear operator R(e) to be linear and continuous in h in a neighborhood e 0 + εh around e 0 = α 0 , u 0 and therefore admit a total first-order G derivative, are as follows: (i) R(e) satisfies a weak Lipschitz condition at e 0 , namely: (ii) R(e) satisfies the following condition In practice, the relation provided in Equations (11) and (12) are seldom used directly since the computation of the expression on the right side of Equation (10) reveals if the respective expression is linear (or not) in h δα, v (1) and, hence, in v (1) (x). Numerical methods (e.g., Newton's method, and variants thereof) for solving Equations (1) and (2) also require the existence of the first-order G derivatives of original model, in which case the components of the operators which appear in these must also satisfy the conditions provided in Equations (11) and (12). Therefore, the conditions provided in Equations (11) and (12) are henceforth considered to be satisfied by the operators which underly the physical system modeled by Equations (1) and (2), as well as by the model responses.
When the first-order G differential δR e 0 ; h satisfies the conditions provided in Equations (11) and (12), it can be written in the following form: where This, the quantities ∂R(u; α)/∂u and ∂R(u; α)/∂α in Equation (13) denote the partial G derivatives of R(e) with respect to u and α, evaluated at the nominal parameter values (and hence also nominal values of the state functions). The notation { } α 0 will be used in this work to indicate that the quantity enclosed within the bracket is to be evaluated at the respective nominal parameter and state functions values. The quantity {∂R(u; α)/∂α} α 0 δα is called the "direct-effect term" because it arises directly from parameter variations δα. The direct-effect term can be computed once the nominal values e 0 = u 0 , α 0 are available. The quantity {∂R(u; α)/∂u} α 0 v (1) (x) is called the "indirect-effect term" because it arises indirectly, through the variations in the state functions (which are caused through the model by parameter variations). The indirect-effect term can be quantified only after having determined the variations v (1) (x) in terms of the variations δα.
The first-order relationship between the vectors v (1) (x) and δα is determined by taking the G differentials of Equations (1) and (2). Thus, applying the definition of the G differential to Equations (1) and (2) yields the following: Carrying out the differentiations with respect to ε in Equations (20) and (21), and setting ε = 0 in the resulting expressions yields the following: In Equations (22) and (23), the superscript "(1)" indicates "1st-Level" and the various quantities which appear in these are defined as follows: The system comprising Equations (22) and (23) is called the "1st-Level Variational Sensitivity System" (1st-LVSS). In order to determine the solutions of the 1st-LVSS that would correspond to every parameter variation δα j 1 , j 1 = 1, . . . , TP, the 1st-LVSS would need to be solved TP times, with distinct right sides for each δα j 1 , as follows: Subsequently, the solutions v (1) (j 1 ; x) could be used, in turn, in Equation (17) to compute the indirect-effect term corresponding to each parameter variation δα j 1 , j 1 = 1, . . . , TP, to obtain the following contribution from the indirect-effect term to the respective partial sensitivity of the response with respect to the parameter α j 1 : Adding the contribution from the indirect-effect term obtained in Equation (31) to the contribution from the direct-effect term provided in Equation (18) yields the following expression for the sensitivity (i.e., partial G derivative) ∂R[u(x); α]/∂α j 1 of the response R[u(x); α] to the parameter α j 1 , j 1 = 1, . . . , TP: The quantities ∂R[u(x); α]/∂α j 1 are independent of the parameter variations δα i 1 and represent the first-order partial sensitivities (first-order partial G derivatives) of the response R(e) with respect to each of the model parameters α j 1 , j 1 = 1, . . . , TP, evaluated at the nominal values e 0 = α 0 , u 0 . Computing the response sensitivities by using the solutions v (1) (j 1 ; x), j 1 = 1, . . . , TP, of the 1st-LVSS requires TP large-scale forward computations in order to determine the functions v (1) (j 1 ; x), j 1 = 1, . . . , TP. Since most problems of practical interest are characterized by many parameters (i.e., α has many components) and comparatively few responses, it becomes prohibitively expensive to solve repeatedly the 1st-LVSS in order to determine the functions v (1) (j 1 ; x), j 1 = 1, . . . , TP. Even though the 1st-LVSS contains first-order parameter and state-functions variations, it is called "firstlevel" (rather than "first-order") in anticipation of determining second-order sensitivities, which will use "second-level" forward and adjoint systems. These "second-level" systems will not be called "second-order" because they will not contain second-order parameter and/or state-function variations, but will also contain only first-order variations, even though they will be used for determining second-order sensitivities. Similar terminology, i.e., "third-level" (as opposed to "third-order") forward/adjoint systems, will be used for determining the third-order sensitivities, and so on.
In most practical situations, the number of model parameters significantly exceeds the number of functional responses of interest, i.e., TR TP, so it would be advantageous to perform just TR (rather than TP) computations. The goal of the "1st-order comprehensive adjoint sensitivity analysis methodology for nonlinear systems (1st-CASAM-N)" is to compute exactly and efficiently the "indirect effect term" defined in Equation (17) without needing to compute explicitly the vectors v (1) (j 1 ; x), j 1 = 1, . . . , TP. The qualifier "comprehensive" is meant to indicate that that the 1st-CASAM-N considers that the internal and external boundaries ∂Ω x (α) of the phase-space domain depend on the uncertain model parameters α and are thereby imprecisely known, subject to uncertainties. Thus, the 1st-CASAM-N represents a generalization of the pioneering works by Cacuci [1,2] that conceived the "adjoint sensitivity analysis methodology", in which the domain boundary was considered to be perfectly well known, free of uncertainties. The fundamental ideas underlying the 1st-CASAM-N are as conceived by Cacuci [1,2], aiming at eliminating the appearance of the vectors v (1) (j 1 ; x) from the expression of the indirect-effect term defined in Equation (17). This elimination is achieved by expressing the right side of Equation (17) in terms of the solutions of the "1st-Level Adjoint Sensitivity System (1st-LASS)", the construction of which requires the introduction of adjoint operators. Adjoint operators can be defined in Banach spaces but are most useful in Hilbert spaces. Since real Hilbert spaces provide the natural mathematical setting for computational purposes, the derivations presented in this section are set in real (as opposed to complex) Hilbert spaces, without affecting the generality of the concepts presented herein. Thus, the spaces E u and E Q are henceforth considered to be self-dual Hilbert spaces and will be denoted as H 1 (Ω x ). The inner product of two vectors u (a) (x) ∈ H 1 and u (b) (x) ∈ H 1 will be denoted as u (a) , u (b) 1 , and is defined as follows: where the dot indicates the "scalar product of two vectors" defined as follows: i (x). It is important to note that the inner product defined in Equation (33) is continuous in α, i.e., it holds at any value particular value of α, including at the nominal parameter values α 0 .
The construction of the 1st-LASS commences by noting that the vector v (1) The next step is to form the inner product of Equation (22) with a vector TD (x) ∈ H 1 , where the superscript "(1)" indicates "1st-Level": Using the definition of the adjoint operator in H 1 (Ω x ), the left side of Equation (35) is transformed as follows: where P (1) u; α; a (1) ; . The symbol [ ] * will be used in this work to indicate "adjoint" operator. In certain situations, it might be computationally advantageous to include certain boundary components of P (1) u; α; a (1) ; The domain of A (1) (u; α) is determined by selecting appropriate adjoint boundary and/or initial conditions, which will be denoted in operator form as: The above boundary conditions for A (1) (u; α) are usually inhomogeneous, i.e., b (1) A (0; 0; α) = 0, and are obtained by imposing the following requirements: 1. They must be independent of unknown values of v (1) (x) and δα; 2. The substitution of the boundary and/or initial conditions represented by Equations (23) and (37) into the expression of Constructing the adjoint initial and/or boundary conditions for A (1) (u; α) as described above and implementing them together with the variational boundary and/or initial conditions represented by Equation (23) into Equation (35) reduces the bilinear con- to a quantity that will contain boundary terms involving only known values of δα, α 0 , u 0 , and ψ (1) ; this quantity will be denoted by In general, the boundary terms represented by do not vanish automatically. In certain cases, however, may vanish automatically or it may be forced to vanish by considering appropriately constructed extensions of the adjoint operator A (1) (α, u); however, such extensions are seldom needed in practice. Since P (1) u; α; a (1) ; δα can be expressed in the following form: Implementing the forward and adjoint boundary and/or initial conditions given in Equations (23) and (37) into Equation (36) transforms the later into the following relation: Replacing the quantity V (1) (α; u)v (1) in the first term on the right side of Equation (38) by the right side of Equation (22) yields the following relation: The definition of the function a (1) (x) will now be completed by requiring that the left side of Equation (39) is the same as the indirect-effect term defined in Equation (17), which is achieved by imposing the following relationship: while satisfying the adjoint boundary conditions represented by Equation (37). The subscript "A" attached to the source term on the right side of Equation (40) indicates "adjoint".
Since the source s (1) A [u(x); α] may contain distributions (e.g., Dirac delta functions and derivatives thereof), the equality in Equation (40) is considered to hold in the weak sense. The well-known Riesz representation theorem ensures that the relationship in Equation (40) holds uniquely.
The results obtained in Equations (38)-(40) are now replaced in Equation (17) to obtain the following expression of the indirect-effect term as a function of a (1) where, for each j 1 = 1, . . . , TP, the contribution of the indirect-effect term to the sensitivity of the response with respect to the parameter α j 1 is given by As the identity on the right side of Equation (41) indicates, the desired elimination of all unknown values of v (1) ind now depends on the adjoint function a (1) (x) ∈ H 1 . Replacing in Equation (16) the result obtained in Equation (41) together with the expression provided (18) for the direct-effect term yields the following expression δR u(x); α; v (1) The expressions of the first-order response sensitivities ∂R[u(x); α]/∂α j 1 α 0 of the response R[u(x); α] with respect to the parameters α j 1 are obtained identifying in Equation (43) the quantities that multiply the respective parameter variations δα j 1 , j 1 = 1, . . . , TP. This identification yields the following expressions for the first-order response sensitivities ∂R[u(x); α]/∂α j 1 α 0 , j 1 = 1, . . . , TP, computed at the model's nominal parameter and state function values: . . .
As indicated by Equation (44), each of the first-order sensitivities R (1) j 1 ; u(x); a (1) with respect to the model parameters α j 1 (including boundary and initial conditions) can be computed inexpensively after having obtained the function a (1) (x) ∈ H 1 , using just quadrature formulas to evaluate the various inner products involving a (1) (x) ∈ H 1 in the expression of the indirect-effect term obtained in Equation (42). The function a (1) (x) ∈ H 1 is obtained by solving numerically Equations (37) and (40), which is the only large-scale computation needed for obtaining all of the first-order sensitivities. Equations (37) and (40) will be called the first-level adjoint sensitivity system (1st-LASS), and its solution, a (1) (x) ∈ H 1 (Ω x ), will be called the first-level adjoint function. It is very important to note that the 1st-LASS is independent of parameter variation δα j 1 , j 1 = 1, . . . , TP, and therefore needs to be solved only once, regardless of the number of model parameters under consideration. Furthermore, since Equation (40) is linear in a (1) (x)ψ (2) 1,i 1 (x), solving it requires less computational effort than solving the original Equation (1), which is nonlinear in u(x).

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

B
(2) The argument "2" which appears in the list of arguments of the vector U (2) (2; x) and the "variational vector" V (2) (2; x) in Equations (57) and (58) indicates that each of these vectors is a two-block column vector, with each block comprising a column vector of dimension TD, defined as follows: Thus, the (column) block vector U (2) (2; x) has a total of 2 × TD components; evidently, the (column) block vector V (2) (2; x) also has a total of 2 × TD components. In the relatively simple case regarding the components of either the vector U (2) (2; x) or the vector V (2) (2; x), the numbers "1" and "2" could also be used as subscripts, but such a subscript notation would become unwieldy for the higher-level (adjoint) functions, which will be introduced in the sections to follow below. The superscript "(2)" which appears in the notation of the vectors U (2) (2; x) and V (2) (2; x) indicates "2nd-level". Henceforth, such "higher-level" (i.e., level higher than first) variational and adjoint functions/vectors will be denoted using bold capital letters. The argument "2" in the expression 0 [2] [0, 0] † indicates that the quantity 0[2] is a two-block column vector comprising two vectors, each of which has TD components, all of which are zero-valued, as defined in Equation (3). Thus, the column vector 0[2] has a total of 2 × TD components, all of which are identically zero.
To distinguish block vectors from block matrices, two capital bold letter are used (and will henceforth be used) to denote block matrices, as in the case of "the second-level" "variational matrix" VM (2) 2 × 2; u (2) (x); α . The "2nd-level" is indicated by the superscript "(2)". Subsequently in this work, levels higher than second will also be indicated by a corresponding superscript attached to the appropriate block vectors and/or block matrices.
The argument "2 × 2", which appears in the list of arguments of VM (2) 2 × 2; u (2) (x); α , indicates that this matrix is a 2 × 2-dimensional block matrix comprising four matrices, each with of dimensions TD × TD, having the following structure: Thus, the matrix VM (2) 2 × 2; u (2) (x); α has a total of (2 × TD) 2 components (or elements). The other quantities which appear in Equations (57) and (58) are also two-block vectors, with the same structure as V (2) (2; x), and are defined as follows: Solving the 2nd-LVSS requires TP 2 large-scale computations, which is unrealistic to perform for large-scale systems comprising many parameters. The 2nd-CASAM-N circumvents the need for solving the 2nd-LVSS by deriving an alternative expression for the indirect-effect term defined in Equation (47), in which the function V (2) (2; x) is replaced by a second-level adjoint function which is independent of variations in the model parameter and state functions. This second-level adjoint function will satisfy a second-level adjoint sensitivity system (2nd-LASS), which will be constructed by using the 2nd-LVSS as the starting point, following the same principles outlined in Section 3.1 The 2nd-LASS will be constructed in a Hilbert space which will be denoted as H 2 (Ω x ) and which comprises as elements block vectors of the same form as V (2) (2; x). Thus, a generic vector in H 2 (Ω x ), denoted as Ψ (2) comprises two components of the form ψ (2) (1; x) ψ 1 (1; x), . . . , ψ 1 (2; x), . . . , ψ The inner product of two vectors Ψ (2) and defined as follows: The inner product defined in Equation (63) is continuous in α in a neighborhood of α 0 . Using the definition of the inner product defined in Equation (63), construct the inner product of Equation (57) with a vector A (2) (2; x) a (2) (1; x), a (2) (2; x) † ∈ H 2 (Ω x ) to obtain the following relation: The inner product on the left side of Equation (64) is now further transformed by using the definition of the adjoint operator to obtain the following relation: (2) ; A (2) ; V (2) ; α where the adjoint matrix-valued operator AM (2) 2 × 2; u (2) (x); α is defined as follows: The matrix AM (2) 2 × 2; u (2) (x); α comprises (2 × 2) block matrices, each with dimensions TD 2 , thus comprising a total of (2 × 2)TD 2 components (or elements).
The system of equations represented by Equations (70)-(72) will be called the secondlevel adjoint sensitivity system (2nd-LASS) and its solution, , will be called the second-level adjoint function. The 2nd-LASS is independent of parameter variations δα and variations V (2) (2; x) in the respective state functions. It is also important to note that the (2 × TD) 2 -dimensional matrix AM (2) 2 × 2; U (2) (2; x); α is independent of the index j 1 . Only the source term Q (2) A 2; j 1 ; U (2) (2; x); α depends on the index j 1 . Therefore, the same solver can be used to invert the AM (2) 2 × 2; U (2) (2; x); α and numerically solve the 2nd-LASS for each j 1dependent source Q A 2; j 1 ; U (2) (2; x); α , for each index j 1 , in order to obtain the corresponding j 1 -dependent 2 × TD-dimensional second-level adjoint function A (2) (2; j 1 ; x) a (2) (1; j 1 ; x), a (2) The two components a (2) (1; j 1 ; x) and a (2) (2; j 1 ; x) of the second-level adjoint function are distinguished from each other by the use of the numbers "1" and, respectively, "2" in the respective list of arguments. In this particularly simple case, the numbers "1" and "2" could also be used as subscripts, in the customary notation for vector components, but such a use would not lend itself to generalizations because the subscript notation would become unwieldy for the higher-level adjoint functions, which will be introduced in the sections that follow below.
Since the adjoint matrix AM (2) 2 × 2; U (2) (2; x); α is block diagonal, solving the 2nd-LASS is equivalent to solving two 1st-LASS, with two different source terms. Thus, the "solvers" and the computer program used for solving the 1st-LASS can also be used for solving the 2nd-LASS. The 2nd-LASS was designated as the "second-level" rather than the "second-order" adjoint sensitivity system, since the 2nd-LASS does not involve any explicit second-order G derivatives of the operators underlying the original system but involves the inversion of the same operators that needed to be inverted for solving the 1st-LASS.
It is important to note that if the 2nd-LASS is solved TP-times, the second-order mixed sensitivities R (2) j 2 ; j 1 ; U (2) (2; x); A (2) (2; j 1 ; x); α ≡ ∂ 2 R/∂α j 2 ∂α j 1 will be computed twice, in two different ways, in terms of two distinct second-level adjoint functions. Consequently, the symmetry property ∂ 2 R[u(x); α]/∂α j 2 ∂α j 1 =∂ 2 R[u(x); α]/∂α j 1 ∂α j 2 enjoyed by the second-order sensitivities provides an intrinsic (numerical) verification that the components of the second-level adjoint function A (2) (2; j 1 ; x), as well as the first-level adjoint function a (1) (x) are computed accurately. The structure of the 2nd-LASS enables full flexibility for prioritizing the computation of the second-order sensitivities. The computation of the second-order sensitivities would logically be prioritized based on the relative magnitudes of the first-order sensitivities: the largest relative first-order response sensitivity should have the highest priority for computing the corresponding second-order mixed sensitivities; then, the second largest relative first-order response sensitivity should be considered next, and so on. The unimportant second-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them. Computing second-order sensitivities that correspond to vanishing first-order sensitivities may also be of interest, since vanishing first-order sensitivities may indicate critical points of the response in the phase space of model parameters.
Concatenating Equations (82) and (83) with the 2nd-LVSS represented by Equations (57) and (58) yields the following system of, which will be called the "3rdorder variational sensitivity system" (3rd-LVSS), for determining the vectors V (2) (2; x) and δA (2)  (91) For subsequent reference, it is important to note that the quantities Solving the 3rd-LVSS would require TP 3 large-scale computations, which is unrealistic for large-scale systems comprising many parameters. The 3rd-CASAM-N circumvents the need for solving the 3rd-LVSS by deriving an alternative expression for the indirect-effect term defined in Equation (81), in which the function V (3) (4; j 1 ; x) is replaced by a third-level adjoint function which is independent of parameter variations. This third-level adjoint function will be the solution of a third-level adjoint sensitivity system (3rd-LASS) which will be constructed by applying the same principles as those used for constructing the 1st-LASS and the 2nd-LASS. The Hilbert space appropriate for constructing the 3rd-LASS will be denoted as H 3 (Ω x ) and comprise element block vectors of the same form as V (3) (4; j 1 ; x). Thus, a generic block vector in H 3 (Ω x ), denoted as Ψ (3) where each of these four components is a TD-dimensional column vector. The inner product of two vectors Ψ (3) (4; x) ∈ H 3 (Ω x ) and Φ (3) (4; x) ∈ H 3 (Ω x ) in the Hilbert space H 3 (Ω x ) will be denoted as Ψ (3) (4; x), Φ (3) (4; x) 3 and defined as follows: The inner product defined in Equation (94) The inner product on the left side of Equation (95) is further transformed by using the definition of the adjoint operator to obtain the following relation: where and where P  (1; x), a (3) (2; x), a (3) (3; x), a (3) (4; x) † ∈ H 3 (Ω x ) satisfies adjoint boundary/initial conditions denoted as follows: The third-level adjoint boundary/initial conditions represented by Equation (98) (85) and (98) into Equation (96) will transform the latter into the following form: The definition of the third-level adjoint function A (3) (4; x) ∈ H 3 (Ω x ) is now completed by requiring that the left side of Equation (100) and the right side of Equation (81) represent the "indirect-effect term" δR (2) j 2 ; j 1 ; U (2) (2; x); A (2) (2; j 1 ; x); α; V (2) (2; x); δA (2) (2; j 1 ; x) ind for each of the indices j 1 = 1, . . . , TP; j 2 = 1, . . . ., j 1 .
The boundary conditions to be satisfied by each of the third-level adjoint functions  (101) is considered to hold in the weak sense. The Riesz representation theorem ensures that the weak equality in Equation (101) holds uniquely.
The matrix AM (3) 4 × 4; U (3) (4; j 1 ; x); α is block diagonal; therefore, solving the 3rd-LASS is equivalent to solving three 1st-LASS, with different source terms. The 3rd-LASS was designated as "the third-level" rather than "third-order" adjoint sensitivity system since the 3rd-LASS does not involve any explicit second-order and/or third-order G derivatives of the operators underlying the original system, but involves only the inversion of the same operators that needed to be inverted for solving the 1st-LASS.
By solving the 3rd-LASS TP(TP + 1)/2 times, the third-order mixed sensitivities ∂ 3 R(α, u)/∂α j 3 ∂α j 2 ∂α j 1 (α 0 ) will be computed three times, in three different ways. Consequently, the multiple symmetries intrinsic to the third-order sensitivities provide an intrinsic numerical verification that the components of the first-, second-, and third-level adjoint functions are computed accurately.
The structure of the 3rd-LASS enables full flexibility for prioritizing the computation of the third-order sensitivities. The computation of the third-order sensitivities would logically be prioritized based on the relative magnitudes of the second-order sensitivities, so that the unimportant third-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them.
The adjoint matrix AM (4) 8 × 8; U (4) (8; j 2 ; j 1 ; x); α is block diagonal; therefore, solving the 4th-LASS is equivalent to solving four 1st-LASS, with different source terms. The 4th-LASS was designated as the "fourth-level" rather than "fourth-order" adjoint sensitivity system since the 3rd-LASS does not involve any explicit second-order, third-order, and/or fourth-order G derivatives of the operators underlying the original system but involves the inversion of the operators similar to those that needed to be inverted for solving the 1st-LASS.
By solving the 4th-LASS TP(TP + 1)(TP + 2)/6 times, the fourth-order mixed sensitivities ∂ 4 R[u(x); α]/∂α j 1 ∂α j 3 ∂α j 3 ∂α j 4 will be computed four times, in four different ways using distinct adjoint functions. Consequently, the multiple symmetries intrinsic to the fourth-order sensitivities provide an intrinsic numerical verification that the components of the first-, second-, third-, and fourth-level adjoint functions are computed accurately.
The structure of the 4th-LASS enables full flexibility for prioritizing the computation of the third-order sensitivities. The computation of the fourth-order sensitivities would be prioritized based on the relative magnitudes of the third-order sensitivities, so that the unimportant fourth-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them.

Discussion
This work has presented the "fourth-order comprehensive sensitivity analysis methodology for nonlinear systems (abbreviated as "4th-CASAM-N"), which enables the hitherto very difficult, if not intractable, exact computation of all of the first-, second-, third-, and fourth-order response sensitivities for large-scale nonlinear systems involving many parameters. The qualifier "comprehensive" indicates that the fourth-CASAM-N methodology enables the exact and efficient computation not only of response sensitivities with respect to the customary model parameters (including computational input data, correlations, initial and/or boundary conditions) but also with respect to imprecisely known material boundaries, which would be caused by manufacturing tolerances.
It has been shown that the first-order sensitivities of the system response under consideration with respect to the model parameters (including boundary and initial conditions) can be computed inexpensively, using just quadrature formulas, after having obtained the first-level adjoint state function. The first-level adjoint state function is obtained by solving once the first-level adjoint sensitivity system (1st-LASS), which is the sole large-scale computation needed for obtaining all of the first-order sensitivities. This is because the 1st-LASS is independent of parameter variations, and therefore needs to be solved only once, regardless of the number of model parameters (denoted as "TP") under consideration. Furthermore, solving the 1st-LASS requires less computational effort than solving the nonlinear underlying the original models, since the 1st-LASS is linear in the first-level adjoint function.
The second-order sensitivities which correspond to each first-order sensitivity are also computed by using inexpensive quadrature formulas, after having obtained the secondlevel adjoint state function. The computation of the second-level adjoint state function requires solving the second-level adjoint sensitivity system (2nd-LASS). Since the 2nd-LASS is block diagonal, solving it is equivalent to solving two 1st-LASS, with two different source terms. Thus, the "solvers" and the computer program used for solving the 1st-LASS can also be used for solving the 2nd-LASS. The 2nd-LASS was designated as the "second-level" rather than the "second-order" adjoint sensitivity system, since the 2nd-LASS does not involve any second-order G derivatives of the operators underlying the original system but involves the inversion of the same operators that need to be inverted to solve the 1st-LASS. The structure of the 2nd-LASS enables full flexibility for prioritizing the computation of the second-order sensitivities. The computation of the second-order sensitivities would logically be prioritized based on the relative magnitudes of the firstorder sensitivities: the largest relative first-order response sensitivity should have the highest priority for computing the corresponding second-order mixed sensitivities; then, the second largest relative first-order response sensitivity should be considered next, and so on. The unimportant second-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them. Computing second-order sensitivities that correspond to vanishing first-order sensitivities may also be of interest, since vanishing first-order sensitivities may indicate critical points of the response in the phase space of model parameters. If the 2nd-LASS is solved TP times, the second-order mixed sensitivities will be computed twice, in two different ways, in terms of two distinct second-level adjoint functions. Consequently, the symmetry property enjoyed by the second-order sensitivities provides an intrinsic (numerical) verification that the components of the first-and secondlevel adjoint functions are computed accurately.
The exact computation of all of the partial third-order sensitivities that correspond to a second-order sensitivities is also accomplished by using quadrature formulas after having obtained the third-level adjoint function. Thus, for each second-order sensitivity, the third-level adjoint function is obtained by solving the third-level adjoint sensitivity system (3rd-LASS) once, which is equivalent to solving three 1st-LASS, with different source terms. The 3rd-LASS is designated as "the third-level" rather than the "third-order" adjoint sensitivity system since the 3rd-LASS does not involve any explicit second-order and/or third-order G derivatives of the operators underlying the original system. By solving the 3rd-LASS TP(TP + 1)/2 times, the third-order mixed sensitivities will be computed three times, in three different ways. Consequently, the multiple symmetries intrinsic to the third-order sensitivities provide an intrinsic numerical verification that the components of the first-, second-and third-level adjoint functions are computed accurately. The structure of the 3rd-LASS enables full flexibility for prioritizing the computation of the third-order sensitivities. The computation of the third-order sensitivities would logically be prioritized based on the relative magnitudes of the second-order sensitivities, so that the unimportant third-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them.
The exact computation of all of the partial fourth-order sensitivities that correspond to a third-order sensitivity also involves the use of quadrature formulas after having obtained the fourth-level adjoint function by solving the fourth-level adjoint sensitivity system (4th-LASS). The 4th-LASS is designated as the "fourth-level" rather than "fourth-order" adjoint sensitivity system since the 3rd-LASS does not involve any explicit second-order, thirdorder, and/or fourth-order G derivatives of the operators underlying the original system but involves the inversion of the operators similar to those that needed to be inverted for solving the 1st-LASS. The 4th-LASS is block diagonal; solving it is equivalent to solving four 1st-LASS, with different source terms. If the 4th-LASS is solved TP(TP + 1)(TP + 2)/6 times, the fourth-order mixed-response sensitivities will be computed four times, in four different ways using distinct adjoint functions. Consequently, the multiple symmetries intrinsic to the fourth-order sensitivities provide an intrinsic numerical verification that the components of the first-, second-, third-, and fourth-level adjoint functions are computed accurately. The structure of the 4th-LASS enables full flexibility for prioritizing the computation of the third-order sensitivities. The computation of the fourth-order sensitivities would be prioritized based on the relative magnitudes of the third-order sensitivities, so that the unimportant fourth-order sensitivities can be deliberately neglected while knowing the error incurred by neglecting them.
In summary, the implementation of the 4th-CASAM-N requires very little additional effort beyond the construction of the adjoint sensitivity system needed for computing the first-order sensitivities. An illustrative application of the 4th-CASAM-N to a paradigm nonlinear heat conduction is presented in the accompanying work [17].