Page MenuHomec4science
Files F93764295

task1.tex
No OneTemporary
Actions

File Metadata

Created
Sun, Dec 1, 07:40

task1.tex
View Options

\documentclass[a4paper,DIV14, BCOR0.5cm]{scrreprt}
\setlength{\parindent}{0mm}
\setlength{\parskip}{12pt}
\pagestyle{empty}
\usepackage{amssymb}
\usepackage{times}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{verbatim}
\newcommand{\bo}[1]{\boldsymbol{#1}}
\usepackage{enumitem}
\setenumerate{leftmargin=*,font=\bfseries}
\begin{document}
\begin{center}
\rule[2mm]{13.5cm}{.1pt} \\
{\Large {\textbf{Computer simulation of physical systems I}}}\\
\rule[2mm]{13.5cm}{.1pt}
\end{center}
\vspace{-1cm}
\begin{center}
{\textbf{Task I: numerical integration algorithms}} \\
\end{center}
Integration of one-dimensional harmonic oscillator
\begin{equation*}
\ddot{x} = -x
\end{equation*}
using Euler, predictor-corrector, and Verlet integrators,
and analysis of errors.
\begin{enumerate}
\item Go to the course homepage.
Download the archive for \verb!Task 1.zip! and unpack it.
The archive contains two files: \verb!solve_osc1d.py! and \verb!osc1d.py!.
Look at the code and try to understand it.
\item You can run python scripts by passing the name of the file to the
python interpreter\\ (ex. \verb!python solve_osc1d.py!). Or you can make the
script executable \verb!chmod +x solve_osc1d.py! and run it \verb!./solve_osc1d.py!.
\item Direct the output to a file
\verb!./solve_osc1d.py > output! and plot using any software of your choice.
You can also used included gnuplot script \verb!gnuplot plotter.p! or
functions defined in \verb!osc1d.py!.
The output contains calculated trajectories, exact trajectories and energy
($E(t) = x^2(t)+v^2(t)$).
\item Now it's your turn.
Using the method of your choice, calculate the following error measures
\begin{itemize}
\item Absolute error, i.e. maximum deviation, of the calculated
trajectory with respect to exact trajectory
${\rm max}_t |x_{calc}(t)-x_{exact}(t)| $
\item Energy drift: fit $E(t)$ to a straight line and drift is
then the slope of the line (you can use gnuplot script \verb!fitter.p!
or function defined in \verb!osc1d.py!)
\item Noise in the conservation of energy: calculate r.m.s. deviation
of $E(t)$ from the fitted line
$\sqrt{\frac{1}{n-2}\sum_n(\Delta E)^2}$
\end{itemize}
For the given example, these numbers should be about
1.009, 0.290, and 0.2.
\item You can see how predictor-corrector method improves Euler integrator. Use
\verb!euler_pred_corr(x,v,t)! function to perform integration. Analyze the
errors and the energy drift.
\item We can improve Euler integrator by
adding the second order term to the Taylor expansion of
the trajectory $\frac{1}{2} (\Delta t)^2 f(t)$ (see e.g. equation (3.8)
in the "Documentation: integrators in MD" in Moodle).
Implement this function in \verb!osc1d.py! and observe what happens.
Next, implement Verlet (optionally velocity Verlet) integrator as presented in the lecture.
Compare the error measures calculated above among the different integrators.
You should find out that Euler is bad, predictor-corrector
has a very small noise, and Verlet wins for absolute error and
for the drift.
\item Plot the above errors for different integrators as
a function of $\Delta t$.
%Note that you can just give the number
%of steps after the executable as e.g. \verb!./harmosc_euler 80!.
%For a few timesteps this can be done manually,
%or a Bash-script \verb!rundt.sh! is also included
%showing how to automatize a loop over a number of steps.
\item [Optional:] For a longer simulation find out the frequency
as a function of $\Delta t$.
\item [Optional:] Now consider damping in the 1D harmonic oscillator
\begin{equation*}
\ddot{x} + \gamma \dot{x} + \omega_0^2 x = 0
\end{equation*}
where $\gamma$ is the damping coefficient,
and $\omega_0$ is the undamped angular frequency of the oscillator.
Implement damping in your code using the Euler, pred/corr, and Verlet algorithms
(take the same initial conditions as that of the earlier problem).
Vary the ratio ($\gamma/\omega_0$) and see the evolution of the damping behavior.
\begin{comment}
\item [Optional:] \textbf{Particle in a box}
\item[*] {Finite difference method} \\
Using the Euler method, show that for $f(x)$, its second derivative can be obtained from
\begin{equation*}
\frac{d^2f(x)}{dx^2}|_{x=x_n} = \frac{1}{a^2}[f(x_{n+1})-2f(x_n)+f(x_{n-1})]
\end{equation*}
where $a$ is the step size $x_n-x_{n-1}$.
\item[*] {Schr\"{o}dinger equation} \\
The one-dimensional time-independent Schr\"{o}dinger equation (in a.u.)
\begin{equation*}
\hat{H}\psi(x) = \varepsilon \psi(x),
\hat{H} = -\frac{1}{2}\frac{d^2}{dx^2} + V(x)
\end{equation*}
can be written in the matrix representation
, where the Hamiltonian $H$ is a $N \times N$ matrix,
and the eigenfunction $\psi$ is a $N \times 1$ vector ($N$ is the number of discrete points).
Using the finite difference method, we show that $\hat{H}$ has the following form
\[ \left( \begin{array}{ccccc}
2t_0+V_1 & -t_0 & \cdots & 0 & 0 \\
-t_0 & 2t_0+V_2 & \cdots & 0 & 0 \\
\vdots & \vdots & \ddots & 0 & 0 \\
0 & 0 & \cdots & 2t_0+V_{N-1} & -t_0 \\
0 & 0 & \cdots & -t_0 & 2t_0+V_{N} \end{array} \right)\]
where $t_0=\frac{1}{2a^2}$.
\item[*] Suppose now a particle is trapped in a box with a width $L = 100$ (a.u.).
The external potential is zero inside the box ($0 \le x \le L$), and is infinite outside.
Use 100 lattice points, and calculate the eigenenergies for the particle by diagonalizing $\hat{H}$,
and compare the numerical results to the exact dispersion $E_k = \frac{\pi^2 k^2}{2L^2}$
($k=1,2,3,\cdots,N$).
Figure out why the finite difference method deviates from the exact solution at higher energies.
Plot the probability distribution ($|\psi_k|^2$) for the two modes ($k=1$ and $k=5$).
\begin{figure}[hb]
\centering
\includegraphics[scale=0.45]{box.eps}
\end{figure}
\end{comment}
\end{enumerate}
\end{document}

Event Timeline

c4science · Help