\documentclass[12pt,letterpaper]{article}
\usepackage{amsmath,amssymb,verbatim}
\usepackage{epsfig}

\def\ds{\displaystyle}
\def\d{\partial}
\newcommand{\dspace}{\vskip 24pt}
\newcommand{\sspace}{\vskip 12pt}
\newcommand{\hfspace}{\vskip 6pt}
\setlength{\textwidth}{6.5in}
\setlength{\oddsidemargin}{0.0in}
\setlength{\textheight}{9in}
\setlength{\topmargin}{-0.5in}
\DeclareMathOperator{\diag}{diag} 


\begin{document}

\pagestyle{empty}
\raggedright

\begin{center}
  {\Huge 
    Numerical Methods for Inverse Problems\\[2.0ex]
   \LARGE HW 3
  }
\dspace

\end{center}

\dspace
The objective of this project is to help students familiarize
themselves with the concepts of paramer estimation problems. In addition,
students will practice using MATLAB to carry out many computations
associated with solving inverse problems.  Read the model description, then
work \underline{either} 1 or 2 ({\bf Bonus: do both!})


\bigskip

\begin{list}{}{}
\item Consider the following mathematical model for a
spring-mass-dashpot system (using a linear spring assumption,
Hooke's law, and viscous air damping):
\begin{equation}
\begin{array}{l}
  {\ds m \frac{d^2 x}{d t^2}(t) + c \frac{d x}{d t}(t) + k x(t) = f(t) } \\
  \noalign{\bigskip}
  {\ds x(0) = x_0 \ , \ \dot{x}(0) = v_0 }
\end{array}
\label{eq1}
\end{equation}
where $m, c$ and $k$ respectively denote the mass, damping and stiffness
coefficients and $x(t)$ is the vertical displacement of the mass about the
equilibrium position.
We know that the analytic solution to (\ref{eq1}) with $f(t) = 0$ is
\begin{equation*}
  x(t) = e^{-ct/2m} [ A \cos(\nu t) + B \sin(\nu t)]
\end{equation*}
where
$$
  \nu = \frac{\sqrt{4km - c^2}}{2m} \ , \ A = x_0 \ , \ 
  B = \left(v_0 + \frac{c}{2m}x_0 \right)/\nu.
$$
To numerically approximate the solution to (\ref{eq1}), re-write as a
first order linear system
\begin{equation*}
  \frac{d \vec{y}}{d t}(t) = A \vec{y} + \vec{F}
\end{equation*}
where
$$
  A = \left[ \begin{array}{cc}  0 & 1 \\ 
                                \noalign{\smallskip} 
                                -K & -C \end{array} \right]
  \quad , \quad
  \vec{F}(t) = \left[ \begin{array}{c} 0 \\ \noalign{\smallskip}
                                       F(t) \end{array} \right] 
$$
and $\vec{y} = [x \, , \, \dot{x}]^T$, with the substitutions $C=c/m$,
$K=k/m$, and $F=f/m$. 

\end{list}

\clearpage

\begin{enumerate}

\item Mass-spring-dashpot

\begin{enumerate}

\item \label{numerical}Using reasonable values for $m,c,k$ solve the homogeneous
system in MATLAB using the routine of your choice (e.g., {\tt ode23, ode45,
ode15s}, or roll your own).  Plot the numerical solution $x^h(t)$ on the
same graph as the analytical solution for some meaningful time interval.
\underline{Comment} on any discrepancies. 

On your graph, you should \underline {label} the horizontal 
axis as {\em time, t}, the vertical axis as {\em x(t)}, and
place in the \underline {title} of the graph a text string
showing values of $m$, $c$, and $k$ ({\tt sprintf} may be helpful).


\item In this exercise, we will create {\em ``simulated''} data to be used
for estimating the unknown parameters.  Use the analytic solution to
generate discrete data $x_j=x(t_j)$ on a meaningful interval $[0,T]$ at $M$
equally-spaced time points  
\begin{equation*}
  t_j = \frac{jT}{M-1}.
\end{equation*}
Use, for example $M=100$.

\item Define a least squares objective function for the inverse problem to
determine parameters $\vec{q}=[C,K]$ given data $\vec{x}$ and appropriate
initial conditions (known).  Use a numerical solution method to compare to
data (although an analytic solution is possible for this model, it is not
available in general).

\item \label{initial}Using several ``nearby'' initial guesses, and several
``other'' initial guesses, apply a black box optimization routine to
minimize the least squares functional.  I recommend {\tt lsqnonlin}, but
note that it wants the residuals, not the sum of the squares of residuals.
Feel free to use any package you are familiar with.  Report in a table the
initial guess, initial cost, final estimate, final cost, and either the
computation time required to get the estimate or the number of function
evaluations required.  \underline{Comment} on any failures.

\item In practice, the data collected is corrupted by noise (for example,
errors in collecting data, instrumental errors, etc.). In the next part
of the exercise, we wish to test the sensitivity of the inverse least
squares method to errors in sampling the data. For this, we will add to
each simulated datum an error term as follows.  Artificially add noise to
your data with noise level (i.e., standard deviation) $b$ by
\begin{equation*}
d_j = x_j + \eta_j
\end{equation*}
where $\eta_j \sim N(0,b^2)$. In MATLAB this is {\tt d=x+b*randn(M,1);}.

\item \label{noise}Choose one initial estimate from \ref{initial}, and a
variety of increasing noise levels (starting at zero, e.g. 0, 0.01, 0.05,
0.1), and solve the inverse problem again, reporting results in a table as
before.  Additionally, make a plot of the minimum objective function value
as a function of noise level.  \underline{Describe} the sensitivity of the
inverse least squares method with respect to the noise level $b$.

\item Finally, compute the confidence interval for each of the estimates
from \ref{noise}.  Plot estimates for $C$ with confidence intervals as a
function of noise level (see {\tt errorbar}).  Repeat for $K$.
\underline{Comment} on what you observe.

\end{enumerate}

\item Vibrating Beam

Download BeamData.mat from the course website and run {\tt load BeamData}
in MATLAB.  This will 
create variables $Td$ and $Ad$ which contain {\em accelerometer} data
collected from a vibrating beam which was excited at its first fundamental
frequency (approximately 11 Hz). The data is truncated to a time point at
which acceleration is zero (use this to determine the initial position!).
See Figure \ref{fig1} for an idealized picture.

\begin{figure}[!hh]
\bigskip
\begin{center}
  \includegraphics[width=400pt]{fig1}
\end{center}

\vspace{-.25truein}
\caption{Beam Excitation: Shown is idealized accelerometer data from a
  vibrating beam.  Sinusoidal input is stopped at $t=10$, data collection
  begins at $t=15$ where $\ddot{x}=0$.
}
\label{fig1}
\end{figure}


We will use the mass-spring-dashpot model in \eqref{eq1} to attempt to
represent the dynamics of the beam data.  

\begin{enumerate}
\item Plot $Ad$ versus $Td$ to visualize the data.

\item Formulate the least squares parameter identification
problem for $\vec{q}=[C,K,v_0]$ using {\em acceleration data}.

\item \label{alt} Solve this inverse problem using the black box
optimization routine of your choice, for example, {\tt lsqnonlin} or {\tt
fminsearch} ({\bf Bonus: use both and compare results!}).  Use a numerical
solution method to compare to data (see \ref{numerical}).  Do not take the
derivative of your solution to the ODE to get the acceleration, but rather,
plug $x$ and $\dot{x}$ into the ODE.  For an initial
guess, try any three numbers of your choice (note $C\geq0$ and $K\geq0$).
Report your initial guess, initial cost, final estimate and final cost.
Plot the acceleration corresponding to your final estimate on the same plot
as the data.  Just for fun, try other random initial guesses.  Does there
seem to be a local minimum to the objective function that has a very large
basin of attraction?  (Feel free to make a surface plot, or better yet, a
volumetric plot (see {\tt slice}) if you're interested.)

\item \label{alt2} Solve the inverse problem again with an initial guess in
the vicinity of $[0.002, 4850, 1.3]$.  Report your initial guess, initial
cost, final estimate and final cost. Plot the acceleration corresponding to
your final estimate on the same plot as the data. \underline{Describe} any
discrepancy in your final solution versus the data.  \underline{Explain}
why there may be this discrepancy. 


\item Repeat \ref{alt2} with the assumption that $C=0$ in the model.

\item Use the $\chi^2(1)$ test to test for the significance of your improved fit
to the data by allowing nontrivial damping $C\neq0$ in the model.  State in
terms of a Null and Alternate Hypothesis.  Does your answer make sense from
looking at the data? \underline{Explain} why or why not.


\end{enumerate}
\end{enumerate}


\end{document}

