Active subspace analysis and uncertainty quantification for a polydomain ferroelectric phase-field model

Abstract

We perform parameter subset selection and uncertainty analysis for phase-field models that are applied to the ferroelectric material lead titanate. A motivating objective is to determine which parameters are influential in the sense that their uncertainties directly affect the uncertainty in the model response, and fix noninfluential parameters at nominal values for subsequent uncertainty propagation. We employ Bayesian inference to quantify the uncertainties of gradient exchange parameters governing 180° and 90° tetragonal phase domain wall energies. The uncertainties of influential parameters determined by parameter subset selection are then propagated through the models to obtain credible intervals when estimating energy densities quantifying polarization and strain across domain walls. The results illustrate various properties of Landau and electromechanical coupling parameters and their influence on domain wall interactions. We employ energy statistics, which quantify distances between statistical observations, to compare credible intervals constructed using a complete set of parameters against an influential subset of parameters. These intervals are obtained from the uncertainty propagation of the model input parameters on the domain wall energy densities. The investigation provides critical insight into the development of parameter subset selection, uncertainty quantification, and propagation methodologies for material modeling domain wall structure evolution, informed by density functional theory simulations.

Keywords

Ferroelectric active subspace analysis uncertainty quantification

1. Introduction

The engineering and development of adaptive structures and intelligent systems is widely supported by the employment of ferroelectric materials as sensors and actuators. Due to their strong electromechanical coupling, high actuation rates, and large relative work densities as solid-state transducers, these materials are ideal for improving the mechanical capabilities and high energy efficiency of many mechatronic devices. Furthermore, applications such as structural health monitoring sensors (Park et al., 2018), energy harvesting circuits (Song et al., 2010), flow control actuators (Bilgen et al., 2011; Kumar et al., 2011), and ambulatory and flying microrobots (Hoffman and Wood, 2011; Pérez-Arancibia et al., 2011; Wood, 2008; Wood et al., 2012) are enabled by the actuator and sensor capabilities of ferroelectric materials, such as lead zirconate titanate (PZT). The wide employment of these materials motivates the need for analytic theoretical methodologies to accurately characterize energy properties ranging in scales from the atomic level to a macroscopic continuum domain.

Among the properties of ferroelectric materials that produce its strong electromechanical behavior are the influences of domain walls separating regions in which polarization is oriented in different directions. These regions separating polarization orientation are fundamental to the governing physical properties and performance of ferroelectric materials. As detailed in Su and Landis (2007), they are, at least partially, responsible for the nonlinearities, hysteresis, and independent dynamics in the constitutive responses at all scales in ferroelectric materials. Moreover, previous investigations have demonstrated a strong coupling between polarization and strain and its dominance over domain wall energies (Cao and Cross, 1991; Meyer and Vanderbilt, 2002; Su and Landis, 2007). This motivates the uncertainty and sensitivity analysis regarding the effects of polarization and strain interactions leading to the structure evolution of ferroelectric domain walls.

Higher order strain gradients, flexoelectricity, and multiferroic coupling may also influence the electromechanical constitutive behavior (Kalinin et al., 2010, 2018; Morozovska et al., 2009, 2017; Oates, 2017; Vasudevan et al., 2012). This can become more prevalent in nanoscale structures where surface area-to-volume ratios are significant. Additional complexities occur in multiferroics such as bismuth ferrite (Balke et al., 2012). The large number of possible material parameters governing polydomain material interactions calls for methods to understand the most important coupling phenomena and their uncertainty. In subsequent paragraphs, we introduce active subspaces to identify the most important parameters that inform ferroelectric continuum phase-field models in light of higher fidelity computations or experiments. Here we will use density functional theory (DFT) calculation to inform the ferroelectric phase-field model.

DFT calculations accurately estimate many structure-dependent properties in the materials at the atomic scale. Nonetheless, DFT calculations are not feasible when solving problems on a continuum scale due to the required large-scale computations. Furthermore, it is difficult to use DFT to fully quantify the properties of PZT. This requires the introduction of phase-field models, characterizing polarization and strain order parameters, to simulate larger scale domain structure evolution, while accurately accounting for underlying atomic and electronic structure properties.

The uncertainty in data and the introduction of the order parameters, which homogenizes effects from the atomic structure, leads to uncertainty in the phase-field models, which can be quantified via the use of Bayesian inference. These methods may be computationally expensive to implement depending on the model evaluation cost and/or the number of unknown model inputs and associated uncertainties. The uncertainty and cost are increased when considering effects across multi-domain structures and domain wall motion.

In this investigation, we consider 180° and 90° polydomain phase-field energy models and apply them to lead titanate. Each phase-field energy model is composed of a stored energy, which is a function of elastic, Landau or polarization, electrostrictive, and polarization gradient energy relations. In each of these relations, attributes of the energy behavior are governed by phenomenological parameters. We investigate the uncertainty and influence of all the unknown parameters except for the elastic coefficients. The elastic coefficient properties influencing the mechanical energy for lead titanate are well known from other investigations (King-Smith and Vanderbilt, 1994) and are therefore fixed in our analysis. We note that phase-field models are often reserved for the time-dependent Ginzburg–Landau equations. In our context, we use the term phase field simply to denote that a differential equation is used to define a diffuse interface between two domain structures. The solution of the 90° domain wall problem does not require the time-dependent term at equilibrium. This is discussed in more detail in Appendix 2.

Prior to employing Bayesian inference for model calibration, it is necessary to determine which model inputs are uniquely inferred by the data, and fix those that are not influential on model responses. For the phase-field model, we focus on two fundamental questions pertaining to the parameters: (1) which gradient energy parameters are identifiable or influential in the sense that they uniquely contribute to 180° and 90° domain wall energies and (2) how influential are Landau polarization and electrostrictive parameters to the polydomain energy responses? The answers to both questions comprise the parameter subset selection, active subspace analysis, and uncertainty quantification and propagation detailed in this article. To answer the first question, we apply these techniques to domain wall energy densities derived based on Euler’s equation. The relations quantify the change in polarization and strain due to the transition across 180° and 90° domain boundaries. To answer the second question, we apply the methods to energy responses in monodomain regions of uniform polarization and strain, in addition to the energy responses in close proximity to the domain walls. The effect of electromechanical coupling of polarization and strain on the domain wall energies is indicated by the sensitivity of Landau energy parameters and electrostrictive coefficients, when considered with the 180° and 90° domain wall energy systems derived in section 2.

Parameter subset selection methods, such as global sensitivity analysis, isolate parameters that are most influential or contribute the most to a model’s output. In addition, Fisher information–based parameter subset selection methods determine identifiable parameters in the sense that they are uniquely inferred by the data. These methods are also advantageous for determining noninfluential or unidentifiable parameters, which can be fixed at nominal values for model calibration and uncertainty propagation. For a more detailed discussion on identifiable and influential parameters, we refer the reader to sources including Leon et al. (2018b) and Smith (2014). In the context of our problem, we use the Fisher information–based parameter subset selection as this requires a smaller number of function evaluations to determine an identifiable subset of parameters.

Parameter subspace methods provide an alternative appealing approach which identifies important relations or linear combinations in the admissible parameter space that most affect the outputs. We refer to these methods as active subspace methods since active variables are determined as weighted linear combinations of the parameters indicating directions of strongest variability. These directions are separated from the ones in which the model response is relatively flat, yielding the inactive variables.

We determine an active subspace for the parameters via the use of gradient approximation methods, ranking the most important directions in the space of all inputs. We then use this active subspace to construct a response surface as a function of the active variables which approximates the original model output response. From this analysis, we then develop global sensitivity metrics, termed activity scores, that rank the relative importance of the model inputs. We compare these results with parameter subset selection results based on the analysis of the Fisher information matrix. Given the influential set of parameters, we compare uncertainty propagation results when only propagating influential parameters with the case when all parameter uncertainties are propagated on 180° and 90° polydomain energy models.

These methods can be generally applied to other phase-field models containing uncertain phenomenological parameters. Recall that the main objectives are to (1) isolate and fix unidentifiable parameters via parameter subset selection techniques, (2) quantify parameter and model uncertainties, and (3) propagate uncertainties to construct credible intervals about quantities of interest. We describe the techniques used to address these objectives in section 3.

This article has the following structure. Section 2 introduces the phase-field energy model, from which we derive governing equations and numerical solutions for the monodomain energy, representing regions of uniform polarization, and the 180° and 90° domain wall energy relations, as given in Appendix 1. Section 3 details the procedures for parameter subset selection, active subspace analysis, Bayesian uncertainty analysis, and uncertainty propagation for the models of interest. We separate our results into two different sections. In section 4, we provide results for the 180° domain wall model analysis, whereas in section 5, we provide the 90° domain wall model analysis results. Finally, in section 6 we provide concluding remarks.

2. Model

In general, phase-field models may be used for characterizing phase transitions in materials, among other phenomena such as ferroelectric hysteresis. We require a model that incorporates polarization and stress-dependent switching behavior across multiple domains, along with the effects of domain wall interactions, motivating the model used in this investigation.

We use a phase-field energy model for characterizing the multi-domain structure evolution of lead titanate, with polarization $P$ and strain $ε$ as the order parameters (see Table 1). In particular, we consider the stored energy

u = u_{M} (ε_{ij}) + u_{P} (P_{i}) + u_{G} (P_{i, j}) + u_{C} (ε_{ij}, P_{i})

(1)

for $i, j = 1, 2, 3$ where

\begin{matrix} u_{M} = \frac{c_{11}}{2} (ε_{11}^{2} + ε_{22}^{2} + ε_{33}^{2}) \\ + c_{12} (ε_{11} ε_{22} + ε_{22} ε_{33} + ε_{11} ε_{33}) \\ + 2 c_{44} (ε_{12}^{2} + ε_{23}^{2} + ε_{13}^{2}) \\ u_{P} = α_{1} (P_{1}^{2} + P_{2}^{2} + P_{3}^{2}) + α_{11} {(P_{1}^{2} + P_{2}^{2} + P_{3}^{2})}^{2} \\ + α_{12} (P_{1}^{2} P_{2}^{2} + P_{2}^{2} P_{3}^{2} + P_{1}^{2} P_{3}^{2}) \\ + α_{111} (P_{1}^{6} + P_{2}^{6} + P_{3}^{6}) \\ + α_{112} [P_{1}^{4} (P_{2}^{2} + P_{3}^{2}) + P_{2}^{4} (P_{1}^{2} + P_{3}^{2}) \\ + P_{3}^{4} (P_{1}^{2} + P_{2}^{2})] + α_{123} P_{1}^{2} P_{2}^{2} P_{3}^{2} \\ u_{C} = - q_{11} (ε_{11} P_{1}^{2} + ε_{22} P_{2}^{2} + ε_{33} P_{3}^{2}) \\ - q_{12} [ε_{11} (P_{2}^{2} + P_{3}^{2}) + ε_{22} (P_{1}^{2} + P_{3}^{2}) \\ + ε_{33} (P_{1}^{2} + P_{2}^{2})] \\ - 4 q_{44} (ε_{12} P_{1} P_{2} + ε_{13} P_{1} P_{3} + ε_{23} P_{2} P_{3}) \\ u_{G} = \frac{g_{11}}{2} (P_{1, 1}^{2} + P_{2, 2}^{2} + P_{3, 3}^{2}) \\ + g_{12} (P_{1, 1} P_{2, 2} + P_{1, 1} P_{3, 3} + P_{2, 2} P_{3, 3}) \\ + \frac{g_{44}}{2} [{(P_{1, 2} + P_{2, 1})}^{2} + {(P_{1, 3} + P_{3, 1})}^{2} \\ + {(P_{2, 3} + P_{3, 2})}^{2}] \end{matrix}

Here, $u_{M}$ , $u_{P}$ , $u_{C}$ , and $u_{G}$ are the mechanical energy, Landau or polarization energy, electrostrictive energy, and gradient energy, respectively, as described by Cao and Cross (1991). In the relations, the unknown parameters $c_{11}, c_{12}$ , and $c_{44}$ are the elastic coefficients, $q_{11}, q_{12}$ , and $q_{44}$ are the electrostrictive coefficients, $α_{1}, α_{11}, α_{12}, α_{111}$ , and $α_{112}$ are the Landau phenomenological parameters, and $g_{11}, g_{12}$ , and $g_{44}$ are the gradient energy exchange parameters. The gradient energy parameters specify the energy potential as a function of the polarization gradient $P_{i, j} = \partial P_{i} / \partial x_{j}$ and are considered only when modeling interactions across domain walls. The convention used in this article for scaling the $(P_{1}^{2} + P_{2}^{2} + P_{3}^{2})$ term in the Landau energy matches Cao and Cross (1991). In the Appendix of Miles et al. (2018), the authors provide a discussion on the impact of this difference for monodomain structures. In the context of the polydomain analysis, the effective coefficients are the same as a result of taking the derivative of the Landau energy with respect to the polarization.

Table 1.

List and definition of symbols used in the article.

Symbol	Definition	Unit
$P$	Polarization order parameter	$C / m^{2}$
$ε$	Strain order parameter	–
$θ_{ϕ}$	$ϕ = 180 °, 90 °$ domain wall energy model parameters	See Tables 4 and 9
$u$ , $u_{ϕ}$	Stored energy, where the subscript $ϕ = 180 °, 90 °$ may indicate $180 °$ or $90 °$ polydomain energy surface, respectively	$mJ / μ m^{3}$
$E_{ϕ}$	$ϕ = 180 °, 90 °$ total domain wall energy	$mJ / m^{2}$

Note that the units are scaled with respect to the values given in Table 5, to facilitate the analysis.

The conditions for static equilibrium for polarization $P$ and strain $ε$ are

\frac{\partial}{\partial x_{j}} (\frac{\partial u}{\partial P_{i, j}}) - \frac{\partial u}{\partial P_{i}} = 0 (i, j = 1, 2, 3)

(2)

σ_{ij, j}^{tot} = \frac{\partial}{\partial x_{j}} (\frac{\partial u}{\partial ε_{ij}}) = 0

(3)

where $σ_{ij}^{tot}$ denotes the Cauchy stress. We employ the equilibrium conditions in the derivation of energy densities for monodomain and polydomain structures. As detailed in Appendix 1, this includes derivation of the monodomain energy, and the 180° and 90° domain wall energy surfaces.

Specifically, we obtain the monodomain energy density $u_{0}$ by substituting uniformly oriented polarization away from the domain walls. In the case of 180° domain walls, we obtain the energy density $u_{180}$ , assuming that polarization changes by 180° as it transitions through the domain wall. In addition, we obtain the 90° domain wall energy density $u_{90}$ , assuming that polarization is oriented perpendicularly with respect to the domain boundary.

We are primarily interested in the total energy responses associated with the domain walls. As these quantities significantly affect material structure–dependent properties, as detailed in the introduction, we denote them as our quantities of interest. In particular, we obtain the total 180° domain wall energy response

E_{180 °} (θ_{180}) = \int_{- \infty}^{\infty} (u_{180} - u_{0}) d x_{1}

(4)

where

θ_{180} = [α_{1}, α_{11}, α_{111}, q_{11}, q_{12}, g_{44}]

(5)

represents the unknown parameters. For the 90° domain wall, we obtain

E_{90 °} (θ_{90}) = \int_{- \infty}^{\infty} (u_{90} - u_{0}) ds

(6)

where

θ_{90} = [α_{1}, α_{11}, α_{12}, α_{111}, α_{112}, q_{11}, q_{12}, q_{44}, g_{11}, g_{12}, g_{44}]

(7)

are the unknown parameters in the 90° domain wall energy. Recall that we provide the derivations and numerical implementation for the energy terms in Appendices 1 and 2. We note that domain wall properties and parameter dependencies may change with respect to different phases in the material (e.g. rhombohedral, orthorhombic). Therefore, this must be accounted for when developing appropriate domain wall energy relations.

2.1. Parameters and distributions

For the rest of this investigation, unless otherwise noted, we use the notation $θ$ to arbitrarily denote either the parameter $θ_{180}$ or $θ_{90}$ given in equations (5) and (7), respectively. Likewise, the notation $θ_{i}$ denotes the $i th$ parameter in $θ_{180}$ or $θ_{90}$ . Note that we will use $θ$ as realizations of the random variable $Θ$ .

The available data for the 180° and 90° total domain wall energies are, respectively, reported to be $E_{180 °}^{*} = 132 mJ / m^{2}$ and $E_{90 °}^{*} = 50 mJ / m^{2}$ from the DFT studies and Landau–Ginzburg theory analysis in Meyer and Vanderbilt (2002) and Stemmer et al. (1995). As both are scalar valued, we employ the characterized uncertainties for the Landau phenomenological parameters and electrostrictive coefficients quantified in previous studies of single-domain structure evolution for PbTiO₃ (Leon et al., 2018b; Miles et al., 2018). The corresponding nominal values and standard deviations from the single-domain uncertainty analysis are presented in Table 2. To determine nominal values for the new exchange parameters $g_{11}, g_{12}$ , and $g_{44}$ , we employ least-squares optimization using synthetic data generated as

E_{180 °}^{data} (i) = E_{180 °}^{*} + ε_{i}, E_{90 °}^{data} (i) = E_{90 °}^{*} + ε_{i}

(8)

where $E^{data}$ are obtained by adding random noise

ε_{i} ~ N (0, σ^{2})

where $σ = 0.05 E_{180 °}^{*}$ or $σ = 0.05 E_{90 °}^{*}$ is selected arbitrarily. The choice of 5% coefficient of variation (CV) was chosen to reflect that these results were obtained from high-fidelity simulations of which we are reasonably confident. It is true that, by changing the amount of variability in our data, one would expect a corresponding change in the uncertainty of the parameter value(s) being calibrated. Whether or not that uncertainty can be attributed to random variability in the experimental or computational procedure is unknown. Furthermore, expanding the CV can lead to poor parameter combinations which causes the numerical solver to not converge. Finally, Meyer and Vanderbilt (2002) state that the domain wall width had 5% uncertainty, which we assume is transferable to the uncertainty in the domain wall energy.

Table 2.

Nominal values and standard deviations for the Landau energy phenomenological parameters and electrostrictive coefficients, obtained from the monodomain uncertainty analysis of Miles et al. (2018).

Parameter	$α_{1}$	$α_{11}$	$α_{12}$	$α_{111}$
Unit	MV m/C	MV m⁵/C³	MV m⁵/C³	MV m⁹/C⁵
Nominal value	−389.2	761.1	414.6	61.61
Standard deviation	10.49	30.01	241.6	19.98
Parameter	$q_{112}$	$q_{11}$	$q_{12}$	$q_{44}$
Unit	MV m⁹/C⁵	GV m/C	GV m/C	GV m/C
Nominal value	−740.5	19.2	3.14	1.40
Standard deviation	499.4	0.258	0.182	0.538

We propose two cases for the consideration of parameter nominal values and distributions employed in section 3. In the first case, we construct independent normal density distributions for each of the parameters in $θ_{180}$ and $θ_{90}$

Θ_{i} ~ N (μ_{i}, σ_{i}^{2})

where $μ_{i}$ and $σ_{i}$ are the nominal values and standard deviations provided in Table 2. In the second scenario, we consider the case where we sample directly from the parameter chains obtained from the Bayesian analysis in Miles et al. (2018). We compare and contrast uncertainty propagation results with respect to both scenarios in sections 4.4 and 5.4, for the $180 °$ and $90 °$ domain wall models, respectively. It is shown in sections 4.4 and 5.4 that the choice of distribution is independent of the parameter subset selection and active subspace analysis results. This supports choosing a normal distribution to facilitate parameter sampling.

Note that nominal values for the elastic coefficients utilized in this investigation are presented in Table 3. These elastic properties are well understood from previous DFT investigations of lead titanate (King-Smith and Vanderbilt, 1994). Furthermore, we neglect uncertainty and sensitivity with respect to the elastic parameters $c_{11}, c_{12}$ , and $c_{44}$ .

Table 3.

Elastic coefficients for PbTiO₃ obtained from King-Smith and Vanderbilt (1994).

Parameter	$c_{11}$	$c_{12}$	$c_{44}$	Unit
Value	3.206e–1	1.397e–1	9.647e–2	$μ J / μ m^{3}$

3. Methodology

Before detailing our parameter subset selection algorithm, Bayesian model calibration, and uncertainty propagation techniques, we provide a summary of the section’s organizational structure. In all of the methods detailed in this section, we employ the parameter scaling techniques described in Appendix 3 to avoid subset selection inaccuracies due to large differences in parameter magnitudes. This yields the normalized parameters $x$ , as detailed in the Appendix.

We employ this information to construct normalized derivatives used in the gradient-based parameter subset selection and active subspace approaches described in sections 3.1 and 3.2, respectively. Finally, in section 3.3, we introduce Bayesian inference and uncertainty propagation methods used in model calibration and verification of parameter subset selection results.

3.1. Parameter subset selection

Here we introduce the parameter subset selection method used to determine the most influential and/or identifiable parameters. We consider first a method based on local sensitivity and then expand this to one which addresses global sensitivity analysis based on the construction of elementary effects for determining sensitivity measures (Morris, 1991).

We start by considering the Fisher information matrix

F = S^{T} S

(9)

which contains local sensitivity information utilized to determine locally noninfluential or nonidentifiable parameters. Now, consider the scalar-valued response

y_{n} = f (x_{1}^{n}; θ), n = 1, \dots, N

where $x_{1}^{n}$ are values centered around $x_{1} = 0$ . This response could be, for example, values of the 180° domain wall energy density $u_{180} (x_{1}) - u_{0}$ defined in equation (4).

For some data $Y_{n}$ , we consider the functional

J (θ) = \frac{1}{N} \sum_{n = 1}^{N} {[Y_{n} - f (x_{1}^{n}; θ)]}^{2}

(10)

Minimization of equation (10) yields an optimal vector $θ^{*}$ of parameter values. Recall that $θ$ is locally identifiable at $θ^{*}$ if the minimum is uniquely determined by the data.

We consider the multivariate Taylor expansion

f (x_{1}^{n}; θ) \approx f (x_{1}^{n}; θ^{*}) + \nabla_{θ} f (x_{1}^{n}; θ^{*}) \cdot Δ θ

where

\nabla_{θ} f (x_{1}^{n}; θ^{*}) = {[\frac{\partial f}{\partial θ_{1}} (x_{1}^{n}; θ^{*}), \dots, \frac{\partial f}{\partial θ_{p}} (x_{1}^{n}; θ^{*})]}^{T}

(11)

and $Δ θ = θ - θ^{*}$ . The cost functional (equation (10)) can then be approximated by

J (θ) \approx \frac{1}{N} \sum_{n = 1}^{N} {[\nabla_{θ} f (x_{1}^{n}; θ^{*}) \cdot Δ θ]}^{2}

based on the assumption that

Y_{n} \approx f (x_{1}^{n}; θ^{*})

at the minimum $θ^{*}$ . Defining the $N \times p$ sensitivity matrix $S$ by

S = [\begin{matrix} \frac{\partial f}{\partial θ_{1}} (x_{1}^{1}; θ^{*}) & \dots & \frac{\partial f}{\partial θ_{p}} (x_{1}^{1}; θ^{*}) \\ ⋮ & ⋱ & ⋮ \\ \frac{\partial f}{\partial θ_{1}} (x_{1}^{N}; θ^{*}) & \dots & \frac{\partial f}{\partial θ_{p}} (x_{1}^{N}; θ^{*}) \end{matrix}]

(12)

we approximate

J (θ) \approx \frac{1}{N} (S Δ θ)^{T} (S Δ θ)

such that $J (θ^{*} + Δ θ) \approx (1 / N) Δ θ^{T} S^{T} S Δ θ$ . Letting $S^{T} S Δ θ = λ Δ θ$ , where we take $Δ θ$ to be an eigenvector of $S^{T} S$ , yields

J (θ^{*} + Δ θ) \approx \frac{1}{N} λ ‖ Δ θ ‖_{2}^{2}

The cost functional perturbations are approximately $0$ when $λ \approx 0$ , making the corresponding parameters locally nonidentifiable. The algorithms in Brun et al. (2002), Burth et al., (1999), Cintrón-Arias et al., (2009), and Quaiser and Mönnigmann (2009) are thus based on this methodology and we implement the algorithm by Quaiser and Monnigmann (2009), as presented in Algorithm 1, in our investigation.

Algorithm 1. Parameter subset selection algorithm to determine locally unidentifiable parameters (Quaiser and Mönnigmann, 2009).
(0) Set $η = p$ , where $p$ is the number of parameters in the model, and construct the sensitivity matrix $S$ following relation (12). We note that the variable $η$ changes with the iterations of the algorithm. (1) Compute the matrix $S^{T} S$ and its eigenvalues, and order their magnitudes as $\| λ_{1} \| \leq \| λ_{2} \| \leq \dots \leq \| λ_{η} \|$ (2) If $\| λ_{1} \| > ϵ$ , where $ϵ$ is some prescribed threshold value, stop. We take all the parameters to be identifiable. (3) If $\| λ_{1} \| < ϵ$ , then one of the parameters is not identifiable. Proceed as follows. (4) Identify the component with the largest magnitude in the eigenvector $Δ θ_{1}$ associated with $λ_{1}$ . This component corresponds to the least identifiable parameter. (5) Remove the column in $S$ corresponding to the component identified in Step 4, set $η = η - 1$ , and repeat Step 1.

Algorithm 1. Parameter subset selection algorithm to determine locally unidentifiable parameters (Quaiser and Mönnigmann, 2009).

(0) Set

η = p

, where

p

is the number of parameters in the model, and construct the sensitivity matrix

S

following relation (12). We note that the variable

η

changes with the iterations of the algorithm.
(1) Compute the matrix

S^{T} S

and its eigenvalues, and order their magnitudes as

| λ_{1} | \leq | λ_{2} | \leq \dots \leq | λ_{η} |

(2) If

| λ_{1} | > ϵ

, where

ϵ

is some prescribed threshold value, stop. We take all the parameters to be identifiable.
(3) If

| λ_{1} | < ϵ

, then one of the parameters is not identifiable. Proceed as follows.
(4) Identify the component with the largest magnitude in the eigenvector

Δ θ_{1}

associated with

λ_{1}

. This component corresponds to the least identifiable parameter.
(5) Remove the column in

S

corresponding to the component identified in Step 4, set

η = η - 1

, and repeat Step 1.

3.1.1. Global parameter subset selection

Recall that Algorithm 1 relies on local sensitivity information for identifying influential parameters. Here, we propose a methodology that explores the parameter space for $θ$ , to obtain a measure for global sensitivity, while employing the underlying theory behind Algorithm 1.

Consider $ρ (θ)$ to be the underlying distribution for parameters $θ$ . This density could be the normal distribution such that

Θ ~ N (μ, V)

where $μ$ is the vector of nominal values and $V$ is the corresponding covariance matrix for $θ$ . Sampling from distribution $ρ (θ)$ , we obtain $θ^{ℓ}$ and define the sampled local sensitivity entries

s_{ij}^{ℓ} = \frac{\partial f}{\partial θ_{i}} (x_{1}^{j}; θ^{ℓ}), ℓ = 1, \dots, K

where $K$ is the number of samples. Similar to the Morris screening global sensitivity analysis (Morris, 1991), these sensitivity entries are computed as elementary effects

\begin{matrix} s_{ij}^{ℓ} = \frac{f (x_{1}^{i}; θ_{1}^{ℓ}, \dots, θ_{j}^{ℓ} + Δ, θ_{j + 1}^{ℓ}, \dots, θ_{p}^{ℓ}) - f (x_{1}^{i}; θ^{ℓ})}{Δ} \\ = \frac{f (x_{1}^{i}; θ^{ℓ} + Δ \cdot e_{j}) - f (x_{1}^{i}; θ^{ℓ})}{Δ} \end{matrix}

(13)

where $Δ$ is the step size and $e_{j}$ is the $j th$ canonical vector.

To obtain the global sensitivity matrix, we average over all the sampled sensitivity entries to obtain the sensitivity measures

s_{ij}^{*} = \frac{1}{K} \sum_{ℓ = 1}^{K} s_{ij}^{ℓ}

This yields the updated global sensitivity matrix

S^{*} = [\begin{matrix} s_{11}^{*} & \dots & s_{1 p .}^{*} \\ ⋮ & ⋱ & ⋮ \\ s_{N 1}^{*} & \dots & s_{Np}^{*} \end{matrix}]

(14)

We employ this global-based sensitivity matrix along with Algorithm 1 for the parameter subset selection investigation in this article.

3.2. Active subspace

We define the active subspace for a scalar quantity of interest with respect to some parameters following the derivation presented in Leon et al. (2018a). Given $Y = f (Θ)$ , we compute the matrix

C = E [(\nabla_{θ} f) (\nabla_{θ} f)^{T}]

(15)

The matrix $C$ is symmetric and semi-positive definite by construction. Thus, it has the real eigenvalue decomposition

C = W Λ W^{T}

for $Λ = diag (λ_{1}, λ_{2}, \dots, λ_{p}), λ_{1} \geq \dots \geq 0$ . Based on a possible gap in the eigenvalue spectrum, we consider the partitions $W = [W_{1}, W_{2}]$ and

Λ = [\begin{matrix} Λ_{1} \\ Λ_{2} \end{matrix}], Λ_{1} = diag (λ_{1}, \dots, λ_{m}), m < p

Using this significant gap, we define the new rotated coordinate variables

y = W_{1}^{T} θ, z = W_{2}^{T} θ

where $y \in R^{m}$ and $z \in R^{p - m}$ . Based on this transformation, on average the output $f$ varies more dominantly due to variability and perturbations in the directions dictated by $y$ than in the directions dictated by $z$ (Constantine, 2015), defining the active and inactive subspaces.

In practice, construction of the matrix (equation (15)) may require computing high-dimensional integrals. Therefore, based on Constantine (2015) and Bang (2012), we approximate the gradient matrix

G = \frac{1}{\sqrt{M}} [\begin{matrix} \nabla_{θ} f_{1} & \nabla_{θ} f_{2} & \dots & \nabla_{θ} f_{M} \end{matrix}]

(16)

We then use the singular value decomposition (SVD) to extract an input active subspace basis. The basis for the active subspace is contained in the matrix of singular vectors $U$ in

G = U Σ V^{T}

(17)

Note that the diagonal matrix $Σ$ contains the squared roots of the eigenvalues of $C$ . This methodology for computing the active subspace is equivalent to the methodology employed by the formation of the matrix $C$ as shown by Russi (2010). The first $m$ columns of $U$ thus contain the most important directions in the parameter space.

3.2.1. Construction

Using the scaling maps and derivative terms from Appendix 3, we approximate the gradient matrix $G$ (16). The main idea is to first select normalized initial vectors $x^{j} \in R^{p}$ , $j = 1, \dots, M$ , sampled from a probability density $ρ (x)$ , where $x$ are the scaled parameters. In our case, we take $ρ (x)$ to be the normal distribution from equation (67). We then construct a set of coarse derivative approximations, termed elementary effects

\begin{matrix} d_{i} (x^{j}) = \frac{h (x_{1}^{j}, \dots, x_{i - 1}^{j}, x_{i}^{j} + Δ, x_{i + 1}^{j}, \dots, x_{p}^{j}) - h (x^{j})}{Δ} \\ = \frac{h (x^{j} + Δ \cdot e_{i}) - h (x^{j})}{Δ} \end{matrix}

given a step size $Δ$ , for the $i th$ parameter and $j th$ sample point drawn from $ρ (x)$ . Here

h (x) = f (Ax + θ^{nom}) = f (θ)

where $f (θ)$ could represent either $E_{180 °} (θ_{180})$ or $E_{90 °} (θ_{90})$ . Recall that $θ$ denotes the original input space parameters.

We use this method, as presented in Algorithm 2, to construct the columns of the gradient matrix for the active subspace construction.

Algorithm 2. Finite difference method gradient approximation (Lewis et al., 2016).
(0) Identify the number of desired columns $M$ , appropriate step size $Δ$ , and probability density $ρ (x)$ . Use the model function $h (x)$ for normalized parameters as defined in Appendix 3. for $j = 1, \dots, M$ (1) Construct the row vector $D_{1 \times p}$ with randomly chosen entries $\pm Δ$ , where $Δ$ is the step size. (2) Select $x^{j}$ from the probability density $ρ (x)$ . In our case this is the normal distribution. Evaluate the model function at the sampled value along with the elementary effect as follows. for $i = 1, \dots, p$ $d_{i} = \frac{h (x_{i}^{j} + D (i) \cdot e_{i}) - h (x^{j})}{D (i)}$ , where $e_{i}$ is the $i th$ standard basis vector. end (3) Let $\hat{G} (:, j) = d$ , where $\hat{G}$ is the gradient matrix approximation. end

Algorithm 2. Finite difference method gradient approximation (Lewis et al., 2016).

(0) Identify the number of desired columns

M

, appropriate step size

Δ

, and probability density

ρ (x)

. Use the model function

h (x)

for normalized parameters as defined in Appendix 3.
for

j = 1, \dots, M

(1) Construct the row vector

D_{1 \times p}

with randomly chosen entries

\pm Δ

, where

Δ

is the step size.
(2) Select

x^{j}

from the probability density

ρ (x)

. In our case this is the normal distribution. Evaluate the model function at the sampled value along with the elementary effect as follows.
for

i = 1, \dots, p

d_{i} = \frac{h (x_{i}^{j} + D (i) \cdot e_{i}) - h (x^{j})}{D (i)}

, where

e_{i}

is the

i th

standard basis vector.
end
(3) Let

\hat{G} (:, j) = d

, where

\hat{G}

is the gradient matrix approximation.
end

3.2.2. Determining the dimension of the active subspace

Several criteria must be considered to determine and verify the active subspace dimension. The most obvious method for identifying the dimension is to find large gaps in the singular value spectrum. This is achieved by employing SVD of the gradient matrix $G$ . Large gaps are usually noticed based on a difference of one or two orders of magnitude between the first and second or second and third singular values, thereby implying a one-dimensional (1D) or a two-dimensional (2D) active subspace, respectively. To verify results of the gap-based method, we employ several dimension selection methods. First, we verify the gap-based method via the construction of a response surface for different possible active subspace dimensions. We also employ an error-based dimension selection method associated with a user-defined tolerance, as well as principal component analysis (PCA), analyzing the gradient matrix $G$ , as described in Appendix 4. We present the outline of the response surface construction next.

3.2.3. Response surface construction

The second method for active subspace dimension selection consists of constructing a response surface that approximates the original function output. We construct these for possible dimensions of the active subspace, as functions of the active variable $y$ defined earlier in this section. We use, for example, multivariate polynomial regression to construct a response surface denoted by $r (y)$ . The procedure we implement for constructing the response surface is presented in Algorithm 3. To avoid overfitting the data, we employ the Akaike information criteria (AIC)

AIC = 2 k + M \cdot \log (\frac{1}{M} SSR)

(18)

where $k$ is the number of parameters, $M$ is the number of data points $q_{i}$ used in fitting to the data, and

SSR = \sum_{i = 1}^{N} {[q_{i} - r (y_{i})]}^{2}

is the residual sum of squares. The AIC penalizes the number of parameters in the polynomial, while rewarding a maximized likelihood function. For the different polynomial degrees considered, we choose the model with the minimum AIC score.

Algorithm 3. Algorithm for constructing response surface based on active subspace.
(0) Sample the training input values $x_{i}$ with respect to its probability density function and construct the corresponding responses $q_{i} = h (x_{i})$ . (1) Project the sampled values $x_{i}$ onto the active subspace using the transformation $y = W_{1}^{T} x_{i}$ . (2) Using regression analysis, construct a response surface $r (y)$ , using $q_{i}$ , such that $q_{i} \approx r (y_{i})$ . (3) Verify the response surface by approximating the original function $h (x) \approx r (W_{1}^{T} x)$ .

Algorithm 3. Algorithm for constructing response surface based on active subspace.

(0) Sample the training input values

x_{i}

with respect to its probability density function and construct the corresponding responses

q_{i} = h (x_{i})

.
(1) Project the sampled values

x_{i}

onto the active subspace using the transformation

y = W_{1}^{T} x_{i}

.
(2) Using regression analysis, construct a response surface

r (y)

, using

q_{i}

, such that

q_{i} \approx r (y_{i})

.
(3) Verify the response surface by approximating the original function

h (x) \approx r (W_{1}^{T} x)

To quantitatively compare the response surface fit to the actual solution, we employ the mean relative error (MRE)

MRE = \frac{1}{M} \sum_{i = 1}^{M} | \frac{q_{i} - r (y_{i})}{q_{i}} |

We choose the appropriate active subspace dimension when there is less than one order of magnitude in the reduction of MRE, as we increase the number of dimensions $m$ . To validate the dimension of the active subspace, we also compare this method with the dimension selection techniques detailed in Appendix 4.

3.2.3. Activity scores

To obtain a sensitivity metric to rank parameter influence, we compute global sensitivity measures based on the active subspace. As detailed in Constantine and Diaz (2017), these metrics are based on the formula

α_{i} = α_{i} (m) = \sum_{j = 1}^{m} λ_{j} w_{i, j}^{2}, i = 1, \dots, p

(19)

where $m$ is the possible dimension of the active subspace, $λ_{j}$ , $j = 1, \dots, m$ are the first $m$ eigenvalues, and $w_{i, j}$ are the corresponding eigenvector entries. We use these to rank the relative importance of each of the model’s inputs and compare with the results of parameter subset selection.

3.3. Bayesian inference and uncertainty propagation

We employ Bayesian inference techniques to assess parameter uncertainty when the domain wall energy models defined in section 2 are calibrated with respect to data obtained from experimental measurements or high-fidelity model simulations. The data values $E_{180}^{*}$ and $E_{90}^{*}$ , provided in section 2.1 and obtained from DFT simulations in Meyer and Vanderbilt (2002) and Stemmer et al. (1995), are considered high-fidelity results in our investigation. To accommodate the underlying uncertainty, model parameters $Θ$ are assumed to be random variables in the Bayesian framework.

Here, we focus on the uncertainty of the polydomain gradient energy exchange parameters $g_{11}, g_{12}$ , and $g_{44}$ . In the case of $g_{44}$ , we compare the total domain wall energy DFT data point $E_{180 °}^{*}$ to the model response $E_{180 °} (θ_{180})$ . To infer the parameter $g_{44}$ , we fix all other parameters in $θ_{180}$ . Namely, we fix the Landau phenomenological parameters $α_{1}, α_{11}, α_{111}$ , and the electrostrictive coefficients $q_{11}, q_{12}$ at the nominal values in Table 2.

In the case of the parameters $g_{11}$ and $g_{12}$ , we compare the total domain wall energy DFT data point $E_{90 °}^{*}$ to the model response $E_{90 °} (θ_{90})$ . Similar to the $180 °$ domain wall energy model, we fix all other parameters in $θ_{90}$ to nominal values presented in Table 2, whereas the parameter $g_{44}$ is fixed at a nominal value obtained through a least-squares optimization of equation (4). We employ the parameter subset selection analysis detailed in section 3.1 to determine whether $g_{11}$ or $g_{12}$ is most influential. Based on this analysis, we infer the uncertainty of the most influential parameter, while fixing the other at a nominal value obtained through a least-squares optimization of equation (6).

We let $M^{obs}$ denote synthetic measurement observations for either the total 180° or 90° domain wall energies. This yields the statistical model

M^{obs} (s) = f (Θ) + ε_{s}, s = 1, \dots, n

(20)

where $ε_{s} ~ N (0, σ_{obs}^{2})$ are assumed independent and identically distributed (iid) errors induced by the continuum model $f (Θ)$ , denoting either $E_{180 °} (Θ_{180})$ (equation (4)) or $E_{90 °} (θ_{90})$ (equation (6)).

For Bayesian model calibration, we employ Bayes’ equation

π (θ | m^{obs}) = \frac{π (m^{obs} | θ) π_{0} (θ)}{\int_{R^{p}} π (m^{obs} | θ) π_{0} (θ) d θ}

(21)

to infer the probability densities of the gradient exchange parameters, based on the synthetic data

m^{obs} (s) = m_{α}^{*} + ε_{s}, s = 1, \dots, n

(22)

In equation (22), $m_{α}^{*}$ is either of the literature values $E_{180 °}^{*}$ or $E_{90 °}^{*}$ , and $ϵ_{s}$ is random noise such that $ϵ_{s} \sim N (0, σ_{o b s}^{2})$ . In this investigation, we take $σ_{obs} = 0.05 m_{α}^{*}$ so that the data are normally distributed with 5% standard deviation about the literature values $E_{180 °}^{*}$ or $E_{90 °}^{*}$ . The posterior density $π (θ | m^{obs})$ quantifies the probability of observing parameter values $θ$ given data denoted by $m^{obs}$ . In addition, the prior density $π_{0} (θ)$ contains prior knowledge about the distribution of parameters $θ$ , before evaluation of the likelihood $π (m^{obs}, θ)$ . The denominator in equation (21) is a normalization constant to obtain a posterior density which integrates to unity.

For $ε_{s} ~^{iid} N (0, σ_{obs}^{2})$ , the likelihood function is

π (m^{obs} | θ) = e^{- \sum_{s = 1}^{N} {[m^{obs} (s) - f (θ)]}^{2} / 2 σ^{2}}

(23)

We employ sampling-based Markov chain Monte Carlo (MCMC) algorithms to compute the posterior distribution. Specifically, we implement the delayed rejection adaptive Metropolis (DRAM) algorithm, developed by Haario et al. (2001, 2006) and detailed in Smith (2014). The code for the corresponding implementation is available on M Laine’s (2013) website http://helios.fmi.fi/∼lainema/mcmc/.

3.3.1. Uncertainty propagation

Once we quantify uncertainties for the gradient exchange parameters, we propagate them by sampling from the posterior densities and computing realizations of the model $f (θ)$ . With these realizations, we construct $50 %, 90 %, 95 %$ , and $99 %$ credible intervals which quantify model fit accuracy. Along with the posterior densities obtained for the exchange parameters, we propagate the uncertainties for the Landau phenomenological parameters and electrostrictive coefficients, for which the nominal values and standard deviations are given in Table 2.

We consider two cases for the sampling-based uncertainty propagation of these parameters. In the first scenario, we assume that the parameters are each normally distributed and

Θ_{i} ~ N (μ_{i}, σ_{i}^{2})

where the means and standard deviations are given in Table 2. In the second scenario, we sample the values for the parameters directly from the chain of values obtained in single-domain structure analysis of Miles et al. (2018). The second scenario incorporates correlation in the monodomain parameters obtained from previous analyses. Since the polydomain energy densities also incorporate spontaneous polarization far from the domain walls, it is important to also consider these effects.

We additionally compare the uncertainty propagation of all the parameters on the model response $f (θ)$ with the case when only the influential parameters are propagated. These subsets of parameters are determined based on the analysis of section 3.1 on parameter subset selection. This analysis illustrates the use of parameter subset selection for determining which parameters may be fixed in anticipation of uncertainty propagation and future Bayesian uncertainty analysis with respect to longitudinal data. To compare the credible intervals obtained when employing uncertainty propagation in both cases, we employ energy statistics (Székely and Rizzo, 2013), as detailed in Appendix 5.

4. 180° domain wall model analysis

As detailed previously, $180 °$ twinned domain walls separate oppositely oriented polarization domains. In section 4.1, we illustrate the solution to the energy density as the $180 °$ domain wall is crossed. In addition, in sections 4.2 and 4.3, we present parameter subset selection and active subspace analysis results. Finally, we employ information obtained from these sections in the uncertainty propagation analysis of section 4.4.

4.1. Model solution

We first establish the accuracy of the solution of the $180 °$ domain wall energy system. Due to the assumptions made in the modified analytic solution of Cao and Cross (1991) (see Appendix 1), we verify its accuracy by comparing the results of the domain wall energy density $u_{180} - u_{0}$ , polarization $P_{3} (x_{1})$ , polarization gradient $P_{3, 1} (x_{1})$ , and strain $ε_{11} (x_{1})$ with the results obtained when implementing the numerical procedure and MATLAB solver bvp4c.m presented in Appendix 2. We plot the results in Figure 1. We conclude that the two methods are consistent for the solution $P_{3} (x_{1})$ and $ε_{11} (x_{1})$ , with the integration (equation (4)) of both solutions yielding $E_{180 °} \approx 132 mJ / m^{2}$ .

Figure 1.

Comparison of analytic, finite difference, and MATLAB bvp4c.m solver for solution of the $180 °$ domain wall energy system. From top left to bottom right, $180 °$ domain wall energy density solution, polarization $P_{3} (x_{1})$ , strain $ε_{11} (x_{1})$ , and polarization gradient $P_{3, 1} (x_{1})$ as we cross the $180 °$ domain wall.

To obtain a nominal value for $g_{44}$ , we employ a least-squares optimization with respect to the synthetic data defined in equation (8). The compiled set of nominal values for $θ_{180}$ is presented in Table 4. The corresponding measurement unit system is presented in Table 5. This system will serve as the one used in the implementation of both the $180 °$ and $90 °$ domain wall energy models.

Table 4.

Nominal values for the parameters $θ_{180}$ with respect to the domain wall energy model $E_{180 °} (θ_{180})$ .

$θ_{180}$	$α_{1}$	$α_{11}$	$α_{111}$	$q_{11}$	$q_{12}$	$g_{44}$
Nominal value	−389.4e6	761.1e6	61.46e6	19.2e6	3.14e6	5.73e–6
Unit	$\frac{mV mm}{μ C}$	$\frac{mV m m^{5}}{μ C^{3}}$	$\frac{mV m m^{9}}{μ C^{5}}$	$\frac{mV mm}{μ C}$	$\frac{mV mm}{μ C}$	$\frac{mV m m^{3}}{μ C}$

Table 5.

Units of measurement used in the total energy density with respect to the polydomain (PD) models.

Measurement	Energy	Length	Charge	Voltage
Unit (PD)	nJ/mm³	mm	µC	mV

4.2. 180° domain wall model parameter subset selection

Next, we perform parameter subset selection for the $180 °$ domain wall model. We implement the method outlined in section 3.1. Note that the density used in the random sampling of parameters $θ_{180}$ for the construction of the global sensitivity matrix (equation (14)) is the one defined in equation (67). For the parameters $α_{1}, α_{11}, α_{111}, q_{11}, q_{12}$ , we use the nominal values and standard deviations given in Table 2, scaled with respect to the unit system given in Table 5. We obtain the corresponding nominal value and standard deviation for $g_{44}$ by employing Bayesian inference and characterizing the parameter’s uncertainty.

Employing Algorithm 1, we find that the most influential parameters in the $180 °$ domain wall energy model are $α_{11}, α_{111}, q_{11}$ , and $g_{44}$ when using a threshold $ε = 1 \times 10^{- 9}$ , one order of magnitude below the square root of machine epsilon. The threshold value was chosen due to the variability of the smallest eigenvalue in the third iteration of the algorithm. This variability occurs as a result of the random sampling process used in the construction of the sensitivity matrix. Thus, choosing $1 \times 10^{- 9}$ provides a better threshold for accurately choosing the identifiable parameters. Note that we used a total of $500$ realizations of the sensitivity entries $s_{ij}^{ℓ}$ in the construction of the global sensitivity matrix (equation (14)). Our results for the analysis are given in Table 6.

Table 6.

Results from Algorithm 1 with the global sensitivity matrix (equation (14)) to determine noninfluential parameters in $θ_{180}$ (equation (5)) for the 180° domain wall energy $E_{180 °}$ (equation (4)).

Iteration	$\| λ_{1} \|$	Eigenvector $Δ θ_{1}$ with the corresponding parameters
		$α_{1}$	$α_{11}$	$α_{111}$	$q_{11}$	$q_{12}$	$g_{44}$
1	6.22e–13	8.01e–2	−8.20e–2	1.13e–2	−3.64e–1	−9.24e–1	−9.14e–3
2	7.55e–12	9.33e–1	−3.08e–1	8.64e–2	−6.01e–2	–	−1.51e–1
3	2.62e–9	–	−5.43e–1	2.01e–2	−7.12e–1	−	−4.45e–1
Result: The parameters $α_{1}$ and $q_{12}$ are not influential.

4.3. 180° domain wall model active subspace determination

To obtain an active subspace for the model $E_{180 °}$ , we first employ Algorithm 2 to construct the gradient matrix (equation (16)). We obtained a solution using $M = 20, 000$ iterations. In the construction of the active subspace, we assumed a Gaussian distribution of the type (equation (67)) for $θ_{180}$ . The corresponding singular values are presented in Figure 2, where the shading corresponds to two standard deviations. Based on the singular value spectrum, it is difficult to identify a dimension for the active subspace. To determine whether this is a 1D or 2D active subspace, we proceed with the other methods detailed in section 3.2.

Figure 2.

Singular values obtained from the active subspace determination for $E_{180 °}$ assuming normal distributions (equation (67)) for parameters $θ_{180}$ . The shading indicates two standard deviations from the sample mean.

We construct response surfaces for $E_{180 °}$ considering 1D or 2D active subspaces. Employing Algorithm 3, we obtain the response surfaces presented in Figure 3. For a quantitative comparison, we consider the MREs with respect to the two response surfaces. We present the MREs in Table 7. Recall that the active subspace dimension is indicated when the decrease in MRE is less than one order of magnitude. In this case, we therefore conclude that the MRE indicates a 1D active subspace.

Figure 3.

Response surfaces for $E_{180 °}$ constructed based on (a) a one-dimensional and (b) a two-dimensional active subspace.

Table 7.

Mean relative errors (MREs) for response surfaces considering one-, two-, and three-dimensional active subspaces for the mode response $E_{180 °}$ .

No. of dimensions	1	2	3
MRE	0.0244	0.0116	0.0116

The next method we consider for dimension selection is PCA. Taking $t^{*} = 0.95$ , the employment of Algorithm 5 yields an active subspace dimension equal to 1. Both the response surface errors and PCA indicate a 1D active subspace.

In computing activity scores to rank the relative importance of the model inputs $θ_{180}$ on the response $E_{180 °}$ , we consider 1D and 2D active subspaces. We obtain the activity scores presented in Figure 4. Based on the activity scores, the parameters $α_{11}, α_{111}$ , and $q_{11}$ are the most sensitive parameters. This is in agreement with the results from the parameter subset selection of section 4.2, except for the gradient exchange parameter $g_{44}$ . Note however that in the present analysis we consider the total energy associated with the $180 °$ domain wall, whereas in section 4.2 we consider the energy density function across the domain wall $u_{180} (x_{1}) - u_{0}$ , as presented in Figure 1. The contribution of the exchange parameter is thus more apparent in the energy density as opposed to the domain wall energy $E_{180 °}$ .

Figure 4.

Activity scores for the model response $E_{180 °}$ assuming (a) a one-dimensional and (b) a two-dimensional active subspace.

4.4. 180° domain wall model uncertainty propagation

From the results of the parameter subset selection and activity scores, we assume that the parameters $α_{11}, α_{111}, q_{11}$ , and $g_{44}$ are the most influential. We then denote the sensitive parameters as

θ_{180}^{sens} = [α_{11}, α_{111}, q_{11}, g_{44}]

(24)

As noted previously, we consider two cases for uncertainty propagation. In the first case, we consider the parameters to be normally distributed

Θ_{180} ~ N (μ_{180}, V_{180})

(25)

where $μ_{180}$ contains the nominal values presented in Table 4 and $V_{180}$ is defined as a diagonal covariance matrix containing the standard deviation squared, from Table 2. In the second case, we sample the parameters strictly from the chains obtained in the monodomain analysis of Miles et al. (2018).

In both cases, we construct $50 %, 90 %, 95 %$ , and $99 %$ credible intervals and compare these results with the credible intervals obtained when only propagating the uncertainties of the sensitive parameters $θ_{180}^{sens}$ , while fixing the other parameters. We present the results of uncertainty propagation considering the parameters to be normally distributed in Figure 5(a) and (b) for all parameters $θ_{180}$ and only sensitive parameters $θ_{180}$ , respectively, propagated in the model.

Figure 5.

Uncertainty propagation of the (a) $θ_{180}$ model inputs and (b) $θ_{180}^{sens}$ influential inputs on the energy density $u_{180} (x_{1}) - u_{0}$ . The nominal solution indicates the model evaluated using parameter values from Table 4.

To analyze the credible intervals shown in Figure 5, we consider the “peak values” (at $x_{1} = 0$ ) and subtract the sensitive parameter uncertainty propagation from the full parameter uncertainty propagation. Sampling from the distribution defined in equation (25), as well as from the chains obtained in Miles et al. (2018), we obtain the distributions for $u_{180} (0) - u_{0}$ plotted in Figure 6. From these plots, we observe that the distributions of the credible intervals are very similar. A quantitative comparison is however needed in order to determine if the distributions of the obtained credible intervals are indeed sufficiently similar.

Figure 6.

Density of the $180 °$ domain wall energy $u_{180} (0) - u_{0}$ distributions with respect to the case where the uncertainties of all parameters $θ_{180}$ are propagated compared with the case where only the uncertainties of parameters $θ_{180}^{sens}$ are propagated. In (a), parameters were sampled from the normal distribution (equation (25)), whereas in (b) parameters were sampled directly from the chains in the monodomain analysis of Miles et al. (2018).

This motivates the use of energy statistics to test the null hypothesis $H_{0} : F_{180} = F_{180}^{sens}$ that the two densities are equal to each other. Implementing the energy statistics methodology outlined in Appendix 5, we obtain, for $M = 499$ and $n_{1} = n_{2} = 4000$ , the energy test statistic $T_{n_{1}, n_{2}} = 0.2413$ . Constructing the replicates in the manner detailed in Appendix 5, we obtain the results presented in Figure 7. For the confidence levels $α = 0.05, 0.10$ , we present the energy statistics results in Table 8. As is observed from Figure 7 and Table 8, the test statistic is clearly below the critical values dictated by the confidence levels. As a result, we reject the alternate hypothesis that the distributions are not equal. This verifies that the parameters $θ_{180}^{sens}$ are the most influential in the model.

Figure 7.

Energy test statistic replicates for the analysis of the $u_{180} (0) - u_{0}$ distribution considering (a) normal distribution and (b) direct sampling from the monodomain chains of Miles et al. (2018) for the Landau and electrostrictive parameters in $θ_{180}$ .

Table 8.

Energy test statistic and critical values for $α = 0.05$ , $α = 0.10$ with respect to normally distributed $θ_{180}$ (normal) and sampled from the chains in the monodomain (MD) analysis of Miles et al. (2018).

Density peak samples (normal)	Test statistic	Critical value $(α = 0.05)$	Critical value $(α = 0.10)$
$θ_{180}$ versus $θ_{180}^{sens}$ sampled	0.2413	1.5761	1.2728
Density peak samples (MD)	Test statistic	Critical value $(α = 0.05)$	Critical value $(α = 0.10)$
$θ_{180}$ versus $θ_{180}^{sens}$ sampled	0.3428	0.6353	0.4919

5. 90° domain wall model analysis

The $90 °$ domain walls separate perpendicularly oriented polarization domains. As detailed in section “90° domain wall energy” of Appendix 1, we facilitate the system analysis by defining a new coordinate system rotated by $45 °$ . In section 5.1, we provide solutions to the system based on this coordinate transformation. We present our parameter subset selection results for determining influential parameters in $θ_{90}$ (equation (7)) in section 5.2. In addition, in section 5.3 we present our obtained active subspace for $E_{90 °}$ . Finally, we verify our results with the uncertainty propagation of the influential subset of parameters in section 5.4.

5.1. Model solution

In this section, we present the solutions of the $90 °$ domain wall energy model system. We compare the results of the root finding problem procedure, given in section “90° domain wall numerical solution” of Appendix 2, with the results obtained when implementing MATLAB’s solver bvp4c.m. Solutions we obtain for the domain wall energy density $u_{90} (s) - u_{0}$ , polarization $P_{r} (s), P_{s} (s)$ , and polarization gradients $P_{r, s}, P_{s, s}$ are presented in Figure 8.

Figure 8.

Comparison of finite difference fsolve.m and bvp4c.m MATLAB solvers for the solution of the $90 °$ domain wall energy system. From top left to bottom right, $90 °$ domain wall energy density solution, polarization $P_{s}$ in the $s$ direction, polarization $P_{r}$ in the $r$ direction, and polarization gradients $P_{s, s}$ and $P_{r, s}$ in the $s$ direction as we cross the $90 °$ domain wall.

To obtain nominal values for $g_{11}$ and $g_{12}$ , we employ a least-squares optimization with respect to the synthetic data defined in equation (8). The compiled set of nominal values for $θ_{90}$ is presented in Table 9.

Table 9.

Nominal values for parameters $θ_{90}$ with respect to the domain wall energy model $E_{90 °} (θ_{90})$ .

$θ_{90}$	$α_{1}$	$α_{11}$	$α_{12}$	$α_{111}$	$α_{112}$	$q_{11}$
Nominal value	−389.4e6	761.1e6	414.1e6	61.46e6	−740.56e6	19.2e9
Unit	mV mm/µC	mV mm⁵/µC³	mV mm⁵/µC³	mV mm⁹/µC⁵	mV mm⁹/µC⁵	mV mm/µC
$θ_{90}$	$q_{12}$	$q_{44}$	$g_{11}$	$g_{12}$	$g_{44}$
Nominal value	3.14e9	1.40e9	1.60e–5	−2.17e–8	5.73e–6
Unit	mV mm/µC	mV mm/µC	mV mm³/µC	mV mm³/µC	mV mm³/µC

5.2. 90° domain wall model parameter subset selection

Next we implement the parameter subset selection methodology detailed in section 3.1, to determine the influential parameters in equation (7). We construct the global sensitivity matrix following the same methodology used in section 4. For the exchange parameters $g_{11}$ and $g_{12}$ , we assume normal distributions about the nominal values $μ_{g_{11}}, μ_{g_{12}}$ provided in Table 9

G_{11} ~ N (μ_{g_{11}}, (0.05 | μ_{g_{11}} |)^{2}), G_{12} ~ N (μ_{g_{12}}, (0.05 | μ_{g_{12}} |)^{2})

For the parameter $g_{44}$ , we use the nominal value and standard deviation identified in section 4 to construct the respective normal distribution.

We present the parameter subset selection results for the $90 °$ domain wall model in Table 10. As $g_{12}$ is the first parameter eliminated in the algorithm, we quantify the uncertainty in parameter $g_{11}$ while fixing all other parameters including $g_{12}$ with respect to the model response (equation (6)).

Table 10.

Results from Algorithm 1 with the global sensitivity matrix (equation (14)) to determine noninfluential parameters in $θ_{90}$ (equation (7)) for the 90° domain wall energy $E_{90 °}$ (equation (6)).

Iteration	$\| λ_{1} \|$	Eigenvector $Δ θ_{1}$ with the corresponding parameters
		$α_{1}$	$α_{11}$	$α_{12}$	$α_{111}$	$α_{112}$	$q_{11}$	$q_{12}$	$q_{44}$	$g_{11}$	$g_{12}$	$g_{44}$
1	3.30e–14	−2.57e–2	−9.67e–2	−8.97e–2	5.02e–3	9.95e–3	−3.25e–1	−6.27e–1	−3.36e–2	−7.87e–2	−6.68e–1	1.69e–1
2	6.91e–14	−5.91e–2	1.27e–1	−7.65e–2	−1.25e–2	−4.84e–2	−2.43e–2	−4.62e–1	2.65e–1	2.70e–1	−	−7.84e–1
3	1.21e–12	−6.64e–1	−1.37e–1	1.15e–1	−3.90e–1	5.66e–2	−5.07e–1	1.15e–1	−1.22e–1	−2.95e–1	−	−
4	5.86e–12	−	5.06e–1	−2.82e–1	4.75e–2	−4.53e–2	7.55e–1	5.49e–3	7.94e–2	2.90e–1	−	−
5	3.31e–11	−	−2.84e–2	5.39e–1	3.41e–1	−2.03e–2	−	4.67e–2	3.03e–1	6.76e–1	−	−
6	6.18e–10	−	1.76e–1	−3.85e–1	−1.10e–2	1.49e–2	−	−9.06e–1	7.16e–3	−	−	−
7	6.77e–9	−	2.39e–1	−3.71e–1	−6.93e–1	6.94e–2	−	−	5.65e–1	−	−	−
Result: The parameters $α_{1}, q_{11}, q_{12}, g_{11}, g_{12}$ , and $g_{44}$ are not influential based on $ϵ = 1 \times 10^{- 9}$ .

5.3. 90° domain wall model active subspace determination

To obtain an active subspace for the model $E_{90 °}$ , we construct the gradient matrix (equation (16)), using $M = 10, 000$ iterations in the implementation of Algorithm 2. In the construction of the active subspace, we assumed a normal distribution of the type (equation (67)) for $θ_{90}$ . The corresponding singular values are presented in Figure 9. The singular value spectrum motivates employing additional methods to accurately determine the dimension of the active subspace. We determine an appropriate active subspace dimension following the methods outlined in Appendix 4.

Figure 9.

Singular values obtained from the active subspace determination for $E_{90 °}$ assuming normal distributions (equation (67)) for parameters $θ_{90}$ . The shading around the singular values corresponds to two standard deviations.

Employing Algorithm 3, we obtain the response surface plotted in Figure 10. To evaluate an active subspace dimension, we consider the MREs presented in Table 11, for possible 1D and 2D response surfaces.

Figure 10.

Response surface for $E_{90 °}$ constructed based on a one-dimensional active subspace.

Table 11.

Mean relative errors (MREs) for response surfaces considering one-, two-, three-, and four-dimensional active subspaces for the model response $E_{90 °}$ .

No. of dimensions	1	2	3	4
MRE	0.0104	0.0075	0.0061	0.0061

In the PCA, the employment of Algorithm 5 yields an active subspace dimension equal to 1. Both the response surface errors and PCA agree on a 1D active subspace.

Considering a 1D active subspace, we obtain the activity scores presented in Figure 11. Based on the activity scores, the parameters $α_{1}, α_{11}, α_{12}, α_{111}, α_{112}, q_{11}, q_{12}$ , and $g_{11}$ are the most sensitive parameters, whereas those corresponding to activity scores less than $1 \times 10^{- 2}$ were determined to be noninfluential. This is in agreement with the results from the parameter subset selection of section 5.2, except for the Landau parameter $α_{1}$ and electrostrictive coefficients $q_{11}$ and $q_{12}$ . Note, however, that in the present analysis we consider the total energy associated with $90 °$ domain wall defined as $E_{90 °}$ , whereas in section 5.2 we consider the energy density function across the domain wall $u_{90} (s) - u_{0}$ as presented in Figure 8. The contribution of the electrostrictive coefficient $q_{11}$ is thus more apparent in the total energy $E_{90 °}$ compared to the energy density function.

Figure 11.

Activity scores for the model response $E_{90 °}$ assuming a one-dimensional active subspace.

5.4. 90° domain wall model uncertainty propagation

From the results of the parameter subset selection and activity scores, we conclude that the parameters $α_{1}, α_{11}, α_{12}, α_{111}, α_{112}, q_{11}, q_{12}, q_{44}$ , and $g_{11}$ are the most influential. We denote the sensitive parameters as

θ_{90}^{sens} = [α_{1}, α_{11}, α_{12}, α_{111}, α_{112}, q_{11}, q_{12}, q_{44}, g_{11}]

(26)

Recall that we consider two cases for uncertainty propagation. In the first case, we consider the parameters to be normally distributed such that

Θ_{90} ~ N (μ_{90}, V_{90})

(27)

where $μ_{90}$ contains the nominal values presented in Table 9 and $V_{90}$ is defined as a diagonal covariance matrix containing the standard deviation squared from Table 2. The standard deviations for $g_{11}$ and $g_{44}$ are obtained by employing Bayesian inference, whereas we assume a standard deviation $σ_{g_{12}} = 0.05 | μ_{g_{12}} | °$ for $g_{12}$ . In the second case, we sample the Landau and electrostrictive parameters strictly from the chains obtained in the monodomain analysis of Miles et al. (2018).

Repeating the analysis of section 4.4, we obtain $50 %, 90 %, 95 %$ , and $99 %$ credible intervals when considering the parameters to be normally distributed, similar to the plots presented in Figure 5. To compare results when propagating all parameters $θ_{90}$ and only sensitive $θ_{90}^{sens}$ , we analyze the difference in peak values of the $90 °$ domain wall uncertainty propagation. This yields the kernel density estimation (kde) of the distributions for $u_{90} (0) - u_{0}$ for both sampling cases, plotted in Figure 12. To quantitatively compare the densities, we again employ energy statistics. We test the null hypothesis $H_{0} : F_{90} = F_{90}^{sens}$ that the two densities are equal to each other. We obtain, for $M = 499$ and $n_{1} = n_{2} = 4000$ , the energy test statistic $T_{n_{1}, n_{2}} = 1.2881 \times 10^{- 2}$ and the replicates plotted in Figure 13. These correspond to the results in Table 12, for confidence levels $α = 0.05, 0.10$ . As is observed, the test statistic is clearly below the critical values dictated by the confidence levels. As a result, we reject the alternate hypothesis that the distributions are not equal. This verifies that the parameters $θ_{90}^{sens}$ are the most influential in the $90 °$ domain wall energy model. For future model calibration, we can fix the insensitive parameters while only inferring the influential ones.

Figure 12.

Density of the $90 °$ domain wall energy $u_{90} (0) - u_{0}$ distributions with respect to the case where the uncertainties of all parameters $θ_{90}$ are propagated as compared with the case where only the uncertainties of parameters $θ_{90}^{sens}$ are propagated. In (a), parameters were sampled from the normal distribution (equation (27)), whereas in (b) we sampled the Landau and electrostrictive parameters directly from the chains in Miles et al. (2018).

Figure 13.

Energy test statistic replicates for the statistical analysis of the $u_{90} (0) - u_{0}$ distribution, with $θ_{90}$ being normally distributed.

Table 12.

Energy test statistic and critical values for $α = 0.05$ and $α = 0.10$ , when sampling the Landau and electrostrictive parameters from a normal distribution (normal) and directly from the monodomain (MD) chains (Miles et al., 2018).

Density peak samples (normal)	Test statistic	Critical value $(α = 0.05)$	Critical value $(α = 0.10)$
$θ_{90}$ versus $θ_{90}^{sens}$ sampled	1.2881e–2	2.9544e–1	2.1138e–1
Density peak samples (MD)	Test statistic	Critical value $(α = 0.05)$	Critical value $(α = 0.10)$
$θ_{90}$ versus $θ_{90}^{sens}$ sampled	2.1899e–2	1.4273e–1	1.1512e–1

6. Concluding remarks

In this investigation, we have performed a parameter subset selection and uncertainty analysis and propagation study of $180 °$ and $90 °$ domain wall energies. We have considered the Ginzburg–Landau phase-field model applied to ferroelectric domain structures for the material lead titanate. Our analysis shows the importance of Landau phenomenological parameters and electrostrictive coefficients relative to polarization gradient energy associated with domain walls. In addition, we have provided the foundations for a study in which longitudinal data could be used to inform domain wall energy densities to more uniquely infer parameters at the polydomain scale.

Future work will consist in employing a simultaneous analysis of $180 °$ and $90 °$ domain walls. Extensions to other phases (e.g. rhombohedral, orthorhombic) and their associated domain structures could provide important insights on domain engineered materials with unique functionality. This is of particular interest due to the complex motion and interactions between the domain walls. A maximum entropy approach, as proposed by Gao et al. (2017), may be used to address uncertainty analysis challenges posed by data fusion from different sources informing the polydomain energy models. The approach is motivated by accurate weighting mechanisms, supported by the covariance structure between datasets.

Footnotes

Appendix 1

Appendix 2

Appendix 3

Appendix 4

Appendix 5 Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The research of L.S.L. was based upon the work supported by the National Science Foundation (NSF) under Grant No. DGE-1633587. The research of L.S.L. and R.C.S. was also supported in part by the NSF Grant CMMI-1306290 Collaborative Research CDS&E, whereas the research of P.R.M. and W.S.O. was supported in part by the NSF Grant CMMI-1306320 Collaborative Research CDS&E.

ORCID iDs

Lider S Leon

William S Oates

References

Balke

Winchester

Ren

et al. (2012) Enhanced electric conductivity at ferroelectric vortex cores in BiFeO₃. Nature Physics 8(1): 81–88.

Bang

Abdel-Khalik

Hite

(2012) Hybrid reduced order modeling applied to nonlinear models. International Journal for Numerical Methods in Engineering 91(9): 929–949.

Bilgen

De Marqui

Jr Kochersberger

et al. (2011) Macro-fiber composite actuators for flow control of a variable camber airfoil. Journal of Intelligent Material Systems and Structures 22(1): 81–91.

Brun

Kühni

Siegrist

et al. (2002) Practical identifiability of ASM2d parameters—systematic selection and tuning of parameter subsets. Water Research 36(16): 4113–4127.

Burth

Verghese

Vélez-Reyes

(1999) Subset selection for improved parameter estimation in on-line identification of a synchronous generator. IEEE Transactions on Power Systems 14(1): 218–225.

Cao

Cross

(1991) Theory of tetragonal twin structures in ferroelectric perovskites with a first-order phase transition. Physical Review B 44(1): 5–12.

Cintrón-Arias

Banks

Capaldi

et al. (2009) A sensitivity matrix based methodology for inverse problem formulation. Journal of Inverse and Ill-Posed Problems 17(6): 545–564.

Constantine

(2015) Active Subspaces: Emerging Ideas for Dimension Reduction in Parameter Studies, Vol. 2. Philadelphia, PA: SIAM.

Constantine

Diaz

(2017) Global sensitivity metrics from active subspaces. Reliability Engineering & System Safety 162: 1–13.

10.

Gao

Oates

Smith

(2017) A maximum entropy approach for uncertainty quantification and analysis of multifunctional materials. In: Proceedings of the conference on smart materials, adaptive structures and intelligent systems, Snowbird, UT, 18–20 September, p. V001T08A013. New York: American Society of Mechanical Engineers.

11.

Haario

Laine

Mira

et al. (2006) DRAM: efficient adaptive MCMC. Statistics and Computing 16(4): 339–354.

12.

Haario

Saksman

Tamminen

et al. (2001) An adaptive Metropolis algorithm. Bernoulli 7(2): 223–242.

13.

Halko

Martinsson

Tropp

(2011) Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Review 53(2): 217–288.

14.

Hoffman

Wood

(2011) Myriapod-like ambulation of a segmented microrobot. Autonomous Robots 31(1): 103.

15.

Jolliffe

(2002) Principal Component Analysis (Springer Series in Statistics). New York, NY, USA: Springer.

16.

Kalinin

Kim

Fong

et al. (2018) Surface-screening mechanisms in ferroelectric thin films and their effect on polarization dynamics and domain structures. Reports on Progress in Physics 81(3): 036502.

17.

Kalinin

Morozovska

Chen

et al. (2010) Local polarization dynamics in ferroelectric materials. Reports on Progress in Physics 73(5): 056502.

18.

King-Smith

Vanderbilt

(1994) First-principles investigation of ferroelectricity in perovskite compounds. Physical Review B 49(9): 5828.

19.

Kumar

Hays

Fernandez

et al. (2011) Flow sensitive actuators for micro-air vehicles. Smart Materials and Structures 20(10): 105033.

20.

Laine

(2013) DRAM website. Available at: http://helios.fmi.fi/∼lainema/dram/

21.

Leon

Smith

Miles

et al. (2018a) Active subspace uncertainty quantification for a polydomain ferroelectric phase-field model. In: Proceedings of the behavior and mechanics of multifunctional materials and composites XII, Vol. 10596, p. 105960T. International Society for Optics and Photonics. Available at: https://www.spiedigitallibrary.org/conference-proceedings-of-spie/10596/105960T/Active-subspace-uncertainty-quantification-for-a-polydomain-ferroelectric-phase-field/10.1117/12.2297207.short?SSO=1

22.

Leon

Smith

Oates

et al. (2018b) Analysis of a multi-axial quantum-informed ferroelectric continuum model: part 2—sensitivity analysis. Journal of Intelligent Material Systems and Structures 29: 2840–2860.

23.

Lewis

Smith

Williams

(2016) Gradient free active subspace construction using Morris screening elementary effects. Computers & Mathematics with Applications 72(6): 1603–1615.

24.

Meyer

Vanderbilt

(2002) Ab initio study of ferroelectric domain walls in pbTIO₃. Physical Review B 65(10): 104111.

25.

Miles

Leon

Smith

et al. (2018) Analysis of a multi-axial quantum informed ferroelectric continuum model: part 1—uncertainty quantification. Journal of Intelligent Material Systems and Structures 29: 2823–2839.

26.

Miles

Oates

Leon

et al. (2017) Uncertainty analysis of ferroelectric polydomain structures. In: Proceedings of the conference on smart materials, adaptive structures and intelligent systems (SMASIS), Snowbird, UT, 18–20 September, p. V001T02A008. New York: American Society of Mechanical Engineers.

27.

Morozovska

Eliseev

et al. (2009) Thermodynamics of nanodomain formation and breakdown in scanning probe microscopy: Landau-Ginzburg-Devonshire approach. Physical Review B 80(21): 214110.

28.

Morozovska

Kalinin

Eliseev

(2017) Flexoelectricity impact on the domain wall structure and polar properties. In: Morozovska

Kalinin

Eliseev

(eds) Flexoelectricity in Solids: From Theory to Applications. Hackensack, NJ: World Scientific, pp. 311–336.

29.

Morris

(1991) Factorial sampling plans for preliminary computational experiments. Technometrics 33(2): 161–174.

30.

Oates

(2017) Flexoelectricity, strain gradients, and singularities in ferroelectric nanostructures. Journal of Intelligent Material Systems and Structures 28(20): 3091–3105.

31.

Park

Palazotto

Hale

et al. (2018) Internal longitudinal damage detection in a steel beam using lamb waves: simulation and test study. Journal of Intelligent Material Systems and Structures 29: 411–422.

32.

Pérez-Arancibia

Galloway

et al. (2011) First controlled vertical flight of a biologically inspired microrobot. Bioinspiration & Biomimetics 6(3): 036009.

33.

Quaiser

Mönnigmann

(2009) Systematic identifiability testing for unambiguous mechanistic modeling—application to JAK-STAT, MAP kinase, and NF-kappaB signaling pathway models. BMC Systems Biology 3(1): 50.

34.

Russi

(2010) Uncertainty Quantification with Experimental Data and Complex System Models. Berkeley, CA: University of California.

35.

Smith

(2014) Uncertainty Quantification: Theory, Implementation, and Applications. Philadelphia, PA, USA: Society for industrial and applied mathematics.

36.

Song

Choi

Wereley

et al. (2010) Energy harvesting devices using macro-fiber composite materials. Journal of Intelligent Material Systems and Structures 21(6): 647–658.

37.

Stemmer

Streiffer

Ernst

et al. (1995) Atomistic structure of 90 domain walls in ferroelectric pbTIO3 thin films. Philosophical Magazine A 71(3): 713–724.

38.

Landis

(2007) Continuum thermodynamics of ferroelectric domain evolution: theory, finite element implementation, and application to domain wall pinning. Journal of the Mechanics and Physics of Solids 55(2): 280–305.

39.

Székely

Rizzo

(2013) Energy statistics: a class of statistics based on distances. Journal of Statistical Planning and Inference 143(8): 1249–1272.

40.

Vasudevan

Morozovska

Eliseev

et al. (2012) Domain wall geometry controls conduction in ferroelectrics. Nano Letters 12(11): 5524–5531.

41.

Wood

(2008) The first takeoff of a biologically inspired at-scale robotic insect. IEEE Transactions on Robotics 24(2): 341–347.

42.

Wood

Finio

Karpelson

et al. (2012) Progress on “pico” air vehicles. The International Journal of Robotics Research 31(11): 1292–1302.