Problem 1 Statement
We show in Example 3.5 that the boundary value problem d2Ts =2(a+b)h[Ts(x)−Tamb]
MA 540 — Project 2 Answer Key February 22, 2018
dTs(0)= Φ , dTs(L)= h[Tamb −Ts(L)]
Copyright By PowCoder代写 加微信 powcoder
models the steady state temperature Ts(x) of an uninsulated rod with source heat flux Φ at x = 0 and ambient air temperature Tamb. The model parameters are q = [Φ,h,k], where h is the convective heat transfer coefficient and k is the thermal conductivity.
The analytic solution is
y(x, q) = Ts(x, q) = c1(q)e−γx + c2(q)eγx + Tamb,
where γ = 2(a+b)h and abk
c1(q)=−kγ e−γL(h−kγ)+eγL(h+kγ)
Φ eγL(h+kγ) Φ
accuracy. If you have time, compute and compare the analytic sensitivities for h and k.
(b) For the parameters q = [Φ,h,k], construct the sensitivity matrix χij(q) = ∂y(xi,q) and the
, c2(q)= kγ +c1(q).
We suppress the parameter dependence of γ to clarify the notation. In Chapter 7, we will estimate the parameter values k = 2.37, h = 0.00191 and Φ = −18.4, which you can use here. The measured ambient room temperature is Tamb = 21.29oC and the rod has cross-sectional dimensions a = b = 0.95 cm and length L = 70 cm.
(a) Use finite-differences to approximate the sensitivity relations ∂ y , ∂ y and ∂ y and plot your ∂Φ ∂h ∂k
solutions at the 15 equally spaced spatial locations xi = x0 + (i − 1)∆x, where x0 = 10 cm and
∆x = 4 cm. Use a discrete line-type so your plot looks like Figure 7.4. Additionally, compute
the analytic sensitivity relation ∂y and plot with your finite-difference solution to compare their ∂Φ
matrix V = χT χ. Compute the rank of V and discuss the identifiability of the complete parameter
set. Discuss why you can deduce this result based on the model.
(c) As detailed in the text, the thermal conductivity k is well-documented for aluminum and copper.
Fix this parameter and repeat your analysis for the parameters q = [Φ, h]. Are they identifiable? Solution
We use a centered-difference approximation to compute the numerical partial derivatives, shown in
Figure 1. Because ∂y is orders of magnitude larger than the other two, we show it on a separate ∂h
Figure 1: Local sensitivities for Φ and k (left) and h (right).
We evaluate the analytical sensitivity for Φ using Maple. Before substituting in the nominal
parameter values, the expression for ∂y spanned several lines. After substituting in the nominal ∂Φ
parameter values k = 2.37, h = 0.00191, and Φ = −18.4, we arrive at the equation
∂y(x, q) = −0.1228 exp(−0.0583x + 4.0776) − 0.1194 exp(0.0583x − 4.0776) . (1)
∂Φ q=qnom
In Figure 2 compare the analytical sensitivity in (1) to the numerical partial derivative evaluated
on the x-grid. We observe that they agree very closely, as we would hope.
Figure 2: The local sensitivity of y to Φ, computed both numerically and analytically.
To construct V = χT χ, we re-use the numerical sensitivities from above. The i-th column of χ will be the sensitivity of the i-th parameter evaluated on the x-grid. In this case, columns 1, 2, and 3
correspond to the sensitivities evaluated on the x-grid for Φ, h, and k, respectively. We find that
4.4380e+01 4.2905e+05 −1.2159e+00 V = 4.2905e+05 4.4464e+09 −2.5236e+05 ,
−1.2159e+00 −2.5236e+05 1.9394e+02
where e ± n represents 10±n. Using the MATLAB rank command, we obtain rank(V ) = 2. To confirm this, we find the eigenvalues of V , which are −4.2613 × 10−9, 1.8260 × 102, 4.4464 × 109. Using the standard tolerance of 10−8, we indeed confirm that one eigenvalue of V is numerically 0, so V has rank 2.
Since rank(V ) < 3, then not all of the parameters are mutually identifiable. Indeed, this is reason- able because none of Φ, k, h appear in the model by themselves. Each appearance of a parameter is either in multiplication, division, addition, or subtraction with another parameter.
We now fix k = 2.37 (for aluminum) and recompute V ̃ using only the numerical sensitivities of Φ and h. Doing so, we find
̃ 4.4380e+01 4.2905e+05 V = 4.2905e+05 4.4464e+09 ,
which has rank 2. Therefore, with k fixed at a nominal value, Φ and h are now identifiable.
Problem 2 Statement
Consider the Helmholtz energy
ψ(P,q)=α1P2 +α11P4 +α111P6,
where P is the polarization and q = [α1, α11, α111] are parameters. You can take nominal values to be α1 = −389.4, α11 = 761.3 and α111 = 61.5. When sampling for global sensitivity analysis, you should sample each parameter from U(0,1) and map to the interval [αl,αr] that is 20% above and below the nominal value. For uniformly sampling on the interval [a, b], one would use the command q = a + (b-a)*rand(1,1).
(a) Plot the energy for P in the interval [−0.8, 0.8]. Do you observe the double-well behavior?
(b) Analytically compute the sensitivity matrix χ and matrix V = χT χ using 17 equally spaced polarization values in the domain [0, 0.8]. Compute the rank of V and discuss the identifiability of the parameters q.
(c) Use Morris screening with forward differences and r = 50,∆ = 1/20, to compute μ∗i and σi2. You can use the scalar response
which you can compute analytically. To check your solution, show that μ∗ ≈ ∂y , ∂y , ∂y .
ψ(P, q)dP, (2)
∂α1 ∂α11 ∂α111 Explain why you would expect σi2 ≈ 0 for a linearly parameterized problem such as this one. Which
parameter is least influential?
(d) Use the Saltelli algorithm 15.10.1 to approximate the Sobol sensitivity indices Si and STi . Show
Si ≈ 1. i=1
What does this indicate about the second-order effects Sij and is this to be expected for a linearly parameterized problem? Use kde.m to plot a kernel density estimate of yA computed in Step 3 of the algorithm. Now fix any noninfluential parameters at their nominal values and recompute yA. Plot the new kde on the same plot as the 3-parameter case and discuss your results and determine if there is a discrepancy with (b). We will revisit this problem when we do Bayesian inference.
Figure 3 shows the double-well behavior that we expect from the Helmholtz energy.
Figure 3: Helmholtz energy on [−0.8,0.8] using nominal parameters α1 = −389.4, α11 = 761.3, and α111 = 61.5.
We compute the analytical sensitivites of ψ as
∂ψ =P2, ∂ψ =P4, ∂ψ =P6.
∂α1 ∂α11 ∂α111
Notice that Pi = 0.05(i − 1), i = 1, 2, . . . , 17. Using the analytical sensitivities and the convention
q1 = α1,q2 = α11,q3 = α111, we get
χij =(0.05(i−1))2j, i=1,2,...,17, j=1,2,3.
It follows that
1.5241 0.7384 0.3891 V = 0.7384 0.3891 0.2154 ,
0.3891 0.2154 0.1232
which has rank 3. We verified the rank by computing the eigenvalues of V , the smallest of which is on the order of 10−4. Since there are three parameters and rank(V ) = 3, then α1, α11, α111 are mutually identifiable.
For Morris screening, since we draw each Qi from a dilated and shifted uniform distribution, we mag- nify ∆ by the same dilation factor. Specifically, elementary effects di, we use ∆ ̃ i = (1/20)(0.4)qnom
since the sampling interval for Qi is [0.8qnom,1.2qnom]. We analytically compute the response ii
function as
0.83 0.85 0.87
y(α1,α11,α111) = 3 α1 + 5 α11 + 7 α111 . (3)
We show our results for r = 50 realizations and a zero-threshold of 10−8 in Table 1.
α1 α11 α111 μ∗i 0.1707 0.0655 0.0300 σi2 0 0 0
Table 1: Morris screening with r = 50 realizations.
To confirm our computed values of μ∗i , we calculate the analytical sensitivities of y:
≈ 0.1707, ≈ 0.0655, ≈ 0.0300,
= 3 ∂y 0.85
= 5 ∂y 0.87
so μ∗ ≈ ∂y , ∂y , ∂y . Furthermore, because (3) depends linearly on each of the parameters,
dji = f(qj + ∆ei) − f(qj) ∆
∂α1 ∂α11 ∂α111
is constant for fixed i. Hence, σi2 = 0 for each i = 1, 2, 3 since a line’s slope does not change based on where it is measured. Since σi2 = 0 for each i and μ∗1 > μ∗2 > μ∗3, then q3 = α111 is the least influential parameter according to Morris screening.
For the Sobol analysis, we use 10,000 realizations since each realization is inexpensive. Using the Saltelli algorithm, we arrive at the approximate Sobol indices in Table 2.
α1 Si 0.6457 STi 0.6309
α11 0.3624 0.3495
α111 0.0084 −0.0089
Table 2: Sobol indices for α1, α11, α111. We find that 3i=1 Si = 1.0166 ≈ 1. Since
Si+ Sij=1, i=1 1≤i