added root finding and BVP

Lab.tex

@@ -83,6 +83,13 @@
 \section{Ordinary differential equations: Initial value problems}
 \section{Ordinary differential equations: Initial value problems}
 \input{ode.tex}
 \input{ode.tex}
 
 
+\section{Root finding}
+\input{root.tex}
+
+\section{Ordinary differential equations: Two-point boundary value problems}
+\input{ode_bvp.tex}
+
+
 %\section{Lab 2}
 %\section{Lab 2}
 %\input{rest.tex}
 %\input{rest.tex}
 
 

ode_bvp.tex

+
+
+Solve the TOV equations described in Sec.~\ref{sec:TOV} as a two-point boundary value problem using the Shooting method. Use the Newton-Raphson method for root finding. The boundary conditions are given in Eq.~(\ref{eq:tov_boundary}).
+

rest.tex

+
+\subsection{Ordinary differential equations: Calculation of gravitational waves from inspiralling compact binaries}
+
+\subsubsection*{Problems:}
+The time evolution of the orbital phase $\varphi(t)$ of a binary of black holes evolving under the gravitational radiation reaction can be computed, in the post-Newtonian approximation, by solving the following coupled ODEs: 
+\begin{eqnarray}
+\frac{dv}{dt} = -\frac{\mathcal{F}(v)}{dE(v)/dv}, ~~~~~ \frac{d\varphi}{dt} = \frac{v^3}{m},
+\label{eq:phasing_formula}
+\end{eqnarray}
+where $E(v)$ is the binding energy of the orbit, $\mathcal{F}(v)$ is the energy flux of radiated gravitational waves, $m:= m_1 + m_2$ is the total mass of the binary, $v = (m \omega)^{1/3}$, $\omega$ being the orbital frequency. (Here we use geometric units, in which $G = c = 1$. This means that in all the expressions $m$ has to be replaced by $Gm/c^3$.)  
+
+The binding energy and gravitational-wave flux are given as post-Newtonian expansions in terms of the small parameter $v$ 
+\begin{eqnarray}
+E(v) = -\frac{1}{2} \mu v^2 \left[1 + \mathcal{O}(v^2) \right], ~~~ \mathcal{F}(v) = \frac{32}{5} \left(\frac{\mu}{m}\right)^2 \, v^{10} \left[1 + \mathcal{O}(v^2) \right], 
+\end{eqnarray}
+where $\mu :=  m_1 m_2/m$ is the reduced mass of the system. 
+
+\begin{enumerate}
+\item Compute $v$ as a function of $t$ by solving the first equation in Eq.~(\ref{eq:phasing_formula}) using Scipy's \href{http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.odeint.html}{\texttt{odeint}} routine. Assume the following parameters $m_1 = m_2 = 5 M_\odot$, $v_0 = 0.3$, $\varphi_0 = 0$. Plot $v(t)$. 
+\item Solve the coupled system in Eq.~(\ref{eq:phasing_formula}) to compute $v$ and $\varphi$. Compute and plot the two gravitational-wave polarizations: 
+\begin{equation}
+\label{eq:GWpolzn}
+h_+(t)  = 4 \frac{\mu}{m} v^2 \, \cos \varphi(t), ~~~ h_\times(t)  = 4 \frac{\mu}{m} v^2 \, \sin \varphi(t). 
+\end{equation}
+\item Repeat the same calculation using the adaptive Runge-Kutta method by Dormand \& Prince using Scipy's \href{http://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.ode.html}{\texttt{integrate.ode}} module (choose method = ``dopri5''). This method is very similar to the Runge-Kutta-Fehlberg method that we learned in the class. 
+\end{enumerate}
+
+%
+\subsection{Ordinary differential equations: Structure of a relativistic, spherically symmetric star}
+\label{sec:TOV}
+%
+The interior structure of a relativistic, spherically symmetric star is described by a metric that has the line element 
+\begin{equation}
+ds^2 = - e^{2\Phi(r)} \, c^2 dt^2 + \left(1 - \frac{2Gm(r)}{rc^2} \right)^{-1} \, dr^2 + r^2 \, d\Omega^2,
+\end{equation}
+where $m(r)$ is called the \emph{mass function} (which encapsulates the {gravitational mass} inside the radius $r$), $e^{2\Phi(r)}$ is the \emph{lapse function} (which relates the proper time with the coordinate time). Outside the ``surface'' of the star (in vacuum), the spacetime is described by the Schwarzschild metric; the lapse function becomes 
+\begin{equation}
+\Phi(r) = \frac{1}{2} \ln \left(1 - \frac{2Gm_\star}{rc^2} \right)
+\label{eq:lapse_schwarz}
+\end{equation}
+where $m_\star$ is the total (gravitational) mass of the star. The structure can be computed by solving the following set of ordinary differential equations, derived by Tolman, Oppenheimer and Volkoff (TOV). 
+\begin{eqnarray}
+\frac{dm(r)}{dr} & = &{4 \pi} r^2 \rho(r), \label{eq:tov1} \\
+\frac{dP(r)}{dr} & = & -\frac{G (\rho + P(r)/c^2)}{r^2} \left[m(r) + \frac{4\pi r^3 P(r)}{c^2} \right] \left[1- \frac{2G m(r)}{c^2 r} \right]^{-1}, \label{eq:tov2} \\
+\frac{d\Phi(r)}{dr} & = & \frac{G m(r) + 4 \pi G r^3 P(r) / c^2}{c^2 r  [r - 2Gm(r)/c^2]}. \label{eq:tov3}
+\end{eqnarray}
+The TOV equations have to be supplemented by an \emph{equation of state} $P = P(\rho)$, that relates the pressure $P$ to the energy density $\rho$. We assume the equating of state to be of polytropic form
+\begin{equation}
+P(r) = K \, \rho(r) ^{\gamma}.
+\end{equation}
+We also need to specify initial conditions for the variables $m, P, \Phi$. The following conditions can be used 
+\begin{equation}
+\label{eq:tov_boundary}
+m ( r = 0) = 0, ~~~~ P(r = 0) = P_c = P(\rho_c), ~~~~ \Phi(r_\star) = \frac{1}{2} \ln \left(1 - \frac{2Gm_\star}{r_\star c^2} \right) 
+\end{equation}
+
+\subsubsection*{Problems:}
+\begin{enumerate}
+\item Compute the structure, i.e., $m(r)$ and $P(r)$,  of a neutron star with central density $\rho_c = 5 \times 10^{17} ~ \mathrm{kg/m^{3}}$ by solving Eqs.~(\ref{eq:tov1}) and (\ref{eq:tov2}). Assume a polytropic equation of state with $\gamma = 5/3$ and $K = 5380.3$ (SI units). What is the mass $m_\star$ and radius $r_\star$ of the neutron star?  (Useful tip: You will need to start the integration at $r = \Delta r$, where $\Delta r$ is a small number. You can assume $m(r = \Delta r) := 4/3 \, \pi \rho_c  (\Delta r)^3$).
+\item Compute the lapse function $e^{2\Phi(r)}$ by solving Eq.(\ref{eq:tov3}) starting from $r = r_\star$ to $r = 0$. On top of that, plot the lapse function for a Schwarzschild black hole (see Eq.\ref{eq:lapse_schwarz}) from $r = r_s$ to $r = 2 r_\star$, where $r_s \equiv 2 G m_\star/c^2$ is the Schwarzschild radius of the star. This exterior solution should match the interior solution at $r = r_\star$.  
+\end{enumerate}
+
+\subsection{Non-linear ordinary differential equations showing chaotic behavior: Lorenz equations}
+The Lorenz equations were originally developed as a simplified mathematical model for atmospheric convection by Edward Lorenz. This was the first set of equations where deterministic chaos was observed. The equations are given by 
+\begin{equation}
+\frac{dx(t)}{dt} = \sigma [y(t)-x(t)], ~~~~ \frac{dy(t)}{dt} = x(t) [\rho - z(t)] - y(t),~~~~ \frac{dz(t)}{dt} = x(t)y(t) - \beta z(t),
+\end{equation}
+where $\rho, \sigma$ and $\beta$ are parameters of the system. 
+\subsubsection{Problems:}
+\begin{enumerate}
+\item Solve the Lorenz system for $\rho = 28, \sigma = 10$ and $\beta = 8/3$ with the following initial conditions $x(t = 0) = y(t = 0) = z(t=0) = 1$. Plot $x(t)$, $y(t)$ and $z(t)$ for $t = 0 ... 100$. Is the solution deterministic or stochastic? 
+\item Repeat the calculation with same parameters except for a tiny change in the initial condition for $x$: i.e., $x(t = 0) = 1 + 10^{-9}$. Plot  $x(t = 0) = y(t = 0) = z(t=0) = 1$ on top of the earlier estimate. Explain the result. 
+\item Make a 3D plot of $x, y, z$. You should see the famous butterfly shaped structure now! 
+\item Optional exercise: Make an animation of the above~\footnote{You can either use the matplotlib \href{http://matplotlib.org/api/animation_api.html}{animation} package or convert a number of PNG files to a gif animation using \href{http://www.imagemagick.org/Usage/anim_basics/}{ImageMagick}.}. 
+\end{enumerate}
+
+\subsection{Stochastic ordinary differential equations: Langevin equation}
+
+The random motion of a particle in a fluid due to collisions with the molecules of the fluid, called the Brownian motion, is described by the Langevin equation: 
+\begin{equation}
+m \frac{d^2 \bx}{dt^2} = - \lambda \frac{d \bx}{dt} + \boldeta(t), 
+\end{equation}
+where $m$ is the mass of the particle, $\bx$ its position vector, $\lambda$ a damping coefficient, and $\boldeta(t)$ (called the \emph{noise term}) describes the stochastic processes affecting the system. 
+
+\subsubsection{Problems:}
+\begin{enumerate}
+\item Using the Euler-Maruyama method, compute the $1-d$ Brownian motion trajectories generated by the Langevin equation
+${d x}/{dt} = \eta(t)$,  
+where $\eta(t)$ is Gaussian noise with zero mean and unit variance. Plot $x(t)$ for $ t \in [0, 100]$ assuming $x(t=0) = 0$. Is the solution deterministic or stochastic? 
+\item Generalize the code so that it can deal with arbitrary number of particles. Plot $x(t)$ for 1000 particles on a single plot. Compute the average displacement $\bar{x}(t)$ of the particles (from $x(t=0)$) as a function of $t$ and plot it against $t$. What is the relation between  $\bar{x}(t)$ and $t$? 
+\item Using the \href{http://matplotlib.org/api/pyplot_api.html#matplotlib.pyplot.hist}{\texttt{hist}} function, plot the probability distribution $P(x)$ of $x(t)$ at $t = 10, 50, 100$.
+\end{enumerate}
+
+\subsection{Two-point boundary value problems: The shooting method}
+
+Solve the TOV equations described in Sec.~\ref{sec:TOV} as a two-point boundary value problem using the Shooting method. Use the Newton-Raphson method for root finding. The boundary conditions are given in Eq.~(\ref{eq:tov_boundary}).
+
+\section{Partial differential equations}
+
+\subsection{First order hyperbolic PDEs: The advection equation}
+
+A few explicit first order numerical methods for solving the advection equation 
+$\frac{\partial u}{\partial t} = -v \frac{\partial u}{\partial x}$
+are given below: 
+\begin{eqnarray}
+\mathrm{Forward~time~centered~space~(FTCS):} & u^{n+1}_j = \frac{-v \Delta t}{2\Delta x} \, \Big[{u^n_{j+1} - u^n_{j-1}}\Big] + u^n_j \\
+\mathrm{Upwind~method~(for~} v > 0): & u^{n+1}_j = \frac{-v \Delta t}{\Delta x} \, \Big[{u^n_{j} - u^n_{j-1}}\Big] + u^n_j \\
+\mathrm{Downwind~method~(for~} v < 0): & u^{n+1}_j = \frac{-v \Delta t}{\Delta x} \, \Big[{u^n_{j+1} - u^n_{j}}\Big] + u^n_j~.
+\end{eqnarray}
+Above, $u^{n+1}_j = u(t_n, x_j)$, where $x_j = x_0 + j \Delta x$ and $t_n = t_0 + n \Delta t$.
+\subsubsection{Problems:}
+\begin{enumerate}
+\item Solve the advection equation with $v = 1$ and initial condition $u(t=0,x) := \exp[-(x-x_0)^2]$ (Gaussian pulse) where $x \in [0, 25]$ and $x_0 = 5$ using the FTCS scheme, using the following condition for the left boundary $u(x=0,t) = 0$. Is the evolution stable?  
+\item Repeat the calculation with upwind methods using Courant factor $\lambda := |v\Delta t|/\Delta x = 1/2, 1, 2$. Which of the evolutions are stable? 
+\item Repeat the calculation using the periodic boundary condition $u(x=0, t) = u(x=L, t)$ where $x=L$ corresponds to the right boundary. 
+\item Plot the order of convergence of $u(x,t)$ as a function of $x$ at $t = 10$ (use $\lambda = 1)$.   
+\end{enumerate}
+
+\subsection{Second order hyperbolic PDEs: The wave equation}
+
+A simple explicit method for solving the wave equation $\frac{\partial^2 u}{\partial t^2} - v^2 \frac{\partial^2 u}{\partial x^2} = 0$ is given below: 
+\begin{eqnarray}
+\mathrm{For~the~first~time~step:} & u^{n+1}_j = \frac{\lambda^2}{2} \left( u^n_{j+1} + u^n_{j-1} \right) + \left(1 - \lambda^2 \right) u^n_j + \Delta t \frac{\partial u}{\partial t} \big |^n_j\\ 
+\mathrm{All~other~time~steps:} & u^{n+1}_j = \lambda^2 \left( u^n_{j+1} + u^n_{j-1} \right) + 2 \left(1 - \lambda^2 \right) u^n_j - u^{n-1}_j\, , 
+\end{eqnarray}
+where $\lambda = v\Delta t/\Delta x$ is the Courant factor.
+
+\subsubsection{Problems:}
+\begin{enumerate}
+\item Solve the wave equation with $v = 1$ and initial conditions $u(t=0,x) := \sin \pi x$ and $\frac{\partial u}{\partial t}(t=0,x) := 0$, and Dirichlet type boundary conditions $u(t,x=0) = u(t,x=L) = 0 $ over the domain $x \in (0,L)$ and $t \in (0,T)$. Choose $L = 4$ and $T = 20$. Be sure to choose $\lambda \leq 1$. This shows oscillations on a string whose end points are fixed.  
+\end{enumerate}
+
+\subsection{Second order parabolic PDEs: The diffusion equation}
+
+The FTSC method provides a simple explicit scheme for solving the diffusion equation $\frac{\partial u}{\partial t} - k \frac{\partial^2 u}{\partial x^2} = 0$:
+\begin{eqnarray}
+u_j^{n+1} = \gamma \left(u^n_{j+1} + u^n_{j-1} \right) + (1 - 2 \gamma) u^n_j
+\end{eqnarray}
+where $\gamma = \frac{k\Delta t}{\Delta x^2}$. The FTSC scheme can produce stable evolutions for $\gamma \leq 0.5$. 
+
+\subsubsection{Problems}
+\begin{enumerate}
+\item Solve the diffusion equation with $k = 0.5$ over the domain $x \in [0, 2],~t \in [0, 10]$ with the following initial and boundary conditions: 
+\begin{eqnarray}
+\label{eq:heat_ini_cond}
+u (t=0, x) & = & 
+\begin{cases}
+100 x,~\mathrm{for}~ 0 \leq x \leq 1 \\
+200-100 x,~\mathrm{for}~ 1 < x \leq 2 \\
+\end{cases}\\
+u(t, x=0) & = & 50 \\
+u(t, x=2) & = & 10 
+\end{eqnarray}
+This shows the evolution of the temperature profile (initial profile given by Eq.~\ref{eq:heat_ini_cond}) of a one-dimensional object connected to hot baths at both ends, maintained at temperatures 50 K and 10 K. 
+\end{enumerate}
+
+\subsection{Second order elliptic PDEs: The Laplace equation}
+
+The Laplace equation $\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0$ in two dimensions can be solved using the successive relaxation method. This explicit method involves successively iterating the following equation 
+\begin{equation}
+u_{i,j} = \frac{u_{i-1,j} + u_{i+1,j} + u_{i,j-1} + u_{i,j+1}}{4}.
+\end{equation}
+Above, $u_{i,j} = u(x_i, y_j)$, where $x_i = x_0 + i h$ and $y_j = y_0 + j h$. The iteration has to be repeated until the difference of $u_{i,j}$ between two successive iterations is less than the required tolerance. That is, $\epsilon := \sum_{i,j} \, |u_{i,j}^{n} - u_{i,j}^{n-1}| \leq \mathrm{tol}$.
+
+If a von Neumann type boundary condition $\frac{\partial u}{\partial x} = A(y)$ is specified at $x = 0$, then the $i=0$ grid points can be evaluated using the following equation  
+\begin{equation}
+u_{0,j} = \frac{2 u_{1,j} - 2 h A_j + u_{0,j-1} + u_{0,j+1}}{4}.
+\end{equation}
+
+\subsubsection{Problems}
+
+\begin{enumerate}
+\item Solve the Laplace equation over a $100 \times 100$ grid  with the following Dirichlet type boundary conditions:
+\begin{equation}
+u(x = 0, y) = u(x, y=0) = u(x, y=L) = 0, ~~ u(x = L, y) = 100.  
+\label{eq:laplace_bc}
+\end{equation}
+Plot $u(x,y)$ once the system has reached equilibrium state $\epsilon \lesssim 10^{-6}$). This shows the steady-state temperature distribution on a two-dimensional plate which is kept at a fixed temperature of zero degree at all boundaries, except at the right boundary (which is kept at 100 K). 
+\item Plot the error $\epsilon$ as a function of iteration. One should see the error monotonically decreasing. 
+\item Repeat the calculation, except that the boundary condition at $x = 0$ in Eq.\,(\ref{eq:laplace_bc}) is replaced by $\frac{\partial u}{\partial x} = 10$.
+\end{enumerate}
+
+\section{Fourier and spectral methods}
+\subsection{Fast Fourier transform (FFT)}
+\subsubsection{Problems:}
+\begin{enumerate}
+\item Hercules X-1 is a high-mass X-ray binary (HMXB) system having a magnetized spinning neutron star which happens to be an X-ray pulsar. Here~\cite{rxte-data} you are given the data of an RXTE~\cite{rxte} observation after cleaning and pre-processing. The file contains two columns, (i) time (in seconds) and (ii) count-rate of X-ray photons (i.e., number of photons detected per unit time). Plot the count-rate as function of time. This is called lightcurve in X-ray astronomy. Compute the FFT of the given time-series and plot its absolute value against the frequency. Are you able to find any periodic signal(s)? The lowest frequency signal corresponds to the spin frequency of the neutron star.  What is the spin period of this pulsar? Can you identify other periodic signals and how they are related to the lowest frequency signal?~\footnote{This problem and data are graciously provided by Dr. Arunava Mukherjee (IUCAA).}
+\end{enumerate}
+
+\subsection{Power spectrum estimation using FFT}
+\subsubsection{Problems:}
+\begin{enumerate}
+\item A sample data set from the LIGO gravitational-wave observatory can be downloaded from here~\cite{ligo-data}. This contains 256 seconds of LIGO data from 2005, sampled at a rate of 4096 Hz. Compute the power spectral density of the data using Welch's modified periodogram method. 
+\end{enumerate}
+
+\subsection{Time-frequency signal detection methods}
+\subsubsection{Problems}
+\begin{enumerate}
+\item 4U 1636--536 is a low-mass X-ray binary (LMXB) system consisting of a rapidly spinning weakly magnetized neutron star which does not exhibit any coherent X-ray pulsation like Hercules X-1. However, X-ray radiation from this source shows signals of high frequency (500--1000 Hz) quasi-periodic oscillations (QPOs). Here~\cite{rxte-data2} you are given another data of RXTE observation. This contains the arrival times of photons. Identify the frequency span of this QPO using a time-frequency spectrogram or a PSD. {Note: Firstly, you need to sample the data with a fixed sampling rate, say 4096 Hz}.  
+\end{enumerate}
+
+\subsection{Computing correlations using FFT: Matched filtering}
+In the case a known signal $h(t)$ buried in stationary Gaussian, white noise, the optimal technique for 
+signal extraction is the \emph{matched filtering}, which involves cross-correlating the data 
+with a \emph{template} of the signal. The correlation function between two time series $x(t)$ and $\hat h(t)$ (both assumed to be real-valued) for a time shift 
+$\tau$ is defined as:
+\begin{equation}
+\label{eq:corrFn}
+R(\tau) = \int_{-\infty}^{\infty}x(t)\,\hat{h}(t+\tau)\,dt. 
+\end{equation}
+Above, $\hat{h(t)} := h(t) / ||h||$, where the norm $||h||$
+of the template is defined by $$||h||^2 = \int_{0}^{t_c} \, |h(t)|^2/\sigma^2 \, dt,$$ where 
+$\sigma^2$ is the variance of the noise. 
+The optimal signal-to-noise ratio (SNR) is obtained when the template exactly matches with the signal. i.e., $\mathrm{SNR}_\mathrm{opt} = ||h||$. 
+If the SNR is greater than a predetermined threshold (which corresponds to an acceptably small false alarm 
+probability), a detection can be claimed. 
+
+\subsubsection{Problems}
+\begin{enumerate}
+\item You are given a time-series data set $d(t)$ here~\cite{mock-gw-data}. This contains a simulated gravitational-wave signal from a black hole binary buried in zero-mean, Gaussian white noise $n(t)$ with standard deviation $\sigma = 10^{-21}$. i.e.,
+\begin{equation}
+d(t) = n(t) + m \, h_+(t) / D_L,
+\end{equation}
+where $h_+(t)$ is given by Eq.~(\ref{eq:GWpolzn}), $m = m_1 + m_2$ is the total mass of the binary, and $D_L$ is the (unknown) luminosity distance to the binary.  Detect the location of the signal in the data by maximizing the correlation of the waveform templates computed in Eq.~(\ref{eq:GWpolzn}) with the data.: 
+\begin{equation}
+R_\mathrm{max} = \mathrm{max}~_{m,\mu,\tau}~ R(\tau),
+\end{equation}
+where $R(\tau)$ is given by Eq.~(\ref{eq:corrFn}). We have the prior information that the parameters of the signal are in the following range: $5 < m/M_\odot < 15$, $0.2 < \mu/m < 0.25$. Estimate the parameters $m,\mu$ of the signal (parameters that maximize the correlation). 
+
+\end{enumerate}
+

root.tex

+
+The Kepler's equation is a  transcendental equation describing planetary motion:
+\begin{equation}
+M - E + e \sin E = 0. 
+\label{eq:keplers_eqn}
+\end{equation}
+This gives the the relation between the mean anomaly $M$ (a parametrisation of time) and the eccentric anomaly $E$ (parametrisation of the polar angle of the planet, given the orbital eccentricity $e$. The mean anomaly can be expressed in terms of the time coordinate $t$ as 
+\begin{equation}
+M = \frac{2\pi}{T} (t - \tau), 
+\label{eq:mean_anomaly}
+\end{equation}
+where $T$ is the period of the orbit and $\tau$ is the time when the planet reaches the periapsis. Combining this with the Kepler's equation, we can find $E$ as a function of $t$. The coordinates of the planet are then given by 
+\begin{equation}
+x = a  \, (\cos E - e), ~~~~~ y = b \, \sin E, 
+\end{equation}
+where $a$ and $b$ are the semi-major and semi-minor axes, respectively. 
+
+\subsubsection{Problems:}
+\begin{enumerate}
+\item Solve the Kepler's equation using the Newton-Raphson's method. Plot $E$ as a function of $M$ for three different values of $e = 0, 0.1, 0.5$. Take $M$ to be in the range $[0, 2 \pi)$. You can use scipy's \href{https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root_scalar.html#scipy.optimize.root_scalar}{\texttt{optimize.root\_scalar}} function to find the roots. 
+\item Plot the orbits of solar system planets. The necessary data can be downloaded from the \href{https://nssdc.gsfc.nasa.gov/planetary/factsheet/}{NASA Planetary Fact Sheet}.  
+\item Optional exercise: Make an animation of the above. 
+
+ 
+\end{enumerate}
+