Skip to content

Parameswaran Ajith / NM2025

Go to a project
Commit 9cbff345 by Parameswaran Ajith
Options

added the problems from last year

parents
Expand all
Inline Side-by-side
Showing with 619 additions and 0 deletions
Lab.bbl 0 → 100644
\begin{thebibliography}{7}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
\def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
\def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname citenamefont\endcsname\relax
\def\citenamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
\def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}
\bibitem[{nrd()}]{nrdata}
\urlprefix\url{http://home.icts.res.in/~ajith/Downloads/nr_data.gz}.
\bibitem[{SXS()}]{SXScatalog}
\urlprefix\url{http://www.black-holes.org/waveforms/}.
\bibitem[{rxt({\natexlab{a}})}]{rxte-data}
\urlprefix\url{http://home.icts.res.in/~ajith/Downloads/extracted_lightcurve_HerX-1.dat.gz}.
\bibitem[{rxt({\natexlab{b}})}]{rxte}
\bibinfo{note}{RXTE is an X-ray timing satellite},
\urlprefix\url{https://heasarc.gsfc.nasa.gov/docs/xte/rxte.html}.
\bibitem[{lig()}]{ligo-data}
\bibinfo{note}{Download the file L1-STRAIN\_4096Hz-815045078-256.txt.gz from},
\urlprefix\url{http://www.ligo.org/science/GRB051103/index.php}.
\bibitem[{snd()}]{sndata}
\urlprefix\url{https://home.icts.res.in/~ajith/Downloads/SCPUnion2.1_mu_vs_z.txt}.
\bibitem[{SNC()}]{SNCosmology}
\urlprefix\url{https://supernova.lbl.gov/}.
\end{thebibliography}
Lab.bib 0 → 100755
This diff is collapsed. Click to expand it.
Lab.tex 0 → 100644
\documentclass[prd,preprintnumbers,eqsecnum,floatfix,a4paper,nofootinbib]{revtex4}
\usepackage{color}
\usepackage{url}
\usepackage{calc}
\usepackage{amsmath,amssymb,graphicx}
\usepackage{amssymb,amsmath}
\usepackage{tensor}
\usepackage{bm}
\usepackage{times}
\usepackage[varg]{txfonts}
\usepackage[colorlinks, pdfborder={0 0 0}]{hyperref}
\definecolor{LinkColor}{rgb}{0.75, 0, 0}
\definecolor{CiteColor}{rgb}{0, 0.5, 0.5}
\definecolor{UrlColor}{rgb}{0, 0, 0.75}
\hypersetup{linkcolor=LinkColor}
\hypersetup{citecolor=CiteColor}
\hypersetup{urlcolor=UrlColor}
\maxdeadcycles=1000
\allowdisplaybreaks
\newcommand{\comment}[1]{\textcolor{red}{\textit{#1}}}
\newcommand{\checkthis}{\textcolor{magenta}{(CHECKTHIS)}}
\newcommand{\bx}{\mathbf{x}}
\newcommand{\boldeta}{\bm{\eta}}
\begin{document}
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\ber}{\begin{eqnarray}}
\newcommand{\eer}{\end{eqnarray}}
\def\bea{\begin{eqnarray}}
\def\eea{\end{eqnarray}}
\newcommand{\etal}{\emph{et al.}}
\newcommand{\Sl}{S_\ell}
\newcommand{\Sigmal}{\Sigma_\ell}
\newcommand{\Flux}{\mathcal{F}}
\newcommand{\LNh}{\hat{\mathbf{L}}_N}
\newcommand{\LN}{\mathbf{L}_N}
\newcommand{\bS}{\mathbf{S}}
\newcommand{\bJ}{\mathbf{J}}
\newcommand{\e}{\mathrm{e}}
\newcommand{\rmi}{\mathrm{i}}
\newcommand{\flow}{f_\mathrm{low}}
\newcommand{\fcut}{f_\mathrm{cut}}
\newcommand{\bchi}{\bm{\chi}}
\newcommand{\blambda}{\bm{\lambda}}
\newcommand{\bLambda}{\bm{\Lambda}}
\newcommand{\bchia}{\bm{\chi}_a}
\newcommand{\bchis}{\bm{\chi}_s}
\newcommand{\chis}{\chi_s}
\newcommand{\chia}{\chi_a}
\newcommand{\chiadL}{\bchia \cdot \LNh}
\newcommand{\chisdL}{\bchis \cdot \LNh}
\newcommand{\chisSqr}{\bchis^2}
\newcommand{\chiaSqr}{\bchia^2}
\newcommand{\chisDchia}{\bchis \cdot \bchia}
\newcommand{\cA}{\mathcal{A}}
\newcommand{\cB}{\mathcal{B}}
\newcommand{\cC}{\mathcal{C}}
\newcommand{\cP}{\mathcal{P}}
\title{ICTS Graduate Course: Numerical Methods (PHY410.5)}
\author{Parameswaran~Ajith}\email{ajith@icts.res.in}
\author{Prayush~Kumar}\email{prayush@icts.res.in}
\affiliation{International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bangalore 560089, India.}
\author{Koustav Narayan Maity (tutor)}\email{koustav.narayan@icts.res.in}
\author{Vinay Kumar (tutor)}\email{vinay.kumar@icts.res.in}
\affiliation{International Centre for Theoretical Sciences, Tata Institute of Fundamental Research, Bangalore 560089, India.}
\bigskip
\date{\today}
\maketitle
\section{Numerical differentiation}
\input{diff.tex}
\section{Numerical integration}
\label{sec:integr}
\input{integration.tex}
\section{Ordinary differential equations: Initial value problems}
\input{ode.tex}
\section{Root finding}
\input{root.tex}
\section{Ordinary differential equations: Two-point boundary value problems}
\input{ode_bvp.tex}
\section{Fourier and spectral methods}
\input{fourier.tex}
\section{Curve fitting}
\label{sec:curve_fitting}
\input{curve.tex}
\section{Spectral methods for ODEs and PDEs}
\input{PDE1.tex}
\section{Statistical inference}
\label{sec:statinf}
\input{statinf.tex}
\section{Monte-Carlo methods}
\input{mc.tex}
%\section{Lab 2}
%\input{rest.tex}
\bibliography{Lab}
\end{document}
PDE1.tex 0 → 100644
\subsection{Elliptic equations}
We consider the example discussed in class,
\begin{equation}
\dfrac{d^2 u}{dx^2} - x^6(3+x^2)u = 0,
\end{equation}
on the interval $ [-1,1] $ with boundary conditions given by $ u(\pm 1) = 1 $.
\subsubsection{Problems}
\begin{enumerate}
\item Discretize this equation using Chebyshev polynomials. This will result in a set of equations for N unknowns, where two of the equations correspond to boundary conditions. Use direct matrix-inversion to find the solution.
\item Provide plots of the solution.
\item Solve for various different resolutions $N$, and provide plots that demonstrate how accurate your solution is, and that your numerical solution converges with $N$ (i.e. error goes to 0 as N increases).
\item How many modes are needed for an accuracy of $10^{-5}$? How many for $ 10^{-12} $?
\end{enumerate}
\subsection{Hyperbolic equations [Optional]}
Consider the wave equation in 1+1 dimensions, as discussed in class, over the interval $ [0,2\pi] $,
\begin{equation}
\Box g = -\partial_t ^2 g + \partial_x ^2 g = 0.
\end{equation}
This can be written as a set of first order differential equations,
\begin{align*}
\dot{g} &= -\Pi, \\
\dot{\Pi} &= -\partial_x \Phi, \\
\dot{\Phi} &= -\partial_x \Pi.
\end{align*}
\subsubsection{Problems}
Solve the above equations using Fourier polynomial basis. Choose an initial condition which looks like a localized pulse. Enforce the scalar field to be 0 at both boundaries, i.e. set
\begin{equation}
g(0) = \Pi(0) = \Phi(0) = g(2\pi) = \Pi(2\pi) = \Phi(2\pi) = 0.
\end{equation}
\begin{enumerate}
\item What do you observe if you evolve for times longer than $2\pi$?
\item Increase the domain size without changing the initial pulse size. Comment on the behaviour of the scalar in this case.
\end{enumerate}
\ No newline at end of file
curve.tex 0 → 100644
\subsection{Linear regression}
We have seen in the class that, given a theoretical model $Y(x, \textbf{a}) = \sum_{j=1}^M a_j Y_j(x)$ and a data set $\{x_i, y_i\}$ with associated errors $\sigma_i$, we can estimate the parameters $\textbf{a}$ of the model by minimising the cost function (chi-square). This amounts to solving the following Normal equations of the least square problem
\begin{equation}
(\textsf{A}^T \textsf{A}) ~ \textbf{a} = \textsf{A}^T \textbf{b},
\end{equation}
where the elements of the matrix $\textsf{A}$ and vector $\textbf{b}$ are given by
\begin{equation}
A_{ij} = \frac{Y_j(x_i)}{\sigma_i}, ~~~ b_i = \frac{y_i}{\sigma_i}.
\end{equation}
\subsubsection{Problems}
\begin{enumerate}
\item Show that the covariance matrix of the estimated parameters $\textbf{a}$ are given by $\textsf{C} = (\textsf{A}^T \textsf{A})^{-1}$
\item The Hubble's law provides a relation between the luminosity distance $d_L$ and cosmological redshift $z$ that is valid in the local universe ($z \lesssim 0.1$).
\begin{equation}
d_L = \frac{c}{H_0} z,
\label{eq:Hubble_law}
\end{equation}
where $c$ is the speed of light and $H_0$ is the Hubble constant. Here~\cite{sndata} you are given a dataset from the Supernova Cosmology project~\cite{SNCosmology}. This contains the redshift $z$, the distance modulus $\mu$, and the error on the distance modulus $\delta \mu$ measured from several Type 1a supernovae. The distance modulus is related to the luminosity distance (in parsecs) by $\mu = 5 \left(\log_{10} d_L - 1 \right)$. Use the linear regression method to fit the distance and redshift data to estimate the Hubble constant and the associated error. Note that the data contain redshifts above the range where the simple Hubble law is valid. Please select the samples with $z \leq 0.1$.
\end{enumerate}
\subsection{Non-linear least square fitting}
In the flat $\Lambda$CDM model of cosmology, the general relation between he luminosity distance $d_L$ and cosmological redshift $z$ is given by
\begin{equation}
d_L = (1+z) \frac{c}{H_0} \int_0^z \frac{dz'}{\sqrt{\Omega_M(1+z')^3 + \Omega_\Lambda)}},
\label{eq:lcdm}
\end{equation}
where $\Omega_M$ is the energy density of matter and $\Omega_\Lambda$ is the energy density of the cosmological constant. Since we are assuming spatially flat universe: $\Omega_M +\Omega_\Lambda = 1$. Thus the only free parameters of the model are $H_0$ and $\Omega_M$.
\subsubsection{Problems}
\begin{enumerate}
\item Use the supernova data from the entire redshift to estimate $H_0$ and $\Omega_M$ and the corresponding covariance. This is a non-linear least square problem. You can use SciPy's \href{https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html}{\texttt{curve\_fit}} function.
\end{enumerate}
diff.tex 0 → 100644
\subsection{Finite differencing, convergence, error estimates}
We derived the following finite-differencing approximants for the derivative of a function $f(x)$:
\begin{eqnarray}
\mathit{Forward~differencing:} & f'(x) \simeq \frac{f(x+h)-f(x)}{h} + \mathcal{O}\,(h) \\
\mathit{Backward~differencing:} & f'(x) \simeq \frac{f(x)-f(x-h)}{h} + \mathcal{O}\,(h) \\
\mathit{Central~differencing:} & f'(x) \simeq \frac{f(x+h)-f(x-h)}{2h} + \mathcal{O}\,(h^2)
\end{eqnarray}
\subsubsection*{Problems:}
\begin{enumerate}
\item Write a Python function to compute derivatives using these three finite differencing methods. Compute the derivative of the function $f(x) = e^x \, \sin(x)$ over the range $x = [0, 2\pi)$. Plot the numerically computed derivative $f_{(h)}'(x)$ for three different values of $h$.
\item Plot the error $\Delta f_{h}'(x) := |f_{h}'(x) - f'(x)|$ for three different values of $h$, and estimate the order of convergence $n$ of each finite-difference approximation. Plot $n(x)$, where
\begin{equation}
n(x) = \log_2 \frac{f'_{4h}(x) - f'_{2h}(x)}{f'_{2h}(x) - f'_{h}(x)}
\label{eq:order_convg}
\end{equation}
\item Reduce the step size $h$ successively. At what value of $h$ does the round-off error dominate the error budget?
%\item Derive a central differencing approximant for the first derivative $f'(x)$ that is accurate to $\mathcal{O}(h^4)$.
\end{enumerate}
\subsection{Richardson extrapolation}
We have seen that, if the order of the error in the numerical estimate of a function $f(x)$ is known, Richardson extrapolation provides a powerful way of improving the accuracy of the estimate. If we have two numerical estimates $f_{h}(x)$ and $f_{2h}(x)$ each having an error of $\mathcal{O}\,h^k$, a better estimate is given by
\begin{equation}
f(x) \simeq \frac{2 ~ 2^k\, f_{h}(x) - f_{2h}(x)}{2\,2^k -1 } + \mathcal{O}\,(h^{l}),
\end{equation}
where $l$ is the next-to-leading-order error term (for e.g., $l = k+2$ for central differencing, while $l = k+1$ for forward/backward differencing).
\subsubsection*{Problems:}
\label{sec:BBH_nr_data_fd}
\begin{enumerate}
\item Gravitational-waves (GWs) have two independent polarization states -- called ``plus'' and ``cross'' states. GW signals from the coalescence of black-hole binaries, in the simplest case, are circularly polarized:
\begin{eqnarray}
h_+(t) & = A(t) \, \cos \varphi(t), ~~~ h_\times(t) & = A(t) \, \sin \varphi(t).
\end{eqnarray}
Download the data file~\cite{nrdata} containing $h_+(t)$ and $h_\times(t)$. (This is the reduced form of the data produced by a numerical-relativity simulation of black-hole binaries performed by the SXS collaboration and is publicly available at~\cite{SXScatalog}). Compute the phase evolution $\varphi(t)$, the frequency evolution $\omega(t) := d\varphi(t)/dt$ and the rate of change of frequency $\dot{\omega}(t) := d\omega(t)/dt$ using second-order central difference approximation.
\item Estimate the order of convergence of the numerical computation of $\omega(t)$ and $\dot{\omega}(t)$.
\item Perform an extrapolation of $\omega(t)$ and $\dot{\omega}(t)$ to the next order using estimates of two different time-resolutions.
\item Derive an explicit expression for $f'(x)$ with error $\mathcal{O}(h^4)$ using Richardson extrapolation. This is the fourth-order finite differencing approximant for the derivative, which we will use later.
\end{enumerate}
fourier.tex 0 → 100644
\subsection{Fast Fourier transform (FFT)}
\subsubsection{Problems:}
\begin{enumerate}
\item Prove the convolution theorem.
\item Prove the correlation theorem.
\item Compute the Fourier transform of the function $h(t) = \sin 2 \pi f_0 t$ where $f_0 = 10$ Hz, $t = [0, 30]$ seconds using numpy.fft. Use different window functions and compare the Fourier transforms obtained using different window functions.
\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 (SINP).}
\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}
integration.tex 0 → 100644
\subsection{Numerical integration}
We derived the following formulas for evaluating the following integral numerically: $I = \int_{x_1}^{x_2} f(x) \, dx.$
\begin{eqnarray}
\mathit{Midpoint~rule:} & I = h \, f(\frac{x_1+x_2}{2}) + \mathcal{O}\,(h^3) \\
\mathit{Trapezoidal~rule:} & I = \frac{h}{2} \,\left[ f(x_1) + f(x_2) \right] + \mathcal{O}\,(h^3) \\
\mathit{Simpson's~1/3~rule:} & I = \frac{h}{3} \,\left[ f(x_0) + f(x_2) + 4 f(x_1) \right] + \mathcal{O}\,(h^5),
\end{eqnarray}
where $h := x_2 - x_1 = x_1 - x_0$ is the width of one slice. Note that the error estimates given above are the local errors. The global error is proportional to the local error $/~h$.
\subsubsection*{Problems:}
\begin{enumerate}
\item Show that Trapezoidal rule provides an estimate of the integral with local error $ \mathcal{O}\,(h^3)$.
\item Show that Simpson's 1/3 rule provides an estimate of the integral with local error $ \mathcal{O}\,(h^5)$.
\item Compute the following integral numerically using the two explicit methods given above, as well as using Romberg's method:
\begin{equation}
I = \int_{-1}^{1} e^{-x^2} dx.
\end{equation}
You may use the following SciPy functions: \href{http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.cumtrapz.html#scipy.integrate.cumtrapz}{\texttt{scipy.integrate.cumtrapz}}, \href{http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.integrate.simps.html#scipy.integrate.simps}{\texttt{scipy.integrate.simps}} and \href{https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.romb.html}{\texttt{scipy.integrate.romb}}.
\item Repeat the calculation using three different values of $h$. Plot the error $\Delta I$ against $h$. Show that the numerical estimates converge to the exact value according to the expected order of convergence.
\item Using the Trapezoidal and Simpson's methods, compute the energy loss due to GW emission from binary black holes using the numerical-relativity data discussed in Sec.~\ref{sec:BBH_nr_data_fd}. The radiated energy can be computed from the GW polarizations as
\begin{equation}
E = \int_{-\infty}^{\infty} \left[\left(\frac{dh_+}{dt}\right)^2 + \left(\frac{dh_\times}{dt}\right)^2 \right] \, dt.
\end{equation}
\item Repeat the calculation with three different values of $h$ (sampling rate $dt$ of the data). Compute the order of convergence [see, e.g., Eq.\eqref{eq:order_convg}].
\end{enumerate}
mc.tex 0 → 100644
\subsubsection{Problems:}
\begin{enumerate}
\item The area of a circle of radius $r$ is $A_c = \pi r^2$. The area of the smallest square that encloses this circle is $A_s = 4 r^2$. Measure the areas of the circle and square using a Monte-Carlo simulation and estimate the value of $\pi = 4 A_c/A_s$.
\item Repeat the problems 3 and 4 from Sec.~\ref{sec:statinf} by stochastically sampling the posterior using Markov-Chain Monte Carlo and then computing the evidence using Monte-Carlo integration. You can either code up the Metropolis-Hastings algorithm or use an existing sampler such as \href{https://dynesty.readthedocs.io/en/stable/}{dynesty.}
\end{enumerate}
\ No newline at end of file
ode.tex 0 → 100644
\subsection{Gravitational waves from inspiralling compact binaries}
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{https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.RK45.html}{\texttt{integrate.RK45}} routine. This uses the adaptive Runge-Kutta method by Dormand \& Prince that is very similar to the Runge-Kutta-Fehlberg method that we learned in the class. 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}
\end{enumerate}
%
\subsection{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}
For the problems in this section, you may use Scipy's high-level interface to various ODE solvers \href{https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html}{\texttt{solve\_ivp}}.
\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. These coupled ordinary differential equations are
\begin{eqnarray}
\frac{dx(t)}{dt} & = & \sigma [y(t)-x(t)], \nonumber \\
\frac{dy(t)}{dt} & = & x(t) [\rho - z(t)] - y(t), \nonumber \\
\frac{dz(t)}{dt} & = & x(t)y(t) - \beta z(t),
\end{eqnarray}
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 forces (e.g., random collissions of molecules) affecting the particle.
\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}
ode_bvp.tex 0 → 100644
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 0 → 100644
\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 0 → 100644
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}
statinf.tex 0 → 100644
Given some data $d$ and a model $\mathcal{M}$, Bayesian parameter estimation involves computing the posterior distribution of the parameters $\theta$ describing the model. Using Bayes theorem, we can write
\begin{equation}
p(\theta | d, \mathcal{M}) = \frac{p(\theta | \mathcal{M}) ~ p(d | \theta , \mathcal{M})}{p(d | \mathcal{M})},
\end{equation}
where $p(\theta | \mathcal{M})$ is the prior distribution of the parameters $\theta$, $p(d | \theta , \mathcal{M})$ is the likelihood of data given the parameters $\theta$ and the model $\mathcal{M}$, while
\begin{equation}
p(d | \mathcal{M}) = \int p(\theta | \mathcal{M}) ~ p(d | \theta , \mathcal{M}) \, d\theta
\label{eq:evidence}
\end{equation}
is the evidence of the model $\mathcal{M}$. Bayesian model selection involves comparing the evidence of different models, say $\mathcal{M}_1$ and $\mathcal{M}_2$, by means of the likelihood ratio (Bayes factor) between the two models.
\begin{equation}
\mathcal{B}^1_2 = \frac{p(d | \mathcal{M}_1)}{p(d | \mathcal{M}_2)}.
\end{equation}
\subsubsection{Problems}
Here we solve the curve fitting problem presented in Sec.~\ref{sec:curve_fitting} using of Bayesian statistical inference.
\begin{enumerate}
\item Compute the posterior distribution $p(H_0 | d, \mathcal{M}_1)$ of the Hubble constant $H_0$ using the Supernova Cosmology project data $z \leq 0.1$. Assume the Hubble's law [Eq.~\eqref{eq:Hubble_law}] as the model $\mathcal{M}_1$. Assume uniform prior for $H_0$ in the interval $(10, 100)$ km/s/Mpc. You can either use a Gaussian likelihood
%with $\sigma_{i} = 1$ or calculate the standard deviations (${\sigma_{i}}$) of $d_{L}$ from the data. In any case, the likelihood of data is:
\begin{equation}
p(d | H_0, \mathcal{M}_1) = \exp \left( - \sum_i \Big |\dfrac{d_i - d_L(z_i, \mathcal{M}_1)}{2\sigma_{i}^2} \Big |^2 \right),
\end{equation}
where $d_i$ and $z_i$ are the samples of the luminosity distance and redshift from the data and $d_L(z_i, \mathcal{M}_1)$ is the relation between luminosity distance and redshift predicted by model $\mathcal{M}_1$ evaluated at redshift $z_i$. $\sigma_{i}$ is the standard deviation of the measurement in $d_L$ which can be calculated using the error in the distance modulus $\mu$ (given in the data file~\cite{sndata}).
\item Repeat the analysis using the full data set. What are the differences that you see in the posterior?
\item Compute the posterior distribution $p(H_0, \Omega_M | d, \mathcal{M}_2)$ of the Hubble constant $H_0$ and matter density $\Omega_M$ using the $\Lambda$CDM model $\mathcal{M}_2$ [Eq.~\eqref{eq:lcdm}]. You can compute the posterior on a 2-dimensional grid. Assume uniform priors for $H_0$ in the interval $(10, 100)$ km/s/Mpc and for $\Omega_M$ in the interval (0, 1).
\item Compute the likelihood ratio (Bayes factor) between the Hubble's law and $\Lambda$CDM model by computing the evidences [Eq.~\eqref{eq:evidence}] of the two models using a numerical integration method that we learned in Sec.~\ref{sec:integr}.
\end{enumerate}
\ No newline at end of file
table2.tex 0 → 100644
This source diff could not be displayed because it is too large. You can view the blob instead.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or sign in to comment