\documentclass[11pt]{article}
\usepackage{color} %used for font color
\usepackage{amssymb} %maths
\usepackage{amsmath} %maths
\usepackage{epsfig}
\topmargin-.5in
\textwidth6.6in
\textheight9in
\oddsidemargin0in
\def\ds{\displaystyle}
\def\d{\partial}
\def\bH{{\bf H}}
\def\bM{{\bf M}}
\def\bB{{\bf B}}
\def\bm{{\bf m}}
\def\bS{{\bf S}}
\def\hQ{{\Phi}}
\def\ve{{\varepsilon}}
%\def\fs{\scriptsize}
\def\fs{\footnotesize}
\def\reals{\mathbb{R}}
\def\complex{\mathbb{C}}
\def\scriptO{{{\it O}\kern -.42em {\it `}\kern + .20em}}

\begin{document}

\centerline{\large \bf Math 540: Project 2}

\bigskip
\noindent
Due Tuesday, February 26


\bigskip
\noindent
1. We show in Example~3.5 that the boundary value problem
$$
\begin{array}{l}
  {\ds \frac{d^2 T_s}{d x^2} = \frac{2(a+b)}{ab} \frac{h}{k} \left[T_s(x) - T_{amb} \right] } \\
  \noalign{\medskip}
  {\ds \frac{d T_s}{d x}(0) = \frac{\hQ}{k} }  \quad , \quad
  {\ds \frac{d T_s}{d x}(L) = \frac{h}{k} [T_{amb} - T_s(L)] }
\end{array}
$$
models the steady state temperature $T_s(x)$ of an uninsulated rod with source heat flux $\hQ$ 
at $x=0$ and ambient air temperature $T_{amb}$.  The model parameters are $q = [\hQ,h,k]$, where $h$ is the convective heat transfer coefficient and $k$ is the thermal conductivity.

The analytic solution is
$$
  y(x,q) = T_s(x,q) = c_1(q) e^{-\gamma x} + c_2(q) e^{\gamma x} + T_{amb},
$$
where $\gamma = \sqrt{\frac{2(a+b)h}{abk}}$ and
$$
  c_1(q) = -\frac{\hQ}{k\gamma}\left[\frac{e^{\gamma L}(h+k\gamma)}{e^{-\gamma L}(h-k\gamma) + e^{\gamma L}(h+k\gamma)}\right] \quad , \quad c_2(q) = \frac{\hQ}{k\gamma} + c_1(q).
$$
We suppress the parameter dependence of $\gamma$ to clarify the notation.  In Chapter~7, we will estimate the parameter values $k = 2.37, h = 0.00191$ and $\Phi = -18.4$, which you can use here.  The measured ambient room temperature is $T_{amb} = 21.29^o$C and the rod has cross-sectional dimensions $a = b = 0.95$~cm and length $L = 70$~cm.

\medskip
\noindent
(a) Use finite-differences to approximate the sensitivity relations $\frac{\d y}{\d \Phi}, \frac{\d y}{\d h}$ and $\frac{\d y}{\d k}$ and plot your solutions at the 15 equally spaced spatial locations $x_i = x_0 + (i-1)\Delta x$, where $x_0 = 10$~cm and $\Delta x = 4$~cm.  Use a discrete line-type so your plot looks like Figure 7.4.  Additionally, compute the analytic sensitivity relation $\frac{\d y}{\d \Phi}$ and plot with your finite-difference solution to compare their accuracy.  If you have time, compute and compare the analytic sensitivities for $h$ and $k$.

\medskip
\noindent
(b) For the parameters $q = [\Phi,h,k]$, construct the sensitivity matrix $\chi_{ij}(q) = \frac{\d y(x_i,q)}{\d q_j}$ and the matrix $V = \chi^T \chi$.  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.

\medskip
\noindent
(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 = [\Phi,h]$.  Are they identifiable?


\bigskip
\bigskip
\noindent
2. Consider the Helmholtz energy
$$
  \psi(P,q) = \alpha_1 P^2 + \alpha_{11}P^4 + \alpha_{111}P^6,
$$
where $P$ is the polarization and $q = [\alpha_1, \alpha_{11}, \alpha_{111}]$ are parameters.  You can take nominal values to be $\alpha_1 = -389.4, \alpha_{11} = 761.3$ and $\alpha_{111} = 61.5$.  When sampling for global sensitivity analysis, you should sample each parameter from ${\cal U}(0,1)$ and map to the interval $[\alpha_\ell, \alpha_r]$ that is 20\% above and below the nominal value.  For uniformly sampling on the interval $[a,b]$, one would use the command \texttt{q = a + (b-a)*rand(1,1).}

\medskip
\noindent
(a) Plot the energy for $P$ in the interval $[-0.8, 0.8]$.  Do you observe the double-well behavior?

\medskip
\noindent
(b) Analytically compute the sensitivity matrix $\chi$ and matrix $V = \chi^T \chi$ 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$.

\medskip
\noindent
(c) Use Morris screening with forward differences and $r = 50, \Delta = 1/20$, to compute $\mu_i^*$ and $\sigma^2_i$.  You can use the scalar response
\begin{equation}
  y(q) = \int_0^{0.8} \psi(P,q)dP,
\label{eq1}
\end{equation}
which you can compute analytically. 
To check your solution, show that $\mu^* \approx \left[\frac{\d y}{\d \alpha_1}, \frac{\d y}{\d \alpha_{11}}, \frac{\d y}{\d \alpha_{111}}\right]$.  Explain why you would expect $\sigma_i^2 \approx 0$ for a linearly parameterized problem such as this one.  Which parameter is least influential?

\medskip
\noindent
(d) Use the Saltelli algorithm 15.10.1 to approximate the Sobol sensitivity indices $S_i$ and $S_{T_i}$.  Show that
$$
  \sum_{i=1}^3 S_i \approx 1.
$$
What does this indicate about the second-order effects $S_{ij}$ and is this to be expected for a linearly parameterized problem?  Use \texttt{kde.m} to plot a kernel density estimate of $y_A $ computed in Step 3 of the algorithm.  Now fix any noninfluential parameters at their nominal values and recompute $y_A$.  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.


\bigskip
\bigskip
\noindent
3. Exercise 15.5: {\em Do only the Saltelli-Sobol analysis.}  You can use the following commands to numerically solve the ode and approximate the integral.  Here {\tt A} is the matrix in the Saltelli algorithm.

\begin{verbatim}
  tf = 5;
  dt = 0.01;
  t_data = 0:dt:tf;
  Y0 = [S0; I0; R0];
  M = 1000;
  alpha = 0.2;
  beta = 15;
  A = rand(M,4);
  A(:,2) = betarnd(alpha,beta,M,1);
  
  ode_options = odeset('RelTol',1e-6);
  for j=1:M
      params = A(j,:);
      [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,params);
      y_A(j) = sum(dt*Y(:,3));
  end
\end{verbatim}

\noindent
The associated function is

\begin{verbatim}
%
%          SIR_rhs
%
  function dy = SIR_rhs(t,y,params);

  N = 1000;

  gamma = params(1);
  k = params(2);
  r = params(3);
  mu = params(4);

  dy = [mu*N - mu*y(1) - gamma*k*y(2)*y(1);
        gamma*k*y(2)*y(1) - (r + mu)*y(2);
        r*y(2) - mu*y(3)];
\end{verbatim}

\end{document}

