# FYS-4096 Computational Physics

## Introduction to the Art of Numerical Experiments

Week 1: Basics of scientific computing
Course contents
Linux 101
Python refresher
Crash course on Git
Numerical experiments workflow
Week 2: Numerical differentiation, integration, and discrete fourier transformation
Python packages for Scientific Computing
Numerical differentiation
Discrete Fourier transformations
Differentiation with discrete fourier transformation
Numerical integration
Week 3: Visualization
Root finding
Visualization
Week 4: Linear algebra
Numerical linear algebra 101
System of linear equations
LU-decomposition
Eigenvalues of a dense matrix
QR-algorithm
Multidimensional root finding
Week 5: Linear algebra part 2
Sparse matrices
Matrix functions
Krylov subspace 101
Matrix exponential in Krylov subspace
Week 6: Quantum mechanics with finite differences
Short note on units
Quantum Mechanics in 25 minutes
Discretizing PDEs with finite differences — time-independent Schrödinger equation
Discretizing time-dependent PDEs with finite differences — time-dependent Schrödinger equation
Structure of finite-difference TDSE solvers
Week 7: Supercomputers and -clusters
Computing clusters
Transferring and storing files
Accessing Taito
Modules
Interactive jobs
Running simulations on the cluster
Tips and tricks to boost your workflow
Week 8: PDEs and weak formulation
Partial differential equations, part 2
Weak formulation of PDEs
Week 9: Finite Element Method (FEM)
FEM basics
Meshing
Approximating the function spaces
FEM with FEniCS
Example: a boundary value problem
Example: an eigenvalue problem
Other FEM software: Comsol and Ansys
Week 10: Ordinary differential equations
ODE types and reduction of order
Implicit and explicit solvers
Runge Kutta methods
Conservation laws in classical mechanics and numerical methods
Example: Molecular dynamics
Week 11: PDEs and initial value problems
Problem formulation
Stiff ODEs
Difficult linear systems with preconditioners
Week 12: Time-series analysis
Basic terminology
Detrending
Spectral analysis
Correlations
Filtering
Searching patterns
Week 13: Testing, debugging, profiling, and optimizing numerical codes
Verifying correctness of codes
Debugging
How to ask for help and report a bug?
Code optimization in HPC
Week 14: Optimization
Optimization
Curve fitting
Feedforward neural networks

# Course contents

Linux basics Working with remote servers Working with supercomputers, Vim
Visualization Perceptually uniform colormaps Publication-quality figures
Setting up your numerical experiment a.k.a.
"Good Enough Practices in Scientific Computing"
"Best Practices for Scientific Computing" Archiving and publishing your numerical experiments
Numerical calculus Multi-dimensional calculus Advanced methods
Numerical solution of ODEs Examples in classical mechanics Algorithm selection
Numerical linear algebra Multi-dimensional linear algebra,
preconditioning,
examples in time-independent
quantum mechanics
Calculating matrix functions
Solving partial differential equations: Finite difference method, weak formulation, Finite element method Examples in thermodynamics, electrostatics,
and time-dependent quantum mechanics
More complex PDEs
Correlations and spectral analysis Signal decomposition Behind the noise
Machine learning Neural networks Pitfalls and catastrophic failures

# Passing the course

### Lessons

Mondays 12-14
Lectures (attendance optional)
Wednesdays 12-15
Hands-on tutorials and help with exercises (attendance optional)
Fridays 10-12
Exercise session (attendance required)

### Exercises (1/2)

When you get them?
On Mondays after lectures
Where you get them?
Course GitLab page at gitlab.com
When you complete them?
Some effort before Wednesdays' tutorial session
During Wednesdays' tutorial session
Rest before Friday

### Exercises (2/2)

When do return them?
Friday before 5 am
How do you return them?
To your own repository at gitlab.com with a git tag.
See the first exercise set for detailed instructions.
Is attendance mandatory at Friday's exercise session?

### Project work

How many?
2
When are topics released?
February 19th
April 16th
March 12th 23:59
May 14th 23:59
More details?
TBA

### Feedback

During the course!
During lecture, after tutorials, before exercises, all the time!

### Final feedback

Leaving feedback in KAIKU-system is mandatory. Empty feedback possible.

# Assessment

Maximum of 280 experience points

 Exercise sets 14 × 10 XP 50 % 1st project 70 XP 25 % 2nd project 70 XP 25 %

XP % of max. XP Grade
140 50 % 1
168 60 % 2
196 70 % 3
224 80 % 4
252 90 % 5

# Installation instructions for this course (1/2)

2. Select tab VMware
3. Select product corresponding to your operating system:
Windows ↔ VMware Workstation 14
macOS ↔ VMware Fusion 10
Linux ↔ VMware Workstation 14

# Installation instructions for this course (2/2)

1. Get the prepared virtual machine image for this course from https://www.tut.fi/fys/fys4096/
2. Import the image to VMware and launch it.
3. You should automatically login to your virtual machine with Ubuntu 17.10
4. In case you need these: the username is 'student' and the password 'compphys'. Change them if you wish.

# Shell 101

• Navigating inside directory hierarchy
• Basic shell commands
• Executables
• Shell scripts
• Environment variables
• Manual pages

Live demo

# IPython

Interactive Python


$ipython3 Python 3.6.2 (default, Jul 17 2017, 16:44:45) Type 'copyright', 'credits' or 'license' for more information IPython 6.1.0 -- An enhanced Interactive Python. Type '?' for help. In [1]: print( "Hello World" ) Hello World In [2]:  # Python source files$\exp(x) \approx \sum\limits_{k=0}^{n}\frac{1}{k!}x^k$# Python testing # Python package structure . ├── README.rst ├── license.txt ├── requirements.txt ├── setup.py └── taylor_approximation ├── __init__.py ├── exp.py └── tests ├── __init__.py └── test_exp.py # File: README.rst • Short description of the package contents # File: license.txt # File: requirements.txt • Details dependencies. One per line. • No dependencies in our trivial example # File: setup.py • Setup for the python package installation. # File: taylor_approximation/__init__.py • Let's Python know that this directory contains a package • Contains top-level documentation of the package # File: taylor_approximation/exp.py • First actual source code file # File: taylor_approximation/tests/__init__.py • Let's Python know that this directory contains a package • This directory contains all the tests of the package. # File: taylor_approximation/tests/test_exp.py • Actual tests • Our setup.py automatically finds all these tests # Version control with git • Tracks changes in your source files • You save snapshots of your project, called commits • Nonlinear history possible # Git example # Automate! ## If you find yourself doing some repetitive task multiple times, automate it! For example, • Conformance to a programming style guide$\rightarrow$yapf, autopep8, pycodestyle, clang-format • Running multiple simulations with slightly different parameters$\rightarrow$a bash script • Installing a new version of in-house software to supercomputers$\rightarrow$CI/CD, containers • Similar project or document structures$\rightarrow$file templates • Working-hour allocations to online systems$\rightarrow$a Python script # Numerical experiments workflow # Organizing your numerical experiment • Use two or more git repositories: • One for software development (if needed) • The rest for your numerical experiments with said package • Automate and document the entire workflow $ ./install_project_dependencies.sh
$./run_my_experiments.py $ ./generate_analysis_data.py
$./generate_intermediate_figures.py $ ./generate_publication_figures.py
• Document the format, data gathering methods etc. for all files generated manually or by someone else.
• For different versions, use git tags and command-line options.
• DO NOT DO THE FOLLOWING:

.
├── data.txt
├── data_v2.txt
├── data_final.xlsx
├── some_data.txt
├── simulation1.py
├── simulation2.py
├── simulation3.py
├── simulation_finalversion2.py
├── simulation_finalversion_minorfix.py
├── simulation_finalversion.py
└── simulation_not_working.py

# Ideal workflow

First figure our what you need to model and how you need to model it. Then

1. Interactive testing and development with IPython.
 >>> history
2. Save working piece of code to a source file
3. Create unit tests for your new piece of code
4. Create larger tests for the entire package (e.g., test agains some analytically solvable cases etc.)
5. Now you have a first version of the software. Use it to run some simulations either interactively via IPython or via Python scripts.
6. Write and test Python scripts for generating your data and analysis/figures.

# Typical workflow

First figure our what you need to model and how you need to model it. Then

1. Interactive testing and development with IPython.
 >>> history
2. Save working piece of code to a source file
3. Create unit tests for your new piece of code
4. Create larger tests for the entire package (e.g., test agains some analytically solvable cases etc.)
5. Now you have a first version of the software. Use it to run some simulations either interactively via IPython or via Python scripts.
6. Write and test Python scripts for generating your data and analysis/figures.
7. GOTO Step 1

# Rules of Thumb (design)

• Modular design makes expansion easy
• Decompose program to smaller and smaller functions
A function should perform only one task
• Enforce interfaces and preconditions
• Encapsulate related data and functions into classes
• No code duplication (use functions/classes)
• No reinventing the wheel: Use existing libraries
• Smart naming
• Pick a coding style and stick to it
• Do not comment/uncomment sections of the code to change the simulation
• Automatic testing

# Modular design and encapsulation

• In Computational Physics, modularize using:
• Physical objects
(particle, car, EM-wave, quantum wavefunction, ...)
• Theoretical framework
[Interactions (forces) + Newton's equations, electric field, magnetic field, vector potential, ...)
• Mathematical structures
(functions, equations, equation solvers, function spaces, inner product, vector space, ...)

# Encapsulated design example

Note: Bad design for other reasons (vectorization).

# Numpy

• Data-types and fast N-dimensional arrays
• Basic functions for operating with NumPy arrays

# Numpy arrays

• Single- or multidimensional arrays
• Indexing starts at 0
• Numpy's and Scipy's functions also operate on numpy arrays

# Scipy

• Algorithms and methods

• 2D plots

# Mayavi

• 3D plots
Link to Mayavi examples (embedding disabled due to mistakes in their SSL certs)

# h5py — saving data

• Python interface for HDF5-fileformat
• Meta-data in the same file

# Numerical Differentiation

$$f^\prime(x_0)=\lim\limits_{x_1\to x_0}\frac{f(x_1)-f(x_0)}{x_1-x_0} \approx \frac{y_1-y_0}{x_1-x_0}$$

## Finite-differences

• $f^\prime(x_0) \approx \dots$ discuss and give me your ideas!
• $f^\prime(x_0) \approx \displaystyle \frac{f(x_0+h)-f(x_0)}{h}$
• $f^\prime(x_0) \approx \displaystyle \frac{f(x_0)-f(x_0-h)}{h}$
• $f^\prime(x_0) \approx \displaystyle \frac{f(x_0+h)-f(x_0-h)}{2h}$
• $f^\prime(x_0) \approx \frac{1}{h} \left[\frac{1}{12} f(x_0-2h)-\frac{2}{3} f(x_0-h) + \frac{2}{3}f(x_0+h) - \frac{1}{12}f(x_0+2h) \right]$
• $f^{\prime\prime}(x_0) \approx \frac{1}{h^2} \left[ f(x_0-h) - 2 f(x_0) + f(x_0+h) \right]$

#### Course book

Ch. 5.1-5.6: Differentiation

## Deriving a finite-difference approximation

Suppose we know the value of function $f(x)$ at $x_0-h, x_0,$ and $x_0+h$. Let's derive a second-order finite-difference approximation for the first derivative.

First expand as Taylor series: \begin{align} f(x+h) &= f(x) + f'(x)h + \frac{1}{2}f''(x) h^2 + \mathcal{O}(h^3)\\ f(x-h) &= f(x) - f'(x)h + \frac{1}{2} f''(x) h^2 + \mathcal{O}(h^3) \end{align}

Multiply the first by $a$ and the second one by $b$, add them up, and solve for $f'(x)$. \begin{align} a f(x+h) + b f(x-h) &= (a+b) f(x) + (a-b) f'(x) h + \frac{a+b}{2} f''(x) h^2 + \mathcal{O}(h^3) \end{align}

\begin{align} \Leftrightarrow f'(x) &= \frac{a f(x+h) - (a+b) f(x) + b f(x-h) }{(a-b) h} - \frac{ a+b}{2(a-b)} f''(x) h + \mathcal{O}(h^2) \end{align}

The first order error term vanishes when $a=-b$, i.e., $$f'(x) = \frac{f(x+h) - f(x-h)}{2h}+\mathcal{O}(h^2)$$

## Discrete Fourier transformation part 1

Consider a $L$-periodic function $f(x)$, $x \in [0,L]$ with Fourier series $$f(x) = \frac{1}{L}\sum\limits_{k=-\infty}^\infty F_k \exp\left(\frac{2\pi \mathrm{i}}{L}k x\right), \text{where}$$ $$F_k = \int\limits_0^L \exp\left( - \frac{2\pi \mathrm{i}}{L}k x \right) f(x)\,\mathrm{d}x, k\in \mathbb{Z}$$

On a computer, we have only a finite number of samples (consider only equidistant samples here) $$\{ f(x_0), f(x_1), \dots, f(x_{N-1})\} \doteq \left\{ f_0, f_1, \dots, f_{N-1}\right\}$$ so we do a Discrete Fourier Transform to approximate the Fourier coefficients:

#### Course book

ch. 12: Fourier Analysis

\begin{align} F_k &=\int\limits_0^L \exp\left( - \frac{2\pi \mathrm{i}}{L}k x \right) f(x)\,\mathrm{d}x \\ & \approx \sum\limits_{n=0}^{N-1} f(x_n) \exp\left( - \frac{2 \pi \mathrm{i}}{L} k x_n \right) \Delta x \\ & = \sum\limits_{n=0}^{N-1} f(x_n) \exp\left( - \frac{2 \pi \mathrm{i}}{N \color{red}{\Delta x}} k \color{red}{\Delta x} \cdot n \right) \Delta x \\ & = \frac{L}{N} \sum\limits_{n=0}^{N-1} f_n \exp\left( - \frac{2 \pi \mathrm{i}}{N} k n \right) \\ & = \frac{L}{N} \mathcal{F}^\mathrm{DFT}_k,\, k \in \left\{0,\dots,N-1\right\} \end{align}

$\mathcal{F}^{DFT}_k$ is a discrete Fourier Transform coefficient which is typically calculated using the Fast Fourier Transform algorithm (see next slide for an example).

## Discrete Fourier transformations part 2

To obtain the original samples, we just need to do an inverse discrete fourier transformation: $$f(x_n) = \frac{1}{L} \sum\limits_{k=0}^{N-1} F_k \exp\left(\mathrm{i} \frac{2 \pi}{N}n k\right) = \frac{1}{N} \sum\limits_{k=0}^{N-1} \mathcal{F}_k^\text{DFT} \exp\left( \frac{2\pi \mathrm{i}}{N} k n\right) = \mathcal{IF}_n^\text{DFT}\left[ \left\{ \mathcal{F}_k^\text{DFT} \right\}_{k=0}^{N-1} \right]$$

## Discrete Fourier transformation, frequencies

We had \begin{align} F_k &=\int\limits_0^L \exp\left( - \frac{2\pi \mathrm{i}}{L}k x \right) f(x)\,\mathrm{d}x \\ & \approx \frac{L}{N} \sum\limits_{n=0}^{N-1} f_n \exp\left( - \frac{2 \pi \mathrm{i}}{N} k n \right) \\ & = \frac{L}{N} \mathcal{F}^\mathrm{DFT}_k,\, k\in \left\{0, \dots, N-1\right\} \end{align}

This means that the frequency component $\omega_k = 2 \pi f_k = 2 \pi \frac{k}{L}$ of the function $f(x)$ is given by $F_k$ or equivantly $\mathcal{F}^\mathrm{DFT}_k$.

## Discrete Fourier transformation, negative frequencies

In the original continuous-time Fourier series we also had negative frequencies: $$f(x) = \sum\limits_{k=-\infty}^\infty F_k \exp\left( -\frac{2\pi \mathrm{i}}{L} k x\right)$$

The negative frequencies in our discrete transformation are given by: \begin{align} F_{-m} &= \frac{L}{N} \sum\limits_{n=0}^{N-1} f_n \exp\left( - \frac{2 \pi \mathrm{i}}{N} (-m) n \right) = \frac{L}{N} \sum\limits_{n=0}^{N-1} f_n \exp\left( - \frac{2 \pi \mathrm{i}}{N} (-m) n \right)\cdot1^n \\ &= \frac{L}{N} \sum\limits_{n=0}^{N-1} f_n \exp\left( - \frac{2 \pi \mathrm{i}}{N} (-m) n \right)\cdot \exp\left(-2\pi \mathrm{i}\right)^n = \frac{L}{N} \sum\limits_{n=0}^{N-1} f_n \exp\left( - \frac{2 \pi \mathrm{i}}{N} (-m) n \right)\cdot \exp\left(-2\pi \mathrm{i} n\right) \\ &= \frac{L}{N} \sum\limits_{n=0}^{N-1} f_n \exp\left( - \frac{2 \pi \mathrm{i}}{N} (-m) n-2\pi\mathrm{i} n \right) = \frac{L}{N} \sum\limits_{n=0}^{N-1} f_n \exp\left( - \frac{2 \pi \mathrm{i}}{N} (N-m) n \right)\\ &= F_{N-m} \end{align}

# Order of Fourier coefficients in NumPy(and FFTW)

Real space

Fourier space:

## Discrete fourier transformation and interpolation

Differentiation needs a continuous function. How to get that from DFT? No unique way.

Straightforward method is stupid: $$f(x) \approx \frac{1}{L}\sum\limits_{k=0}^{N-1} F_k\exp\left(\frac{2\pi \mathrm{i}}{L}k x\right)$$

## Discrete fourier transformation and interpolation, part 2

For the best approximation use the property that $F_{N-k}=F_{-k}$:

$$f(x) \approx F_0 + \frac{1}{L}\sum\limits_{k=1}^{\left\lfloor \frac{N}{2}\right\rfloor} \left[ F_k\exp\left(\frac{2\pi \mathrm{i}}{L}k x\right) + F_{N-k} \exp\left(-\frac{2\pi \mathrm{i}}{L} k x\right) \right] + F_{\frac{N}{2}} \cos\left( \frac{2 \pi}{L} \cdot 2N x \right)$$

where the last term, $F_{N/2}$, is dropped when $N$ is even.

## Discrete fourier transformation and differentiation

Use our interpolation to calculate the derivative: \begin{align} \frac{\mathrm{d}}{\mathrm{d}x} f(x) & \approx \frac{1}{L}\sum\limits_{k=1}^{\left\lfloor \frac{N}{2}\right\rfloor} \left[ \left( \frac{2 \pi \mathrm{i}}{L} k \right) F_k\exp\left(\frac{2\pi \mathrm{i}}{L}k x\right) + \left( -\frac{2 \pi\mathrm{i}}{L}k \right) F_{N-k} \exp\left(-\frac{2\pi \mathrm{i}}{L} k x\right) \right] \\ &- \left( \frac{\pi}{L} N \right) F_{\frac{N}{2}} \sin\left( \frac{2 \pi}{L} \cdot \frac{N}{2} x \right) \end{align} where the last term, $F_{N/2}$, is dropped when $N$ is even.

Next, evaluate at sample points (since we're on a computer): \begin{align} \frac{\mathrm{d}}{\mathrm{d}x} \Bigg\vert_{x=x_n} f(x) & \approx \frac{1}{L}\sum\limits_{k=1}^{\left\lfloor \frac{N}{2}\right\rfloor} \left[ \left( \frac{2 \pi \mathrm{i}}{L} k \right) F_k\exp\left(\frac{2\pi \mathrm{i}}{L}k x_n\right) + \left( -\frac{2 \pi\mathrm{i}}{L}k \right) F_{N-k} \exp\left(-\frac{2\pi \mathrm{i}}{L} k x_n\right) \right] \\ &\,\,\,\,\,\,\,- \left( \frac{\pi}{L} N \right) F_{\frac{N}{2}} \sin\left( \frac{2 \pi}{L} \cdot \frac{N}{2} x_n \right) \\ & = \frac{1}{L}\sum\limits_{k=1}^{\left\lfloor \frac{N}{2}\right\rfloor} \left[ \left( \frac{2 \pi \mathrm{i}}{L} k \right) F_k\exp\left(\frac{2\pi \mathrm{i}}{L}k \cdot \frac{L}{N}n\right) + \left( -\frac{2 \pi\mathrm{i}}{L}k \right) F_{N-k} \exp\left(-\frac{2\pi \mathrm{i}}{L} k \cdot \frac{L}{N} n\right) \right] \\ & = \frac{1}{N}\sum\limits_{k=0}^{N-1} c_k \mathcal{F}^\text{DFT}_k\exp\left(\frac{2\pi \mathrm{i}}{N}k \cdot n\right) \end{align}

\begin{align} \frac{\mathrm{d}}{\mathrm{d}x} \Bigg\vert_{x=x_n} f(x) & \approx \frac{1}{N}\sum\limits_{k=0}^{N-1} c_k \mathcal{F}^\text{DFT}_k\exp\left(\frac{2\pi \mathrm{i}}{N}k \cdot n\right) \\ & = \mathcal{IF}_n^\text{DFT}\left[ \left\{ c_k \mathcal{F}_k^\text{DFT} \right\}_{k=0}^{N-1}\right], \end{align} where \begin{align} c_k = \left\{ \begin{array}{ll} \frac{2\pi \mathrm{i}}{L} k & k < \frac{N}{2} \\ \frac{2\pi \mathrm{i}}{L} \left(k - N\right) & k > \frac{N}{2} \\ 0 & k=\frac{N}{2} \text{ if N even} \end{array} \right. \end{align}

# Numerical Integrals

## Trapezoidal rule

$$\int\limits_{x_0}^{x_0+N\cdot h}\!\!\!\!\!\! f(x)\,\mathrm{d}x \approx \frac{h}{2} \sum\limits_{i=0}^{N-1} \left[ f(x_0) + f(x_0+i\cdot h) \right]$$

#### Course book

Ch 5.9 Trapezoid rule

## Simpson's rule

$$\int\limits_{x_0}^{x_0+N\cdot h}\!\!\!\!\!\! f(x)\,\mathrm{d}x \approx \frac{h}{3} \sum\limits_{i=1}^{\frac{N-1}{2}} \left[ f(x_0+(2i-2)h) + f(x_0+(2i-1)\cdot h) + f(x_0+2i) \right],\,\text{where }N\text{ is odd}$$

#### Course book

Ch 5.10 Simpson's rule

## Improper integrals

$$\int\limits_{x_0}^\infty f(x)\,\mathrm{d}x = \, \text{#\@&%*!💣}$$

Change of variables, e.g., $x \to \arctan(x)$ yields: $$\int\limits_{x_0}^\infty f(x)\,\mathrm{d}x =\!\!\!\!\!\!\int\limits_{\arctan(x_0)}^{\pi/2}\!\!\!\! \frac{f[\tan(x)]}{\cos^2(x)}\,\mathrm{d}x = \!\!\!\!\!\!\int\limits_{\arctan(x_0)}^{\pi/2}\!\!\!\! g(x) \,\mathrm{d}x$$ The last form we can already integrate with standard techniques.

## Fubini's theorem and the Curse of Dimensionality

2D trapezoidal rule: \begin{align} \int\limits_{x_0}^{x_0+Nh}\int\limits_{y_0}^{y_0+Nh} f(x,y)\,\mathrm{d}x\,\mathrm{d}y & \approx \int\limits_{x_0}^{x_0+Nh} \frac{h}{2}\sum\limits_{j=0}^{N-1} [f(x, y_0+j h)+f(x, y_0+(j+1)h)]\,\mathrm{d}x\\ & \approx \frac{h^2}{4}\sum\limits_{i=0}^{N-1}\sum\limits_{j=0}^{N-1} [f(x_0 + ih, y_0+j h) +f(x_0+i h, y_0+(j+1)h) \\ & \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+f(x_0+(i+1)h, y_0+j h)+f(x_0+(i+1)h, y_0+(j+1)h)] \end{align} $\sim N^2$ operations

M-dimensional integral: $\sim N^M$ operations

## Monte Carlo integration

• Monte arlo = using random numbers
• High-dimensional integrals = ❤ Monte Carlo
• Many algorithms, e.g.,

#### VEGAS Monte Carlo

• Evaluate function where it contributes to the integral most
• Update sampling distribution $g( \overline{\mathbf x})$ every now and then
• Approximation: $g(\overline{\mathbf x}) = g_1(x_1) g_2(x_2)\cdots g_N(x_N)$
• If valid, VEGAS works well

#### Course book:

Ch 5.14-5.15: Monte Carlo Integration
Ch 5.20: Importance Sampling

$$\int\limits_0^3\int\limits_5^6\int\limits_{-5}^7 \frac{1}{1+\exp\left(-\Vert x\Vert^2\right)}\,\mathrm{d} \overline{\mathbf x}$$

# Root finding

## Bisection method= brute force

• Function must always have different signs at end points of each interval

## Newton's method

$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$

# Matplotlib

• 2D and simple 3D plots

# 2D figure design: Important parts

1. Axis labels
2. Axis units
3. Ticks
4. Data
5. Axis ratio

### 1. Axis labels

• Short but descriptive
• Prefer words over symbols (if possible)

### 2. Axis units

• Use standard units of the field (eV vs joules)

### 3. Ticks

• Not too many
• Prefer exact values if possible
• Use symbols if necessary ($\pi$,$e$, $10^5$)

### 4. Data

• Choose linestyle well (solid, dashed, dotted,...)
• Plot types: line, scatter, histogram, error-bars

### 5. Axis ratio

• Best default choise: $\frac{L_x}{L_y} = \frac{\sqrt{5}+1}{2} =$ golden ratio

# Some plot types

### Scatter plot

• Discrete data points
• Color coding
• Markers
• Error bars

### Line plot

• Near continuous data points
• Color coding
• Linestyles

### Scatter plot + errorbars

• Discrete data with errorbars
• Axis scales: lin/log
• Color coding
• Markers

### Line plot + shaded trust region

• Continuous data with errorbars
• Axis scales: lin/log
• Color coding
• Line-styles
• Opacity

### 1D histograms

• Samples from a distribution
• Axis scales: lin/log
• Color coding
• Line/box
• Normalized?

### Contour plots

• 2D function values
• Axis scales: lin/log
• Colormaps
• Colorscale (min/max, lin/log)
• Contour lines
• Filling

### Density plots

• 2D function values
• Axis scales: lin/log
• Colormaps
• Colorscale (min/max, lin/log)

### Quiver plots

• 2D vector fields
• Axis scales: lin/log
• Arrow magnitude (lin/log, min/max)
• Arror colors

### Streamlines

• 2D vector fields
• Axis scales: lin/log
• number of streamlines
• Streamline colors, widths, colormap, number of streamlines

### Surfaces

• Scalar field/data that depends on two variables

### Isosurfaces

• Scalar field of 3 variables
• 2D surfaces of constant function value
• Cf. contour plots

### Contour/density slices

• Scalar field in 3D

### Animations

• Animation adds an extra dimension (best for time, but ok also for others)


mencoder mf://videofiles/*.png -mf type=png:w=1920:h=1034:fps=25 -ovc lavc -lavcopts vcodec=mpeg4 -oac copy -o outfile.mp4


# Colormaps

• What does your data convey?
• Prefer perceptually uniform colormaps
• Design for color-blindness
• Design for printing in grayscale

# Figures in papers

• Same font as in main text
• Matching font size
• Fileformat for figures should be PDF. Only convert to EPS if publisher requires it.
• Avoid bitmaps at all cost
• All aspects of figures should be vector graphics
• Except when you have images or image-like objects. Then those parts may be rasterized, but axes, labels, ticks, etc. still vector graphics.
• Fileformat doesn't make it vector graphics!

# Numerical linear algebra 101

## Linear algebra refresher in $\mathbb{R}^N$

A vector space over the field $\mathbb{R}$ is a set $\mathbb{R}^N$ with two operations,

• $+"$ (addition) and
• $\cdot"$ (multiplication by scalar)

that satisfy
For any $\mathbf{\overline x},\mathbf{\overline y},\mathbf{\overline z}\in \mathbb{R}^N$ and $\alpha,\beta \in \mathbb{R}$ we set the following axioms

• $\mathbf{\overline x}+\mathbf{\overline y} = \mathbf{\overline y}+\mathbf{\overline x}$
• $\mathbf{\overline x}+ (\mathbf{\overline y}+\mathbf{\overline z}) = (\mathbf{\overline x}+\mathbf{\overline y}) +\mathbf{\overline z}$
• $\exists\, \mathbf{\overline 0} \in \mathbb{R}^N : \mathbf{\overline x} + \mathbf{\overline 0} = \mathbf{\overline x}$
• $\exists\, \mathbf{-\overline x}\in \mathbb{R}^N : \mathbf{\overline x} + (\mathbf{-\overline x}) = 0$
• $\alpha \cdot (\beta \cdot \mathbf{\overline x}) = (\alpha \beta)\mathbf{\overline x}$
• $1\cdot\mathbf{\overline x} = \mathbf{\overline x}$
• $(\alpha+\beta)\mathbf{\overline x} = \alpha \mathbf{\overline x} + \beta \mathbf{\overline x}$
• $\alpha ( \mathbf{\overline x} + \mathbf{\overline y} ) = \alpha \mathbf{ \overline x} + \alpha \mathbf{\overline y}$

## Normed vector spaces

• A normed space is a vector space $V$ equipped with a (vector) norm $\Vert \cdot \Vert$.
• Example: $\mathbb{R}^N$ with the norm $\Vert \mathbf{\overline x} \Vert=\sqrt{\sum\limits_{i=0}^{N-1} x_i^2}$ is a normed vector space.

## Inner product space

• An inner product space is a vector space $V$ equipped with an inner product $\langle \cdot \mid \cdot \rangle:V\times V\to\mathbb{R}"$.
• Example: $\mathbb{R}^N$ with the inner (dot) product $\langle \mathbf{\overline x} \mid \mathbf{\overline y} \rangle = \sum\limits_{i=0}^{N-1} x_i y_i$ is an inner product space.
• Note: All inner product spaces are normed vector spaces with the (natural) norm $\Vert \mathbf{\overline x} \Vert = \sqrt{\langle \mathbf{\overline x} \mid \mathbf{\overline x} \rangle }$

## Linear maps and matrices

• A linear map $L$ is a map between two vector spaces (e.g., from $\mathbb{R}^n \to \mathbb{R}^p$) that satisfies:
• $L(\mathbf{\overline x}+\mathbf{\overline y}) = L(\mathbf{\overline x}) + L(\mathbf{\overline y})\,\forall \mathbf{\overline x},\mathbf{\overline y}\in\mathbb{R}^n$
• $L(\alpha \mathbf{\overline x}) = \alpha L(\mathbf{\overline x})\,\forall \mathbf{\overline x} \in \mathbb{R}^n$ and $\forall \alpha \in \mathbb{R}$
• Every linear map (between finite spaces) can be represented by a corresponding matrix: $$\left[ \begin{array}{cccc} L_{00} & L_{01} & L_{02} & L_{03} \\ L_{10} & L_{11} & L_{12} & L_{13}\\ L_{20} & L_{21} & L_{22} & L_{23} \end{array}\right] \left[ \begin{array}{c} x_{0} \\ x_{1} \\ x_{2} \\ x_3 \end{array}\right] = \left[ \begin{array}{c} L_{00}x_{0} + L_{01}x_1 + L_{02}x_2 + L_{03}x_3 \\ L_{10}x_0 + L_{11}x_1 + L_{12}x_2 + L_{13}x_3\\ L_{20}x_0 + L_{21}x_1 + L_{22}x_2 + L_{23}x_3 \end{array}\right]$$

# Numerical linear algebra with Python

## Vectors

$$\mathbf{v}_1 = \left[ \begin{array}{ccc} 1 \\ 3 \\ -5 \end{array} \right], \,\,\,\,\,\,\, \mathbf{v}_2 = \left[ 1,\, 3,\, -5 \right]$$


# No distinction between column and row vectors
# when initialized like this:
v1 = np.array([ 1, 3, -5 ])

v2 = np.array([ 1, 3, -5 ])

# Distinction possible (if absolutely required)
v1 = np.array([
[ 1 ],
[ 3 ],
[-5 ]
])

v2 = np.array([ [1, 3, -5] ])


## Matrices

$$\mathbf{M} = \left[ \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 5 & -1 \\ 7 & 1+\mathrm{i} & -\pi \end{array} \right]$$


M = np.array([
[ 1,    2,  3 ],
[ 1,   -5, -1 ],
[ 7, 1+1j, -np.pi]
])

# or

M = np.matrix([
[ 1,    2,  3 ],
[ 1,   -5, -1 ],
[ 7, 1+1j, -np.pi]
])



## Matrix-vector product

$$\left[ \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 5 & -1 \\ 7 & 1+\mathrm{i} & -\pi \end{array} \right] \left[ \begin{array}{c} 1 \\ 3 \\ -5 \end{array} \right] = \left[ \begin{array}{c} -8 \\ 21 \\ 10+5\pi+3\mathrm{i} \end{array} \right]$$ $$\left[ 1, 3, -5 \right] \left[ \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 5 & -1 \\ 7 & 1+\mathrm{i} & -\pi \end{array} \right] = \left[ -31,\, 12-5\mathrm{i},\, 5\pi\right]$$


M = np.array([
[ 1,    2,  3 ],
[ 1,    5, -1 ],
[ 7, 1+1j, -np.pi]
])
v = np.array([ 1, 3, -5])

res = M @ v
print( res )
# [ -8.+0.j  21.+0.j  25.71+3.j ]

res = v @ M
print( res )
# [ -31.+0.j  12.-5.j  15.708+0.j ]



$$\left[ \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 5 & -1 \\ 7 & 1+\mathrm{i} & -\pi \end{array} \right] \circ \left[ \begin{array}{ccc} 1 & 0 & 1 \\ 1 & 0 & -1 \\ -2 & \mathrm{i} & -\pi \end{array} \right] = \left[ \begin{array}{ccc} 1 & 0 & 3 \\ 1 & 0 & 1 \\ -14 & -1+\mathrm{i} & \pi^2 \end{array} \right]$$


A = np.array([
[ 1,    2,  3 ],
[ 1,    5, -1 ],
[ 7, 1+1j, -np.pi]
])

M = np.array([
[ 1,    0,  1 ],
[ 1,    0, -1 ],
[ -2, 1j, -np.pi]
])

res = np.multiply(A, M)

print( res )

# [[  1. +0.j   0.+0.j    3.+0.j  ],
#  [  1. +0.j   0.+0.j    1.-0.j  ],
#  [ -14.+0.j  -1.+1.j  9.87-0.j  ]]


## Scalar products

$$\mathbf{\overline v}_1 = \left[ \begin{array}{c} 1+\mathrm{i} \\ 2 \\ 3 \end{array}\right], \,\, \mathbf{\overline v}_2 = \left[ \begin{array}{c} 1 \\ 2 \\ -6 \end{array}\right]$$

Dot product $\mathrm{\overline v}_1 \cdot \mathrm{\overline v}_2 = 1+\mathrm{i} +2\cdot2+3\cdot(-6) = -13+\mathrm{i}$

Inner product: $\langle \mathrm{\overline v}_1 \mid \mathrm{\overline v}_2 \rangle = \mathrm{\overline v}_1^\dagger \mathrm{\overline v}_2 = \left[ 1-\mathrm{i},\, 2,\, 3 \right] \left[ \begin{array}{c} 1 \\ 2 \\ -6 \end{array}\right] = -13-\mathrm{i}$

In Numpy they are named just the other way around:


v1 = np.array([ 1 +1j , 2, 3 ])

v2 = np.array([ 1, 2, -6 ])

np.inner(v1, v2)
# -13+1j

np.vdot(v1, v2)
# -13-1j


## Outer product

$$\mathbf{\overline v}_1 = \left[ \begin{array}{c} 1+\mathrm{i} \\ 2 \\ 3 \end{array}\right], \,\, \mathbf{\overline v}_2 = \left[ \begin{array}{c} 1 \\ 2 \\ -6 \end{array}\right]$$ $$\mathrm{\overline v}_1 \otimes \mathrm{\overline v}_2 = \mathrm{\overline v}_1 \mathrm{\overline v}_2^\mathrm{T} = \left[ \begin{array}{c} 1+\mathrm{i} \\ 2 \\ 3 \end{array}\right] \left[1, \, 2, \, -6 \right] = \left[ \begin{array}{ccc} 1+\mathrm{i} & 2+2\mathrm{i} & -6-6\mathrm{i} \\ 2 & 4 & -12 \\ 3 & 6 & -18 \end{array}\right]$$


v1 = np.array([ 1 +1j , 2, 3 ])

v2 = np.array([ 1, 2, -6 ])

np.outer(v1, v2)
# [[  1.+1.j,   2.+2.j,  -6.-6.j],
#  [  2.+0.j,   4.+0.j, -12.+0.j],
#  [  3.+0.j,   6.+0.j, -18.+0.j]]


## (Some) norms

$$\mathrm{M}= \left[ \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 5 & -1 \\ 7 & 1+\mathrm{i} & -\pi \end{array} \right] \,\,\,\, \mathrm{\overline v} = \left[ \begin{array}{c} 1 \\ 3 \\ -5+\mathrm{i} \end{array} \right]$$ $$\Vert \mathrm{\overline v} \Vert = \sqrt{ \sum\limits_{i=0}^2 \vert v_i \vert^2 } = \sqrt{1+9+25+1} = \sqrt{36} = 6$$ $$\Vert \mathrm{M} \Vert_\mathrm{f} = \sqrt{ \sum\limits_{i=0}^2\sum\limits_{j=0}^2 \vert M_{ij} \vert ^2 } = \sqrt{92+\pi^2}$$ $$\Vert \mathrm{M} \Vert_\infty = \max\limits_{i\in\{0,1,2\}} \sum\limits_{k=0}^2 \vert M_{ij} \vert = 7 + \sqrt{2} + \pi$$


M = np.array([
[ 1,    2,  3 ],
[ 1,    5, -1 ],
[ 7, 1+1j, -np.pi]
])
v = np.array([ 1, 3, -5+1j])

np.linalg.norm(v)
# 6.0

np.linalg.norm(M, 'fro')
# 10.093 ≈ sqrt(92+π^2)

np.linalg.norm(M, np.inf)
# 11.556 ≈ 7 + sqrt(2) + π


# Solving a system of linear equations

## The equation

$$\left[ \begin{array}{cccc} A_{00} & A_{01} & A_{02} \\ A_{10} & A_{11} & A_{12} \\ A_{20} & A_{21} & A_{22} \end{array}\right] \left[ \begin{array}{c} x_{0} \\ x_{1} \\ x_{2} \end{array}\right] = \left[ \begin{array}{c} b_0 \\ b_1 \\ b_2 \end{array}\right],$$ i.e., $$\mathbf{A} \mathbf{\overline x} = \mathbf{\overline b},$$ where $\mathbf{A}$ and $\mathbf{\overline b}$ are known.

Unique non-trivial solution exists when $\det(A) \neq 0$.

## Finding the solution with SciPy

$$\mathbf{A} \mathbf{\overline x} = \mathbf{\overline b}$$

import numpy as np
import scipy.linalg as sla

A = np.array([
[ 1,    2,  3 ],
[ 1,    5, -1 ],
[ 7,    2,  3 ]
])

b = np.array([ 1, 3, -5 ])

# Check if a good solution exists
assert not np.isclose( sla.det(A), 0)
# -102.0 != 0

x = sla.solve(A,b)
# [-1. 0.824  0.118 ]


## Under the hood of scipy.linalg.solve

### = LU-decomposition

Some $n$-by-$n$ matrices $\mathbf{A}$ can be written in a form $\mathbf{A} = \mathbf{L}\mathbf{U},$ where $\mathbf{L}$ is a lower triangular matrix and $\mathbf{U}$ an upper triangular matrix $$\mathbf{L} = \left[ \begin{array}{ccccc} 1 & 0 & & \cdots & 0 \\ l_{10} & 1 & 0 & & \\ l_{20} & l_{21} & 1 & \ddots & \vdots \\ \vdots & & & \ddots & 0 \\ l_{n-1,0} & & \cdots & & 1 \end{array}\right],\,\,\, \mathbf{U} = \left[ \begin{array}{ccccc} u_{00} & u_{01} & u_{02} & \cdots & u_{0,n-1} \\ 0 & u_{11} & u_{12} & & \\ 0 & 0 & u_{22} & & \vdots \\ \vdots & & \ddots & \ddots & \\ 0 & \cdots & & 0 & u_{n-1,n-1} \end{array}\right]$$ And moreover, we have explicit formulae for calculating $l_{ij}$ and $u_{ij}$ \begin{align} u_{ij} = l_{ij} - \sum\limits_{k=0}^{i-2}l_{ik}u_{kj}, j\geq i\\ l_{ij} = \frac{1}{u_{jj}}\left(a_{ij}-\sum\limits_{k=0}^{j-2} l_{ik}u_{kj}\right), j<i \end{align}

## LU-decomposition with partial pivoting

\begin{align} u_{ij} = l_{ij} - \sum\limits_{k=0}^{i-2}l_{ik}u_{kj}, j\geq i\\ l_{ij} = \frac{1}{u_{jj}}\left(a_{ij}-\sum\limits_{k=0}^{j-2} l_{ik}u_{kj}\right), j<i \end{align} if $u_{jj}=0$ we are in big trouble.

The solution is to permute the rows of $\mathbf A$ in order to avoid divition by small numbers (partial pivoting).

This corresponds to multiplying with a permutation matrix $\mathbf P$. Every $n$-by-$n$ matrix has a LU decomposition after suitable row permutations: \begin{align} \mathbf{P} \mathbf{A} &= \mathbf{L} \mathbf{U} \\ \mathbf{A} &= \mathbf{P}^{-1} \mathbf{L} \mathbf{U} \end{align}

## scipy.linalg.solve

$$\mathbf{A} \mathbf{\overline x} = \mathbf{\overline b}$$

1. First a LU-decomposition with partial pivoting to get $\mathbf{A} = \mathbf{P}^{-1}\mathbf{L}\mathbf{U}$
Now $$\mathbf{P} \mathbf{A} \mathbf{\overline x} = \mathbf{L} \mathbf{U} \mathbf{\overline x} = \mathbf{P} \mathbf{\overline b} \doteq \mathbf{\overline c}$$
2. Solve $\mathbf{\overline y}$ from $\mathbf{L} \mathbf{\overline{y}} = \mathbf{\overline c}$
3. Solve $\mathbf{\overline x}$ from $\mathbf{U}\mathbf{\overline x} = \mathbf{\overline y}$

All steps here have explicit formulae!

## ProTip on inversion of matrices

Sometimes we want to compute the action of matrix inverse on a vector, i.e., $$\mathbf{A}^{-1} \mathbf{\overline v}.$$

This is equivalent to solving $\mathbf{\overline x}$ from $$\mathbf{A} \mathbf{\overline x} = \mathbf{\overline v}$$

We need not explicitely invert the matrix $\mathbf A$ (which is expensive)!

# Eigenvalues of (dense) matrices

## Eigenvalue problem

The eigenvalue equation of an $N$-by-$N$ matrix $\mathbf A$ is $$\mathbf A \mathbf{\overline x}_\lambda = \lambda \mathbf{\overline x}_\lambda,$$ i.e., we try to find those vectors $\mathbf{\overline x}_\lambda$ that only scale by the corresponding eigenvalue $\lambda$ under action of $\mathbf A$.

With pen and paper, we'd solve the eigenvalues $\lambda_0,\,\lambda_1,\,\cdots,\,\lambda_{N-1}$ by finding the roots of $\det\left( \mathbf{A}-\lambda \mathbb{1} \right)$.

On a computer, we do things differently…

## Finding eigenvalues and eigenvectors with SciPy

$$\mathbf A \mathbf{\overline x}_\lambda = \lambda \mathbf{\overline x}_\lambda,$$

import numpy as np
import scipy.linalg as sla

A = np.array([
[ 1,    2,  3 ],
[ 1,    5, -1 ],
[ 7,    2,  3 ]
])

evals, evecs = sla.eig(A)



## QR decomposition

Any square matrix $\mathbf A$ has a decomposition $$\mathbf A = \mathbf Q \mathbf R,$$ where $\mathbf Q$ is a unitary matrix ($\mathbf Q ^\dagger Q = \mathbf 1$) and $\mathbf R$ is an upper triangular matrix.

Several different algorithms to compute this.

### Inside scipy.linalg.eig

1. Calculate the Hessenberg factorization of $\mathbf A$: $$\mathbf A = \mathbf U^\dagger \mathbf H \mathbf U.$$ $\mathbf H$ is similar to $\mathbf A$ and hence they have the same eigenvalues.
2. Calculate the eigenvalues of $\mathbf H$.
1. Calculate the QR decomposition of $\mathbf H_i$ ($\mathbf H_0=\mathbf H$), i.e., $$\mathbf H_i = \mathbf Q_i \mathbf R_i$$
2. $H_{i+1} = \mathbf R_i \mathbf Q_i$ (opposite order)
3. Goto step 2.1 until we have a upper triangular matrix.
3. Finally, we get the relation \begin{align} A &= \mathbf U^\dagger \mathbf R_0^{-1} \mathbf R_1^{-1}\cdots \mathbf R_{N-1}^{-1} \mathbf Q_N \mathbf R_N \mathbf R_{N-1}\cdots \mathbf R_0 \mathbf U \\ & = \mathbf{P}^{-1} \mathbf H_N \mathbf P, \end{align} where $\mathbf H_N$ is (almost) triangular and has the same eigenvalues as $\mathbf A$. Eigenvalues of $\mathbf H_N$ on the diagonal.
4. Eigenvectors of $\mathbf H_N$ e.g. via $$\left(\mathbf H_N-\lambda_i\mathbf{1}\right)\mathbf{\overline x}_i = \mathbf{\overline 0}$$
5. Eigenvectors of $\mathbf A$: $\mathbf P^{-1} \mathbf{\overline x}_i$

# Finding roots of multidimensional nonlinear functions

## Equation to solve

$$\mathbf{\overline f}(\mathbf{\overline x}) = \mathbf{\overline 0},$$ where $\mathbf{\overline f} : \mathbb{R}^N \to \mathbb{R}^N$ is a continuous and differentiable function.

## Broyden's method

#### Broyden's method

1. Pick initial quess $\mathbf{\overline x}_0$
2. Calculate $\mathbf{\overline f}(\mathbf{\overline x}_0)$
3. Calculate first Broyden-Jacobian matrix
$\mathbf B_0 = \left[ \begin{array}{ccc} \frac{\partial f_0}{\partial x_0} & \cdots & \frac{\partial f_{0}}{\partial x_{N-1}} \\ \vdots & & \vdots \\ \frac{\partial f_{N-1}}{\partial x_0} & \cdots & \frac{\partial f_{N-1}}{\partial x_{N-1}} \end{array}\right]$
4. Next quess: $\mathbf{\overline x}_{i+1} = \mathbf{\overline x}_i - \mathbf{B}_i^{-1} \mathbf{\overline f}(\mathbf{\overline x}_i)$
5. Calculate: \begin{align} B_{i+1} &= B_i \\ &+ \frac{\mathbf{\overline f}(\mathbf{\overline x}_{i+1})-\mathbf{\overline f}(\mathbf{\overline x}_i)-\mathbf B_{i}(\mathbf{\overline x}_{i+1}-\mathbf{\overline x}_i)}{\Vert \mathbf{\overline x}_{i+1} - \mathbf{\overline x}_i \Vert^2}\\ &\hphantom{\frac{\mathbf{\overline f}(\mathbf{\overline x}_{i+1})-\mathbf{\overline f}(\mathbf{\overline x}_i)-}{}} \otimes (\mathbf{\overline x}_{i+1} - \mathbf{\overline x}_i) \end{align}
6. If $\Vert \mathbf{\overline f}(\mathbf {\overline x}_{i+1}) \Vert > \epsilon$, i++ and goto step 4

#### 1D Newton's method

1. Pick initial quess $x_0$
2. Calculate $f(x_0)$
3. Calculate $f'(x_0)$
4. Next quess: $x_{i+1} = x_{i}-\frac{1}{f'(x_i)}f(x_i)$
5. Calculate:
$f'(x_{i+1})$ from known (analytic) expression.
6. If $\vert f(x_{i+1}) \vert > \epsilon$, i++ and goto step 4

Other sources:
R. Mäkinen, lecture notes

## Broyden's good method

#### Broyden's method

1. Pick initial quess $\mathbf{\overline x}_0$
2. Calculate $\mathbf{\overline f}(\mathbf{\overline x}_0)$
3. Calculate first Broyden-Jacobian matrix
$\mathbf B_0 = \left[ \begin{array}{ccc} \frac{\partial f_0}{\partial x_0} & \cdots & \frac{\partial f_{0}}{\partial x_{N-1}} \\ \vdots & & \vdots \\ \frac{\partial f_{N-1}}{\partial x_0} & \cdots & \frac{\partial f_{N-1}}{\partial x_{N-1}} \end{array}\right]$
4. Next quess: $\mathbf{\overline x}_{i+1} = \mathbf{\overline x}_i - \mathbf{B}_i^{-1} \mathbf{\overline f}(\mathbf{\overline x}_i)$
5. Calculate: \begin{align} B_{i+1} &= B_i\\ + &\frac{\mathbf{\overline f}(\mathbf{\overline x}_{i+1})-\mathbf{\overline f}(\mathbf{\overline x}_i)-\mathbf B_{i}(\mathbf{\overline x}_{i+1}-\mathbf{\overline x}_i)}{\Vert \mathbf{\overline x}_{i+1} - \mathbf{\overline x}_i \Vert^2}\\ &\hphantom{\frac{\mathbf{\overline f}(\mathbf{\overline x}_{i+1})-\mathbf{\overline f}(\mathbf{\overline x}_i)-}{}} \otimes (\mathbf{\overline x}_{i+1} - \mathbf{\overline x}_i) \end{align}
6. If $\Vert \mathbf{\overline f}(\mathbf {\overline x}_{i+1}) \Vert > \epsilon$, i++ and goto step 4

#### Broyden's good method

1. Pick initial quess $\mathbf{\overline x}_0$
2. Calculate $\mathbf{\overline f}(\mathbf{\overline x}_0)$
3. Calculate inv of the first Broyden-Jacobian matrix
$\mathbf B_0^{-1} = \left[ \begin{array}{ccc} \frac{\partial f_0}{\partial x_0} & \cdots & \frac{\partial f_{0}}{\partial x_{N-1}} \\ \vdots & & \vdots \\ \frac{\partial f_{N-1}}{\partial x_0} & \cdots & \frac{\partial f_{N-1}}{\partial x_{N-1}} \end{array}\right]^{-1}$
4. Next quess: $\mathbf{\overline x}_{i+1} = \mathbf{\overline x}_i - \mathbf{B}_i^{-1} \mathbf{\overline f}(\mathbf{\overline x}_i)$
5. Calculate: \begin{align} \mathbf B_{i+1}^{-1} &= \mathbf B_{i}^{-1} \\ + &\frac{\mathbf{\overline x}_{i+1}-\mathbf{\overline x}_i-\mathbf B_i^{-1}\left[\mathbf{\overline f}(\mathbf {\overline x}_{i+1})-\mathbf{\overline f}\left(\mathbf {\overline x}_{i}\right)\right]}{\left\langle \mathbf{\overline x}_{i+1}-\mathbf{\overline x}_i \mid \mathbf B_i^{-1}\left[\mathbf{\overline f}(\mathbf {\overline x}_{i+1})-\mathbf{\overline f}(\mathbf {\overline x}_{i})\right]\right\rangle } \\ &\hphantom{\frac{\mathbf{\overline f}(\mathbf{\overline x}_{i+1})-\mathbf{\overline f}(\mathbf{\overline x}_i)-}{}} \otimes (\mathbf{\overline x}_{i+1}-\mathbf{\overline x}_i)^\mathrm{T} B_i^{-1} \end{align}
6. If $\Vert \mathbf{\overline f}(\mathbf {\overline x}_{i+1}) \Vert > \epsilon$, i++ and goto step 4

# Matrix functions

• $f(\mathbf A)$, where $\mathbf A \in \mathbb{F}^{n\times n}$ (for us $\mathbb F = \mathbb R$ or $\mathbb F = \mathbb C$)
• Matrix function = yet another matrix, i.e., $f(\mathbf A)\in \mathbb{F}^{n\times n}$
• Can be defined, e.g., by a series expansion:
• Let $f(x) = \sum\limits_{i=0}^\infty \frac{\mathrm d^i f}{\mathrm d x^i}\Big\vert_{x=0}x^i$ with convergence radius $\vert x \vert \leq R$.
• Then $f(\mathbf A)=\sum\limits_{i=0}^\infty \frac{\mathrm d^i f}{\mathrm d x^i}\Big\vert_{x=0}\mathbf A^i$, and it converges whenever $\Vert \mathbf A \Vert \leq R$.
• Example: $\exp(\mathbf A) = 1 + \mathbf A + \frac{1}{2}\mathbf A^2 + \cdots$ which converges for all $\mathbf A$

# Importance of $\exp(\mathbf A)$

• Vectorial differential equations of type \begin{align} \frac{\mathrm d}{\mathrm d t}\mathbf{\overline u}(t) &= \mathbf A \mathrm{\overline u}(t),\\ \mathrm{\overline u}(0)&=\mathbf{\overline u}_0,\,\mathbf{A}\text{ is a matrix} \end{align} have the solution $$\mathbf{\overline u}(t) = \color{LimeGreen}{\exp(t \mathbf A)}\mathbf{\overline u}_0$$
• Differential equations of type \begin{align} \frac{\mathrm d}{\mathrm d t}\mathbf{\overline u}(t) &= \mathbf A \mathrm{\overline u}(t)+\mathbf{\overline f},\\ \mathrm{\overline u}(0)&=\mathbf{\overline u}_0,\,\mathbf{A}\text{ is a matrix} \end{align} have the solution $$\mathbf{\overline u}(t) = \color{LimeGreen}{\exp(t \mathbf A)}\mathbf{\overline u}_0 + \left[\color{LimeGreen}{\exp(t \mathbf A)} -1\right]\mathbf A^{-1} \mathbf{\overline f}$$

• Time-dependent Schrödinger equation \begin{align} \frac{\partial}{\partial t} \psi(\mathbf{\overline r}, t) &= -\mathrm{i}\mathbf H(t) \psi(\mathbf{\overline r}, t) \\ \Rightarrow\psi(\mathbf{\overline r}, t) &= \mathcal{T}\left\{\exp\left[-\mathrm{i} \int\limits_0^t\,\mathrm{d}\tau \mathbf{H}(\tau)\right]\right\} \psi(\mathbf{\overline r}, 0) \end{align}
• Wave-equation \begin{align} \frac{\partial^2}{\partial t^2}u(\mathbf{\overline r},t) &= c^2 \nabla^2 u(\mathbf{\overline r},t) \\ \Rightarrow \frac{\partial}{\partial t}\left[\begin{array}{c} v(\mathbf{\overline r}, t) \\ u(\mathbf{\overline r}, t) \end{array}\right] &= \left[\begin{array}{c} c^2 \nabla^2 u(\mathbf{\overline r, }t) \\ v(\mathbf{\overline r}, t) \end{array}\right] \\ \Rightarrow\left[\begin{array}{c} v(\mathbf{\overline r}, t) \\ u(\mathbf{\overline r}, t) \end{array}\right] &= \exp\left( t \left[\begin{array}{cc} 0 & c^2\nabla^2 \\ 1 & 0 \end{array}\right] \right) \left[\begin{array}{c} v(\mathbf{\overline r}, t=0) \\ u(\mathbf{\overline r}, t=0) \end{array}\right] \end{align}
• (In)homogeneous Maxwell's equations
• $\ldots$

# Krylov subspace

• Consider a vector space $\mathbb F^{n}$, a vector $\mathbf{\overline v}\in\mathbb F^n$, and a matrix $\mathbf{A}\in \mathbb F^{n\times n}$
• By definition $$\mathbf{\overline w} =\exp(\mathbf A)\mathbf{\overline v} = \mathbf{\overline v} + \mathbf A \mathbf{\overline v} + \frac{1}{2}\mathbf A^2 \mathbf{\overline v} + \cdots$$
• Tempting to look for an approximate solution $\mathbf{\overline w}_\mathrm{approx}$ in a subspace $$\mathcal{K}_m(\mathbf A,\mathbf{\overline v}) = \mathrm{span}\left\{\mathbf{\overline v}, \mathbf{A}\mathbf{\overline v}, \ldots, \mathbf{A}^{m-1}\mathbf{\overline v}\right\}$$ called the $m$th order Krylov subspace generated by the matrix $\mathbf A$ and the vector $\mathbf{\overline v}$.
• The core idea: For sparse matrices A@v is cheap.
$\rightarrow$ Can we compute exp(A)@v just by doing matrix-vector products
• $\mathcal K_m(\mathbf A, \mathbf{\overline v})$ is said to be $\mathbf A$-invariant if and only if $\mathbf A^m \mathbf{\overline v} \in \mathcal K_m(\mathbf A, \mathbf{\overline v})$.
$\Rightarrow$ Possibility for exact solution

## Working with $\mathcal{K}_m(\mathbf A, \mathbf{\overline v})$ numerically

• Natural basis, $\mathbf{\overline v}, \mathbf A \mathbf{\overline v}, \ldots, \mathbf A^{m-1}\mathbf{\overline v}$ is numerically terrible (not orthogonal, not normalized)
• Find a better basis:
• Arnoldi method for general square matrices
• Lanczos method for hermitian matrices $\mathbf A^\dagger = \mathbf A$
• Important to be able to grow the size of $\mathcal K_m(\mathbf A, \mathbf{\overline v})$ incrementally.

## Arnoldi method

For finding an orthonormal basis $\mathbf{\overline v}_0,\mathbf{\overline v}_1,\ldots,\mathbf{\overline v}_{m-1}$ of $\mathcal{K}_m(\mathbf A, \mathbf{\overline v})$.
The following is the simplest numerically stable version of the Arnoldi method:

1. $\beta = \Vert\mathbf{\overline v} \Vert$
2. $\mathbf{\overline v}_0 = \frac{\mathbf{\overline v}}{\beta}$
3. For $j=0,\ldots, m-1$ (modified Gram-Schmidt orthonormalization)
1. $\mathbf{\overline p} = \mathbf{A} \mathbf{\overline v}_{j}$
2. For $i=0,\ldots,j$
1. $\mathbf H_{ij} = \langle \mathbf{\overline v}_j \mid \mathbf{\overline p}\rangle$
2. $\mathbf{\overline p} \leftarrow \mathbf{\overline p} - H_{ij} \mathbf{\overline v}_j$
3. $\mathbf H_{j+1, j} = \Vert \mathbf{\overline p} \Vert$
4. If $\mathbf H_{j+1, j} < \epsilon$ (small number),
$\mathcal K_{j+1}$ is $\mathbf{A}$-invariant. Save $m_\mathrm{invariant}=j+1$. Goto step 4.
5. $\mathbf{\overline v}_{j+1} = \frac{\mathbf{\overline p}}{H_{j+1,j}}$
4. If found an $\mathbf{A}$-invariant subspace,
return $\mathbf V_{m_\mathrm{invariant}} = \left[ \mathbf{\overline v}_0, \ldots, \mathbf{\overline v}_{m_\mathrm{invariant}-1}\right]\in \mathbb F^{n\times m}$ and $\mathbf{ H}_{m_\mathrm{invariant}} \in F^{m\times m}$
else:
return $\mathbf V_{m+1} = \left[ \mathbf{\overline v}_0, \ldots, \mathbf{\overline v}_{m}\right]\in \mathbb F^{n\times m+1}$ and $\mathbf{\overline H}_m \in F^{m+1\times m+1}$

## What we get from $m$ Arnoldi iterations?

• An orthonormal basis for the Krylov subspace $\mathcal K_{m}(\mathbf A, \mathbf{\overline v})$: $$\mathbf{\overline v}_0,\mathbf{\overline v}_1,\ldots,\mathbf{\overline v}_{m-1}$$ Note: We can also use the "extra" vector $\mathbf{\overline v}_{m}$ to get more accuracy for our matrix exponential!
• Projection from $\mathbb F^{n}$ to $\mathcal K_{m}(\mathbf A, \mathbf{\overline v})$ is accomplished with the matrix $$\mathbf V_m = \left[ \mathbf{\overline v}_0, \ldots, \mathbf{\overline v}_{m-1}\right]$$
• We also get an upper Hessenberg matrix $$\mathbf{H}_m = \left[\begin{array}{ccccc} h_{0,0} & h_{0,1} & h_{0,2} & \cdots & h_{0,m-1} \\ h_{1,0} & h_{1,1} & h_{1,2} & & h_{1,m-1} \\ 0 & h_{2,1} & h_{2,2} & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & h_{m-2, m-1} \\ 0 & \cdots & 0 & h_{m-1,m-2} & h_{m-1, m-1} \end{array}\right]$$

## What we get from $m$ Arnoldi iterations?

• By inspecting the Arnoldi scheme, we get the relation $$\mathbf H_m = \mathbf V_m^\mathrm{T}\mathbf A\mathbf V_m.$$ $\Rightarrow\,$ $\mathbf H_m$ is projection of $\mathbf A$ onto $\mathcal{K}_m(\mathbf A, \mathbf{\overline v})$
• Another useful relation shows how $\mathbf A$ operates on the basis vectors of $\mathcal K_m(\mathbf A, \mathbf{\overline v})$ $$\mathbf A \mathbf V_m = \mathbf V_m \mathbf H_m + H_{m,m-1}\mathbf{\overline v}_{m}\otimes \mathbf{\overline e}_{m}$$

## Approximating $\exp(A)\mathbf{\overline v}$

The standard approximation is $$\exp(t \mathbf A)\mathbf{\overline v} \approx \beta \mathbf V_m \exp(t \mathbf H_m) \mathbf{\overline e}_0,$$ where $\mathbf{\overline e}_0=\left[1, 0, \ldots, 0\right]^\mathrm{T}\in \mathbb F^m$.

We can use $\mathbf{\overline v}_{m}$ to improve the above approximation with practically no extra cost: $$\exp(t \mathbf A)\mathbf{\overline v} \approx \beta \mathbf V_{m+1} \exp(t \mathbf{ \overline{H}}_{m}) \mathbf{\overline e}_0,$$ where $$\mathbf{\overline H}_{m} = \left[\begin{array}{ccccc} h_{0,0} & h_{0,1} & \cdots & h_{0,m-1} & 0 \\ h_{1,0} & h_{1,1} & \cdots & h_{1,m-1} & 0 \\ 0 & h_{2,1} & \ddots & \vdots & \vdots \\ \vdots & \ddots & \ddots & h_{m-1,m-1}& 0 \\ 0 & \cdots & 0 & h_{m,m-1} & 0 \end{array}\right],$$ and $\mathbf{\overline e}_0 = [1, 0, \ldots,0]^\mathrm{T}\in \mathbb F^{m+1}$.

## Error estimates for $\exp(t \mathbf A)\mathbf{\overline v}$

The improved method has a rough error estimate $$\Vert \exp(t \mathbf A)\mathbf{\overline v} - \beta \mathbf V_{m+1} \exp(t \mathbf{ \overline{H}}_{m}) \mathbf{\overline e}_0\Vert \approx \Vert \mathbf{\overline v}\Vert \left\vert \left\langle \mathbf{\overline e}_{m} \mid \exp\left(t \mathbf{\overline H}_{m}\right) \mathbf{\overline e}_0\right\rangle\right\vert,$$ where $\mathbf{\overline{e}}_m = \left[0,0, \ldots, 1\right]^T \in \mathcal F^{m+1}$.

More error estimates behind the link in the title.

## Approximating $\exp(A)\mathbf{\overline v}$, part 2

$\exp(t \mathbf A)\mathbf{\overline v} \approx \beta \mathbf V_{m+1} \exp(t \mathbf{ \overline{H}}_{m}) \mathbf{\overline e}_0$
How to compute $\exp\left(t \mathbf{ \overline{H}}_{m}\right)$?

• Many methods evailable
• Typically $m \ll n$
$\Rightarrow$ Even if $\mathbf{\overline H}_m$ is dense, we can use ordinary matrix exponential methods.
• Padé approximants are typically used
• The core idea of Krylov subspace methods: Transform the operation of a sparse matrix $\exp(t \mathbf A)\mathbf{\overline v}$ to a small dense operation $\exp(t \mathbf{\overline H}_{m+1})$

# Short note on units

• In computing, use natural units
• What this means depends on the equation we're dealing with
• In Atomic Physics we use (Hartree) atomic units (a.u.):
• Electron mass $m_e=1$
• Elementary charge $q_e=1$
• Reduced Planck's constant $\hbar=1$
• Coulomb constant $1 / (4\pi\epsilon_0) = 1$
• Dimensions should be determined from context:
• Pulse duration $T=1000$ a.u. $\approx 24$ femtoseconds
• Peak laser electric field $\epsilon_\mathrm{max} = 0.1$ a.u. $\approx 0.5 \frac{\text{GV}}{\text{cm}}$
• Nanostructure width $L= 10$ a.u. = 10 bohr radii = $0.53$ nm.
• Often "a.u." is omitted

# Quantum Mechanics in 25 minutes

#### Basics of (single-particle) quantum mechanics

• Particle's status/configuration defined by wavefunction/state \begin{align} \psi(\mathbf{\overline r}) \in &\, W^{2,2}(\mathbb R^3) = \left\{ \psi : \mathbb R^3 \to \mathbb C \text{ for which } \int\limits_{\mathbb R^3}\mathrm{d}\mathbf{\overline r}\, \left\vert \psi(\mathbf{\overline r}) \right\vert^2 < \infty \text{ and }\\ \int\limits_{\mathbb R^3}\mathrm{d}\mathbf{\overline r}\, \left\vert \frac{\partial^s}{\partial r_i^s}\psi(\mathbf{\overline r}) \right\vert^2 < \infty\,\, \forall s\in\{1,2\}\text{ and } i \in\{0,1,2\} \right\} \end{align}
• Inner product between two states $\psi_1(\mathbf{\overline r})$ and $\psi_2(\mathbf{\overline r})$ $$\left\langle \psi_1 \mid \psi_2 \right\rangle = \int\limits_{\mathbb R^3}\,\mathrm{d}\mathbf{\overline r}\, \psi_1^*(\mathbf{\overline r}) \psi_2(\mathbf{\overline r})$$
• Important postulate of mathematical quantum mechanics is that all states that differ merely by a complex factor correspond to the same physical state.
$\Rightarrow$ always use normalized states: $\langle \psi \vert \psi \rangle = 1$
• Probability of finding the particle in a region $\Omega \subset \mathbb R^3$ is $\int\limits_\Omega \,\mathrm{d} \mathbf{\overline r}\, \vert \psi(\mathbf{\overline r}) \vert^2,$ where $\psi$ should be normalized (see the previous bullet point)

#### Measurements in quantum theory

• All measurable quantities ("observables") have corresponding self-adjoint ("hermitian") operators $$\hat{O} : W^{2,2}(\mathbb R^3) \to W^{2,2}(\mathbb R^3) \text{ for which }\\ \langle \psi \mid \hat{O} \phi \rangle = \langle \hat{O} \psi \mid \phi \rangle,\,\forall \psi,\phi \in W^{2,2}(\mathbb R^3)$$
• Possible outcomes of the measurement can be continuous or discrete, but they are always given by the eigenvalues $\lambda$ of the operator in question: $$\hat{O}\psi_\lambda(\mathbf{\overline r}) = \lambda \psi_\lambda(\mathbf{\overline r})$$
• If the system is in a quantum state described by the wavefunction $\psi(\mathbf{\overline r})$, the expectation value of the measurement is $$\langle \psi \mid \hat{O} \psi \rangle = \int\limits_{\mathbb{R}^3}\,\mathrm{d}\mathbf{\overline r}\, \psi^*(\mathbf{\overline r}) \left[\hat{O} \psi\right](\mathbf{\overline r})$$
• If the system is in a quantum state described by the wavefunction $\psi(\mathbf{\overline r})$, the probability density function for getting $\lambda$ out of the measurement is given by $$P(\lambda) = \vert \langle \psi_\lambda \vert \psi \rangle \vert^2$$

#### More material on quantum mechanics

• Example in 1D: A mysterious stranger tells us that the electron we're considering is described by the wavefunction $$\psi(x) = \frac{\sqrt[4]{2}}{\sqrt[4]{9\pi}}\exp\left[-\frac{(x-5)^2}{9} \right]$$ If we measure its position, the probability of finding it within the interval $\left[-10,3\right]$ is $$\int\limits_{-10}^{3}\,\mathrm{d}x\, \psi^*(x) \psi(x) = \cdots \approx 9\,\%$$
• Example in 1D: A mysterious stranger tells us that the electron we're considering is described by the wavefunction $$\psi(x) = \frac{\sqrt[4]{2}}{\sqrt[4]{9\pi}}\exp\left[-\frac{(x-5)^2}{9} \right]$$ If we measure its momentum, the expectation value of the measurement will be $$\left\langle \psi \middle\vert \left(-\mathrm{i}\frac{\mathrm{d}}{\mathrm{d}x}\right) \psi \right\rangle = \int\limits_{-\infty}^{\infty}\,\mathrm{d}x\, \psi^*(x) \left(-\mathrm{i}\frac{\mathrm{d}}{\mathrm{d}x}\right) \psi(x) = \cdots = 0$$

#### Time-evolution in QM

• Time-evolution of a quantum system follows the time-dependent Schrödinger equation (TDSE): $$\mathrm{i}\frac{\partial}{\partial t}\psi(\mathbf{\overline r}, t) = \hat{H}(t) \psi(\mathbf{\overline r}, t)$$ with the initial condition $$\psi(\mathbf{\overline r}, t=0) = \phi(\mathbf{\overline r})$$
• $\hat{H}(t) = -\frac{\nabla^2}{2} + V(\mathbf{\overline r}) + W(\mathbf{\overline r}, \nabla, t)$ is the time-dependent Hamiltonian (energy) operator
• Often the initial state is the lowest eigenstate of the time-independent Hamiltonian $\hat{H}_0=-\frac{\nabla^2}{2} + V(\mathbf{\overline r})$ , $$\hat{H}_0 \phi(\mathbf{\overline r}) = E_\mathrm{gs}\phi(\mathbf{\overline r})$$
• If $\hat{H}(t) \equiv \hat{H}_0$, i.e., it has no time-dependence, the time-evolved wavefunction is given by $$\psi(\mathbf{\overline r}, t) = e^{-\mathrm{i} t \hat{H_0} } \psi(\mathbf{\overline r}, t=0)$$
• If $\hat{H}(t)$ depends on time, the time-evolved wavefunction is given by $$\psi(\mathbf{\overline r}, t) = \lim\limits_{n\to\infty} e^{-\mathrm{i} \Delta t\, \hat{H}[(n-1) \Delta t] } \cdots e^{-\mathrm{i} \Delta t\, \hat{H}(t=0) } \psi(\mathbf{\overline r}, t=0),$$ where $\Delta t = \frac{t}{n}$.
• We write $\psi(\mathbf{\overline r}, t) = \hat{U}(t, t_0) \psi(\mathbf{\overline r}, t_0)$, where $\hat{U}(t, t_0)$ is the time-evolution operator.

#### Light-matter interaction 101

• $\hat{H}(t) = -\frac{\nabla^2}{2} + V(\mathbf{\overline r}) + W(\mathbf{\overline r}, \nabla, t)$ was the time-dependent Hamiltonian (energy) operator
• $-\frac{\nabla^2}{2}$ is the Kinetic energy operator
• $V(\mathbf{\overline r})$ is the static potential (atomic potential, quantum dot potential, ...). What the electron feels from its surroundings.
• $W(\mathbf{\overline r}, \nabla, t)$ is the interaction operator. E.g., interaction with a laser field. Easiest approximation for electron in an atom is the length-gauge dipole approximation $$W = \mathbf{\overline r}\cdot \mathbf{\overline E}(t),$$ where $\mathbf{\overline E}(t)$ is the electric field of the laser.

#### Example: Electron in 1D hydrogen under laser pulse

• The kinetic energy operator: $\hat{T} = -\frac{1}{2} \frac{\partial^2}{\partial x^2}$
• The atomic potential $V(x) = -\frac{1}{\sqrt{x^2+1}}$
• Field interaction: $W = x \epsilon(t)$ where $\epsilon(t)$ is the laser electric field at time $t$
• Total time-dependent Hamiltonian: $$H(t) = -\frac{1}{2} \frac{\partial^2}{\partial x^2} -\frac{1}{\sqrt{x^2+1}} + x \epsilon(t)$$

# Discretizing PDEs with finite differences

• Discretize = continuous $\to$ discrete
• Let's first discretize 1D TISE $$\left[ -\frac{1}{2}\frac{\mathrm{d}^2}{\mathrm{d}x^2} + V(x) \right] \psi(x) = E \psi(x),$$ where $x\in[a,b]\subset \mathbb R$. We often use Dirichlet boundary conditions $$\psi(a) = 0 \text{ and } \psi(b) = 0.$$
• Even if we're interested in unbounded domain, i.e., $x\in\mathbb R$, the easiest method is the truncate to some finite region $[a,b]\subset \mathbb R$ and set up, e.g., Dirichlet boundary conditions.
• Define a uniform coordinate grid $$[a,b] \to \{a = x_0-\Delta x,\, x_0 ,\, x_0+\Delta x,\,\ldots,\, x_0+(N-1)\Delta x,\, x_0+N\Delta x=b\}$$
• Denote the values of the wavefunction at these coordinate by $$\{\psi_a, \psi_0, \psi_1, \ldots, \psi_k = \psi[x_a+(k+1)\Delta x],\ldots,\psi_{N-1},\psi_N=\psi_b\}$$

• Approximate derivatives with finite differences: $$\left.\frac{\mathrm{d}^2\psi}{\mathrm{d}x^2}\right\vert_{x=x_k} \approx \frac{\psi_{k-1}-2\psi_k+\psi_{k+1}}{\Delta x^2}$$
• The TISE must hold for all points in our coordinate grid, and thus $$-\frac{1}{2\Delta x^2}\psi_{k-1}+\left[\frac{1}{\Delta x^2} + V(x_k)\right]\psi_k-\frac{1}{2\Delta x^2}\psi_{k+1} = E \psi_k$$
• Remember that $\psi_a = \psi_b = 0$, and we get $$\left[ \begin{array}{ccccc} \frac{1}{\Delta x^2} + V(x_0) & -\frac{1}{2\Delta x^2} & 0 & \cdots & 0\\ -\frac{1}{2\Delta x^2} &\frac{1}{\Delta x^2} + V(x_1) & -\frac{1}{2\Delta x^2} & \ddots & \vdots\\ 0 & -\frac{1}{2\Delta x^2} &\frac{1}{\Delta x^2}+ V(x_2) & \ddots & 0\\ \vdots & \ddots & \ddots & \ddots & -\frac{1}{2\Delta x^2} \\ 0 & \cdots & 0&-\frac{1}{2\Delta x^2} & \frac{1}{\Delta x^2} + V(x_{N-1}) \end{array} \middle] \middle[\begin{array}{c} \psi_0 \\ \psi_1 \\ \psi_2 \\ \vdots \\ \psi_{N-2}\\ \psi_{N-1} \end{array}\middle] = E \middle[\begin{array}{c} \psi_0 \\ \psi_1 \\ \psi_2 \\ \vdots \\ \psi_{N-2}\\ \psi_{N-1} \end{array}\right]$$

• The previous can be written more clearly as $$\mathbf H \overline \psi = E \overline \psi,$$ where $\overline \psi = [\psi_0,\psi_1,\ldots,\psi_{N-1}]^\mathrm{T}$
• Notice that $\overline \psi$ doesn't contain the wavefunction at the endpoints, $\psi_a$ and $\psi_b$
• Now we have an equation we can solve with linear algebra routines on a computer

## Discretizing time-dependent PDEs with finite differences

• Let's have a look at 1D TDSE $$\frac{\partial}{\partial t}\psi(x,t) = -\mathrm{i} \hat{H}(t) \psi(x, t)$$ with the initial condition $$\psi(x, t=0) = \phi(x)$$
• Even if we're interested in unbounded domain, i.e., $x\in\mathbb R$, the easiest method is the truncate to some finite region $[a,b]\subset \mathbb R$ and set up, e.g., Dirichlet boundary conditions.
• Discretizing coordinate space the same way as with TISE gives us $$\frac{\partial}{\partial t}\overline \psi(t) = - \mathrm{i} \mathbf{H}(t) \overline \psi(t)$$ with the initial condition $$\overline \psi(t=0) = \overline \phi$$ where $\mathbf H$ is the Hamiltonian matrix and $\overline \psi$ is the vector of wavefunction values at the gridpoints.

• The space-discretized TDSE has the formal solution $$\overline \psi(t) = \lim\limits_{n\to\infty} e^{-\mathrm{i} \Delta t\, \mathbf{H}[(n-1) \Delta t] } \cdots e^{-\mathrm{i} \Delta t\, \mathbf{H}(t=0) } \overline \psi(t=0),$$ where $\Delta t= \frac{t}{n}$.
• Time-stepping means that we solve the time-dependent equation by solving first $\overline \psi(\Delta t)$, then $\overline \psi(2\Delta t)$ and so on until we reach the final time.
• Time-stepping can be done, e.g., with the exponential mid-point rule: \begin{align} \overline \psi(t+\Delta t) &= e^{-\mathrm{i}\Delta t \mathbf H(t+\Delta t/2)}\overline \psi(t), \end{align} where the matrix-exponential $\times$ a vector can be computed, e.g., with Krylov subspace methods from last week.

### Modeling unbounded domains

• With wave-like equations such as the TDSE truncating an unbounded domain to a bounded one gives arise to unphysical scattering from the boundaries.

• We can solve this problem, e.g., by:
• Including a purely imaginary potential at the boundaries ("complex absorbing potential"): $$\mathrm{CAP}(x, w) = \left\{\begin{array}{rl} 0, & x\in]a+w,b-w[ \\ -\mathrm{i} C_\mathrm{max}\cos^2\left[\frac{\pi}{2 w} \min(\vert x-a \vert, \vert x-b \vert)\right],& x\in [a,a+w]\cup [b-w,b], \end{array}\right.$$
• Multiplying the wavefunction by a mask function between every time-step: $$\mathrm{M}(x, w) = \left\{\begin{array}{rl} 1, & x\in]a+w,b-w[ \\ 1-\cos^2\left[\frac{\pi}{2 w} \min(\vert x-a \vert, \vert x-b \vert)\right],& x\in [a,a+w]\cup [b-w,b] \end{array}\right.$$

### Modeling unbounded domains

• Other schemes for solving the issue with unbounded domains include, e.g.,:
• Partitioning the space into multiple coupled domains which can be solved with different methods.
• Perfectly matched layers
• Exterior complex scaling (the best for TDSE)

What the solution to TDSE looks like without (left) and with (right) an absorbing boundary:

# What does a TDSE code look like?

#### 1D TDSE solver with finite differences: general structure

$\mathrm{i}\frac{\partial t}{\partial t} \psi(x, t) = \hat{H}(t)\psi(x, t),\,\mathbf{\overline r}\in \mathbb D\subseteq \mathbb R,\,\text{with initial condition } \psi(x, t=0) = \phi(x)$
1. Unbounded or bounded domain $\mathbb D$?
• Truncate if needed
• Decide boundary conditions if truncating
$\Rightarrow$ Finite domain in coordinate space, $\Omega = [a, b]$
$\Rightarrow$ Boundary conditions, e.g., $\psi(a, t) = \psi(b, t)= 0$.
2. Discretize the coordinate space $[a,b] \to \{a = x_0-\Delta x,\, x_0 ,\, x_0+\Delta x,\,\ldots,\, x_0+(N-1)\Delta x,\, x_0+N\Delta x=b\}$
• Remember that (at least for Dirichlet boundary conditions) we don't need to solve for $\psi(a,t)$ and $\psi(b, t)$ so we take the coordinate grid as $\overline x = [x_0, \ldots, x_0+(N-1)\Delta x]^\mathrm{T}$
3. Obtain the time-independent part of the discretized Hamiltonian matrix (which is sparse), $\mathbf{H}_0$
4. Get the initial wavefunction values on the coordinate grid $\overline \psi(t=0) = \phi(\overline{x})$ This is often an eigenstate of $\mathbf H_0$. Remember to normalize the initial state.
1. Add the complex absorbing potential to $\mathbf H_0$ to get $\mathbf H_\mathrm{ti}$, the time-independent part of the Hamiltonian matrix for the time-evolution,
2. while $t$ < $T_\mathrm{max}$
1. Store $\overline \psi(t)$
2. Calculate observables from $\overline \psi(t)$ here or do it in post-processing
3. Calculate the Hamiltonian matrix at time $t+\frac{\Delta t}{2}$, e.g., the laser-electron interaction in dipole approximation $\mathbf H\left(t+\frac{\Delta t}{2}\right) = \mathbf H_\mathrm{ti} + \epsilon\left(t+\frac{\Delta t}{2}\right) \left[\begin{array}{cccc} x_0 & 0 & \cdots & 0 \\ 0 & x_1 & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \cdots & 0 & x_{N-1} \end{array}\right]$
4. Propagate to the next time-step $\overline \psi(t+\Delta t) = e^{-\mathrm{i}\Delta t \mathbf{H}\left(t+\frac{\Delta t}{2}\right)} \overline \psi(t)$
5. $t \leftarrow t+\Delta t$

For more complex systems (2D, 3D, many-particles), you should encapsulate data and methods inside classes to

1. keep the main simulation code easy to read
2. make unit-testing easier

# Supercomputers and superclusters

Most of this material is from csc.fi

### Cluster autopsy: Nodes and connections

• Node = a single "computer" with typically 2-4 CPUs and 0-4 GPUs
• Nodes can discuss with each other via interconnects (high-speed, low-latency links)

### Cluster autopsy: Computing node

• Node has a few CPUs, each with their own RAM

### Cluster autopsy: Processor

• Each CPU has several cores
• Different levels (L1, L2, L3) of memory inside the CPU

# Using Taito-cluster for computing

## Moving files to/from Taito

For a few odd files, you can use

local_pc:~/ $scp myfile.py username@taito.csc.fi:~/  For lots of data, use  local_pc:~/$ rsync -r -v -e ssh myfolder/ username@taito.csc.fi:~/location/


### Logging into Taito


local_pc:~/ $ssh username@taito.csc.fi ... username@taito-login3:~$

When graphics are required:

local_pc:~/ $ssh username@taito.csc.fi -YC ... username@taito-login3:~$ xeyes


Generate ssh-keypair

local_pc:~/ $ssh-keygen  Copy public key to remote  local_pc:~/$ ssh-copy-id -i ~/.ssh/id_rsa.pub username@taito.csc.fi

Create ssh-config

local_pc:~/ $cat .ssh/config Host taito HostName taito.csc.fi User username IdentityFile /home/student/.ssh/id_rsa  Enjoy  local_pc:~/$ ssh taito
...
username@taito-login3:~$ ### Pre-installed software / modules • Lots of preinstalled modules • Can be loaded with module load whatever username@taito-login3$ module avail
---------------------------------------------- /homeappl/home/solanpaa/appl_taito/modules -----------------------------------------------
autotools/12.2015                            compphys/libs/intel/17.0.4/pythonlibs/default
bazel/0.4.4                                  compphys/libs/intel/17.0.4/tclap/1.2.1
compphys/defaults                            compphys/libs/nvcc/8/fftw/3.3.6
compphys/gpu                                 compphys/libs/nvcc/8/gsl/2.3
compphys/gpuflags/fast                       compphys/libs/nvcc/8/gsl/2.4                      (D)
compphys/intelflags/fast                     compphys/libs/nvcc/8/hdf5/1.8.11
compphys/libs/intel/17.0.4/arpack/default    compphys/libs/nvcc/8/libxc/3.0.0
compphys/libs/intel/17.0.4/boost/1.64.0      compphys/libs/nvcc/8/openblas/0.2.20
compphys/libs/intel/17.0.4/fftw/3.3.6-pl2    compphys/libs/nvcc/8/parmetis/default
compphys/libs/intel/17.0.4/gmp/6.1.2         compphys/sw/boost.build/1.64.0
compphys/libs/intel/17.0.4/gsl/2.3           compphys/sw/intel/17.0.4/bill2d/softliebgas_v0.91
compphys/libs/intel/17.0.4/hdf5/1.10.1       compphys/sw/intel/17.0.4/bill2d/8a042c3           (D)
compphys/libs/intel/17.0.4/libxc/3.0.0       compphys/sw/intel/17.0.4/bill2d/47b0022
compphys/libs/intel/17.0.4/mpc/1.0.3         compphys/sw/intel/17.0.4/itp2d/1.0.0.2
compphys/libs/intel/17.0.4/mpfr/3.1.5        compphys/sw/intel/17.0.4/octopus/4eb3136
compphys/libs/intel/17.0.4/nfft/3.3.2        compphys/sw/nvcc/8/octopus/4eb3136
compphys/libs/intel/17.0.4/nlopt/2.4.2       compphys/sw/nvcc/8/octopus/6360360                (D)
compphys/libs/intel/17.0.4/parmetis/4.0.3    itp2d/2e7ab21
compphys/libs/intel/17.0.4/pfft/5b0be20      itp2d/41938ca                                     (D)
compphys/libs/intel/17.0.4/pnfft/e34d373     python/2.7.13

------------------------------------------- /appl/modulefiles/MPI/intel/16.0.0/intelmpi/5.1.1 -------------------------------------------
amber/16         elmer/permafrost        esmf/7_1_0_beta      gromacs/5.1.1-mic        mumps/4.10.0       roms/roms_fisoc
cpmd/4.1         elmer/release           fftw/3.3.4           gromacs/5.1.2-mic (D)    netcdf4/4.3.3.1    roms/roms       (D)
elmer/a7b00af    elmer/release81         fisoc/fisoc          hdf5-par/1.8.15          nwchem/6.8
elmer/fisoc      elmer/release82         flexpart-wrf/3.3     hypre/2.9.0b             openifs/38r1v04
elmer/latest     elmer/release83  (D)    gromacs/5.0.7-mic    molpro/2015.1            parmetis/3.2


username@taito-login3:~$module spider python ----------------------------------------------------------------------------------------------------------------------------------------- python: ----------------------------------------------------------------------------------------------------------------------------------------- Versions: python/serial-2.7.13 python/2.7.3 python/2.7.6 python/2.7.9 python/2.7.10-intel python/2.7.10 python/2.7.13 python/3.4.0 python/3.4.1 python/3.4.5 python/3.5.3 Other possible modules matches: biopython-env compphys/libs/intel/17.0.4/pythonlibs geopython python-env ----------------------------------------------------------------------------------------------------------------------------------------- To find other possible module matches do: module -r spider '.*python.*' ----------------------------------------------------------------------------------------------------------------------------------------- To find detailed information about python please enter the full name. For example:$ module spider python/2.7.10-intel


username@taito-login3:~$module list Currently Loaded Modules: 1) intel/16.0.0 2) mkl/11.3.0 3) intelmpi/5.1.1 4) StdEnv username@taito-login3:~$ python3 --version

username@taito-login3:~$module load python-env/intelpython3.5 Loading application Intel Distribution for Python 2017 update 1 username@taito-login3:~$ python3 --version
Python 3.5.2 :: Intel Corporation

username@taito-login3:~$python3 -c "import GPy" Traceback (most recent call last): File "<string>", line 1, in <module> ImportError: No module named 'GPy' username@taito-login3:~$ pip install --user GPy
...
8776
username@taito-login3:~$python3 -c "import GPy" username@taito-login3:~$


### Interactive simulations

DO NOT RUN SIMULATIONS ON THE LOGIN NODES!

local_pc:~/ $ssh taito -YC ... username@taito-login3:~$ ssh taito-shell -YC
...
username@c307:~$module load python-env/intelpython3.5 ... username@c307:~$ python3 run_my_super_heavy_interactive_simulation.py


### Non-interactive simulations

1. Write a job file (e.g., job.sub) describing
• Resources (CPUs, RAM, runtime)
• Commands to execute
2. Submit job to scheduler for queueing
 $sbatch job.sub 3. Monitor your job via output files or $ squeue -u whoami

### Batch job files

#!/bin/bash -l
#SBATCH -J YOUR_JOBNAME
#SBATCH -e OUTPUTFILE_NAME
#SBATCH -o OUTPUTFILE_NAME
#SBATCH --mem-per-cpu=4000
#SBATCH -t 12:00:00
#SBATCH -N 1
#SBATCH -n 1
#SBATCH -p serial
ml python-env/intelpython3.5

srun python3 run_my_script.py


### How much does it cost?

• CSC uses billing units to measure service usage
• 1 billing unit = 0.01 €
• 1 BU = 0.5 hours with one CPU core on Taito and Sisu
• A short 1 hour quantum simulation with 40 full nodes ≈ 4000 BU = 40 €
• A week long simulation with 40 full nodes ≈ 135 €
• Most of this is free for Finnish academic institutions (i.e., someone else pays for it)
• Some other providers:
• Amazon AWS
• Penguin Computing POD

### FENICS Example: Eigenmodes

Visualized with ParaView (note: the ParaView version in Ubuntu package repository is broken, install the one from their webpage if you want to use ParaView within the VM)

## Other FEM software suites

• Comsol (TUT has a student license available)
• Ansys (TUT has some licenses?)
• Matlab has a rudimentary PDE toolbox
• Mathematica has a rudimentary PDE toolbox

# Ordinary Differential Equations

## ODE

• ODE = an equation for a single-variable function containing first and possibly higher-order derivatives of the function $y'(x) + 3 y''(x) -1/y'''(x) = \sin^3 (x)$ One has to add suitable set of
• initial conditions [ $y(0) = 1,\,y'(0) = 4, y''(0) = -2$ ], or
• boundary conditions [ $y(-1) = -4, \, y(4) = 4$ ]
• The necessary/sufficient amount of boundary/initial conditions in order to have a unique (and existing) solution depends on the ODE in question
• ODEs can be used to model multiple phenomena:
• Classical mechanics (particles, vehicles, rigid body motion, ...)
• Nuclear decay
• ...

### Technique: Reduction of order

#### Motivation:

Most ODE solvers can solve only 1st order ODEs
$\Rightarrow$ Higher-order ODEs have to be transformed to lower order ODEs

#### Example: $y''(x) + 4 \cos[ y'(x) ] = x^3-4$

1. Define a new function $v(x) \equiv y'(x)$
2. Now we get two coupled first order ODEs: $\left\{ \begin{array}{rcl} y'(x) & = & v(x) \\ v'(x) & = & -4\cos[v(x)] + x^3 - 4 \end{array}\right.$
3. Profit.

## Implicit and explicit ODE solvers

Our ODE is of type $\frac{\mathrm{d}\mathbf{\overline y}}{\mathrm{d}x} = \mathbf{\overline f}(x,\,\mathbf{\overline y})$ with the initial condition $\mathbf{\overline y}(x=0) = \mathbf{\overline y}_0$.

#### Explicit methods

• Explicit methods are of type $\mathbf{\overline y}(x + \Delta x) \approx G[x, \mathbf{\overline y}(x), \mathbf{\overline f}]$
• $\Delta x$ is the step size. It determines the accuracy of our numerical solution.
• For example Euler's method: $$\begin{split} \mathbf{\overline f}[x,\,\mathbf{\overline y}(x)] &= \frac{\mathrm{d}\mathbf{\overline y}}{\mathrm{d}x} = \frac{\mathbf{\overline y}(x+\Delta x) - \mathbf{\overline y}(x)}{\Delta x} + \mathcal{O}(\Delta x) \\ \Rightarrow\mathbf{\overline y}(x+\Delta x) &= \mathbf{\overline y}(x) + \Delta x \cdot \mathbf{\overline f}(x,\,\mathbf{\overline y}) + \mathcal{O}(\Delta x^2) \end{split}$$ which is a first order method (whose error is of second order)

#### Implicit methods

• Implicit methods are of type $F[x, \mathbf{\overline y}(x), \mathbf{\overline y}(x+\Delta x),\mathbf{\overline f}] = 0,$ from which we cannot derive an explicit formula for $\frac{\mathrm{d}\mathbf{\overline y}}{\mathrm{d}x}$
• For example implicit Euler's method: $$\begin{split} \mathbf{\overline f}(x,\,\mathbf{\overline y(x+\Delta x)}) &= \left.\frac{\mathrm{d}\mathbf{\overline y}}{\mathrm{d}x}\right\vert_{x+\Delta x} = \frac{\mathbf{\overline y}(x+\Delta x) - \mathbf{\overline y}(x)}{\Delta x} + \mathcal{O}(\Delta x) \\ \Rightarrow\mathbf{\overline y}(x+\Delta x) - \mathbf{\overline y}(x) &- \Delta x \cdot \mathbf{\overline f}[x,\,\mathbf{\overline y}(x+\Delta x)] = 0 + \mathcal{O}(\Delta x^2) \end{split}$$
• Implicit methods require extra work for solving the (non)linear equation for $\mathbf{\overline y}(x+\Delta x)$
• But implicit methods are numerically more stable than explicit

## Explicit Runge-Kutta methods

• Explicit RK methods of order $s$ are numerical solvers for $y'(x) = f[x, y(x)]$ of type $y(x+\Delta x) = y(t) + \Delta x \sum\limits_{i=1}^s b_i k_i + \mathcal{O}(\Delta x^{s+1}),$ where $b_i$ are real coefficients with $\sum\limits_{i=1}^s b_i=1$, $k_i = f\left[x + c_i \Delta x, y(x)+\Delta x \sum\limits_{j=1}^{i-1}a_{ij}k_j\right],\;i=1,\ldots,s$ are estimates for the derivative, $c_i$ real coefficients with $\sum\limits_{i=1}^{s}c_i = 1$, and $a_{ij}$ some real coefficients.

### Example: RK2

1. Taylor expansion: $y(x+\Delta x) = y(x) + y'(x) \Delta x+ \frac{1}{2}y''(x)\Delta x^2 + \mathcal{O}(\Delta x^3)$
Remember also
1. $y'(x) = f[x, y(x)] \doteq f$
2. $y''(x) = \frac{\mathrm{d}}{\mathrm{d}x}f[x, y(x)] = \partial_x f + f \partial_y f$
$\Rightarrow y(x+\Delta x) = y(x) + \Delta x f+\frac{\Delta x^2}{2}\left[ \partial_x f + f \partial_y f\right] + \mathcal{O}(\Delta x^3)$
2. Our RK expansion was $y(x+\Delta x) = y(x) + \Delta x b_1 k_1 + \Delta x b_2 k_2 + \mathcal{O}(\Delta x^3).$ $\Rightarrow$ expand $k_1$ and $k_2$ around $\Delta x = 0$:
1. $k_1 = f[x, y(x)]$
2. $k_2 = f[x+ c_2 \Delta x, y(x) + \Delta x a_{21} k_1]$
$= f[x, y(x)] + c_2\Delta x \partial_x f[x, y(x)] + a_{21}\Delta x f \partial_y f + \mathcal{O}(\Delta x^2)$
3. $\Rightarrow y(x+\Delta x) = y(x)+\Delta x (b_1+b_2) f +\Delta x^2 b_2 (c_2 \partial_x f + a_{21}f \partial_y f) + \mathcal{O}(\Delta x^3)$

### Example: RK2

1. Require that the Taylor expansion and RK expansion are equivalent up to $\mathcal{O}(\Delta x^3)$ to get $\left\{ \begin{array}{l} b_1+b_2 = 1 \\ b_2 c_2 = \frac{1}{2} \\ b_2 a_{21} = \frac{1}{2} \end{array}\right.$

No unique solution!
Possible solution, e.g., $(b_1,b_2,c_2,a_{21}) = \left(\frac{1}{2}, \frac{1}{2}, 1, 1\right)$

#### How to compute a single RK step?

1. $k_1 = f[x, y(x)]$
2. $k_2 = f[x+ c_2 \Delta x, y(x) + \Delta x a_{21} k_1]$
3. $y(x+\Delta x) = y(x) + \Delta x (b_1 k_1 + b_2 k_2)$

## Classical (Hamiltonian) mechanics and conservation laws

• Many systems in classical mechanics have conserved quantities, e.g.,
total energy, total angular momentum, …
• Most ODE solvers do not conserve these quantities
• E.g., standard Runge-Kutta methods don't conserve the total angular momentum of Sun-Earth-Moon system:

### Some problems if violating conservation laws

• System starts to spin without external forces
• Total energy increases (system heats up)
• Momentum is not conserved (weird collision physics)

### Remedy: Symplectic integrators

• Integrator = method for solving ODEs = scheme = algorithm
• Symplectic integrator = integrator that conserves certain geometric structures of the phase space (keeps the symplectic structure intact)
• More details in the course Advanced Classical Mechanics or the book V. I. Arnold: Mathematical Methods of Classical Mechanics
• Example: Velocity Verlet
1. $\mathbf{\overline r}(t+\Delta t) = \mathbf{\overline r}(t) + \mathbf{\overline v}(t) \Delta t + \frac{1}{2}\mathbf{\overline a}(t) \Delta t^2$
2. $\mathbf{\overline a}(t+\Delta t) = \frac{1}{m}\mathbf{\overline F}[\mathbf{\overline r}(t+\Delta t) ]$
3. $\mathbf{\overline v}(t+\Delta t) = \mathbf{\overline v}(t) + \frac{\Delta t}{2}[\mathbf{\overline a}(t+\Delta t) + \mathbf{\overline a}(t) ]$
• Note: you have to compute first step 1 for all particles, then steps 2 and 3

## Molecular dynamics (MD)

• Modeling the dynamics of atoms and molecules using classical mechanics.
• Interactomic interactions due to quantum mechanics and classical electron-electron interactions taken into account e.g.,
• via empirical potentials, e.g., the Lennard-Jones potential $U(r) = 4 \epsilon \left[ \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6 \right]$
• by coupling QM electrons to (semi)classical nuclei (Ehrenfest dynamics, coupled quantum-classical methods)

### Modeling gases with MD

• Computers can not handle infinitely many particles
• Take a big number of particles and apply periodic boundary conditions (in a big box):
• unit-cell: $[0,L]^3$
• When a particle goes outside the box, it's injected back in from the opposite side (periodicity).
• Interparticle interactions must be calculated ALSO OVER PERIODIC BOUNDARIES
• Multiple techniques (Direct summation, Ewald summation, ...)
• Box-size determined by preset average density $\rho$ and number of particles $N$ we're willing to model: $N \left\langle m_{atom}\right\rangle = \rho L^3$
• Initial positions uniformly random (unless we start from a lattice configuration)
• Initial speeds drawn from the Maxwell-Boltzmann distribution. Directions uniformly random.
• Even better: include friction (damping) forces and heat bath (random forces)
$\Rightarrow$ Langevin equation and stochastic differential equations

## Beware of chaos

• Chaos = ultrasensitivity to initial conditions in a deterministic system
• Inherent property of many systems (and ODEs), e.g., in chemistry, biology, physics, and mathematics
• Arises often in nonlinear dynamics

### Example: Chua's circuit

• Described by the coupled ODE $$\begin{split} \frac{\mathrm{d}x}{\mathrm{d}t} &= \alpha [y-x-f(x)] \\ \frac{\mathrm{d}y}{\mathrm{d}t} &= x-y+z \\ \frac{\mathrm{d}z}{\mathrm{d}t} &= - \beta y, \end{split}$$ where $x,y,z$ essentially describe the voltages of the capacitors $C_1,C_2$ and the current over $L$. $f(x) = m_1 x + (m_0 - m_1) ( \vert x+1\vert - \vert x-1 \vert )$ describes the nonlinear resistance of the resistor $R_{nl}$.

### Chua's circuit continues

• For some parameters $\alpha, \beta, m_0,$ and $m_1$ the system is chaotic:
• Take two trajectories that are initially very close to each other
$\rightarrow$ their distance will grow exponentially

### Chua's circuit continues

• For some parameters $\alpha, \beta, m_0,$ and $m_1$ the system is chaotic:
• Take two trajectories that are initially very close to each other
$\rightarrow$ their distance will grow exponentially

## Initial value problems with PDEs (recap)

Ingredients

1. Initial value $u(\mathbf{\overline r}, t=0) = u_0(\mathbf {\overline r})$
2. Time-dependent PDE, e.g., $\partial_t u(\mathbf{\overline r}, t) = -\alpha \nabla^2 u(\mathbf{\overline r}, t) \,\forall \mathbf{\overline r} \in \Omega$
3. With some appropriate boundary conditions $u(\mathbf{\overline r}, t) = D(\mathbf{\overline r})\,\forall \mathbf{\overline r} \in \partial \Omega$

### Example: Waves in a box

• Initial condition: High wave in one of the corners.
• Boundary condition: Hard walls (Neumann boundary condition)
• Governing equation: $\partial_t^2 u(\mathbf {\overline r}, t) = -c^2 \nabla^2 u(\mathbf {\overline r}, t)$

## How to solve initial value PDEs?

$\partial_t^2 u + \partial_t u + 0.1\, u\, \partial_x^2 u = e^{-\frac{x^2}{20} \sin\left(2\pi \frac{t}{7}\right)},\; x\in[-20, 20]\\ u(x, t=0) = \frac{1}{1+x^2},\; u(x+40) = u(x)$
1. Reduce the order of time derivatives to 1 (see also ODEs: reduction of order): $\partial_t u \doteq w\\ \Rightarrow \left\{ \begin{array}{rcl} \partial_t u & = & w \\ \partial_t w & = & -w - 0.1\,u\, \partial_x^2 u + e^{-\frac{x^2}{20} \sin\left(2\pi \frac{t}{7}\right)} \end{array}\right.$
2. Weak formulation for the spatial part. We take test functions $v_1,\, v_2 \in H^1(\Omega)$ and trial functions $u,\, w \in \{ \nu \in H^1(\Omega): \nu(x) = \nu(x+40)\}$.
1. $$\left\{ \begin{split} \partial_t \int v_1 u \,\mathrm{d}x &= \int v_1 w \,\mathrm{d}x\\ \partial_t \int v_2 w \,\mathrm{d}x &= - \int v_2 w \,\mathrm{d}x + 0.1 \int \partial_x (v_2 u) \partial_x u \,\mathrm{d}x + \int v_2 e^{-\frac{x^2}{20} \sin\left(2\pi \frac{t}{7}\right)}\,\mathrm{d}x \end{split}\right.$$ Add them up!
1. Implicit Euler: $\partial_t u = f[ t, u(t+\Delta t) ] \Rightarrow \frac{u(t+\Delta t) - u(t)}{\Delta t} = f[t+\Delta t, u(t+\Delta t)]$ $$\begin{split} &\int v_1 \frac{u(t+\Delta t)-u(t)}{\Delta t} \,\mathrm{d}x + \int v_2 \frac{w(t+\Delta t) - w(t)}{\Delta t} \,\mathrm{d}x \\ =& \int v_1 w(t+\Delta t) \,\mathrm{d}x - \int v_2 w(t+\Delta t) \,\mathrm{d}x \\ &+ 0.1 \int \partial_x [ v_2 u(t+\Delta t)] \partial_x u(t + \Delta t) \,\mathrm{d}x + \int v_2 e^{-\frac{x^2}{20} \sin\left(2\pi\frac{t+\Delta t}{7}\right)}\,\mathrm{d}x \end{split}$$

Note: Implicit Euler becomes unstable very early!

1. Implicit Euler
2. Pretty much every other solver:
$$\begin{split} \partial_t \underbrace{\int (v_1 u + v_2 w) \,\mathrm{d}x}_{\doteq L \mathbf{\overline c}(t)} &= \underbrace{\int v_1 w \,\mathrm{d}x - \int v_2 w \,\mathrm{d}x + 0.1 \int \partial_x (v_2 u) \partial_x w \,\mathrm{d}x + \int v_2 e^{-\frac{x^2}{20} \sin\left(2\pi \frac{t}{7}\right)}\,\mathrm{d}x}_{\doteq \mathbf{\overline b}[t, \mathbf{\overline c}]} \end{split}$$
$$\begin{split} \Leftrightarrow& \partial_t L \mathbf{\overline c}(t) = \mathbf{\overline b}[t, \mathbf{\overline c}(t)] \\ \Leftrightarrow& \partial_t \mathbf{\overline c}(t) = L^{-1} \mathbf{\overline b}[t, \mathbf{\overline c}(t)] \end{split}$$ This is just a (large) system of coupled ODEs! We already solved these last week, and can use any solver we wish.

Note: Other solvers might be more stable and faster, but they require more manual work.

## Stiff ODEs

• Sometimes regular ODE solvers (including RK) have to use extremely small time-steps to get accurate results even if the solution behaves nicely.
= stiff ode
• Stiff ODEs often arise in initial value FEM
• Stiff ODEs require specialized methods for efficiency
• Multiple different methods exist

### Backward differentiation formula

• A class of popular numerical integrators for stiff ODEs
• Basic idea is to use multiple previous time-steps to estimate the next value, e.g., $y(t+\Delta t) -\frac{4}{3}y(t) + \frac{1}{3}y(t-\Delta t) = \frac{2}{3}\Delta t f[t+\Delta t, y(t+\Delta t)]$
• Superior stability
• Initialization (first $n-1$ steps of $n$th order method) with, e.g., high-order RK with extra small $\Delta t$

## Solving difficult linear systems with preconditioning

• In initial value problems of FEM (and many other applications) we may encounter extra difficult linear systems of equations $L \mathbf{\overline x} = \mathbf{\overline b}$
• Condition number of a linear map $L$ (and $L^{-1}$) $\kappa (L) = \Vert L \Vert \Vert L^{-1} \Vert$ tells us roughly speaking "how much the result changes when the input changes slightly".
• If $\kappa(L)$ is big, then even small numerical errors will increase the inaccuracy of our solution $\mathbf{\overline x}$.
• Preconditioning = find equivalent problem with smaller condition number $L P^{-1} \mathbf{\overline y} = \mathbf{\overline b} \\ P \mathbf{\overline x} = \mathbf{\overline y},$ where $\kappa(LP^{-1}) < \kappa(L)$ for our preconditioner (matrix) $P$.
• Good preconditioner depends on the problem at hand.

# Time-series analysis

## Basic terminology

Stationary time series
Mean and variance don't depend on time
Trend
General non-oscillatory increase/decrease in data
Cyclic/seasonal component
Low-frequency oscillatory part of the time-series
Irregular component
Signal - trend - cyclic component
Noise
Unphysical irregular component of the data that arises e.g., from measurement errors

## Detrending

• Removing trend from the time series
• Important part in stationarizing data
• Parametrized detrending = fit a function to data
• Adaptive methods such as Empirical Mode Decomposition

### Detrending with a functional form

1. Visualize data
2. Quess the functional form (if no preknowledge available), e.g., a power law $f(t) = a t^b+c$
3. Fit the function to data
4. Remove the fitted function from the original data
1. #### Typical functions

• Linear
• Cubic
• Power law
• Exponential
• Piecewise linear

## Empirical mode decomposition (EMD)

• Method to decompose the signal into multiple components called intrinsic mode functions (IMF)
• IMFs are oscillatory functions which represent the high- and low-frequency oscillatory components of the signal
• IMF do not have a single absolute frequency: their frequencies adapt locally to the signal

## EMD

1. Extract an intrinsic mode function (IMF) via a sifting process:
1. Find local maxima (minima) of the signal and connect the maxima (minima) with an interpolating function, e.g., a cubic spline
1. Extract an intrinsic mode function (IMF) via a sifting process:
1. Substract the average of the splines from the signal.
1. Extract an intrinsic mode function (IMF) via a sifting process:
1. Repeat steps 1.I and 1.II until sufficient convergence
1. Extract an intrinsic mode function (IMF) via a sifting process:
1. Repeat steps 1.I and 1.II until sufficient convergence
2. Substract the latest IMF from the original signal and repeat step 1 for the new signal without the highest frequency IMF.
3. Stops when we have just a monotonic function left.

## Ensemble Empirical Mode Decomposition (EEMD)

• Improved version of EMD, fixes certain problems
• Main difference:
1. Create ensemble of noisy signals
2. IMFs are the ensemble average of the IMFs of noisy signals

## Spectral analysis

• Analysis based on frequencies in the data
• Detrend the data first!

(Power) spectral density:
$S(f) = \left\vert \int\limits_{-\infty}^{\infty} s(t) e^{-2\pi \mathrm{i} f t}\,\mathrm{d}t \right\vert^2$ or $S(f) = \left\vert \mathrm{FFT}[s(t), t\to f] \right\vert^2$
Short Time Fourier Tranfrom (STFT):
$s(\tau, f) = \int\limits_{-\infty}^{\infty} s(t) w(t-\tau) e^{-2\pi \mathrm{i} f t}\,\mathrm{d}t$ or $s(\tau, f) = \mathrm{FFT}[s(t)w(t-\tau), t\to f]$

### Spectral density

• Describes how much the entire signal contains certain frequencies
$S(f) = \left\vert \mathrm{FFT}[s(t), t\to f] \right\vert^2$

### Short Time Fourier Transform

$s(\tau, f) = \int\limits_{-\infty}^{\infty} s(t) w(t-\tau) e^{-2\pi \mathrm{i} f t}\,\mathrm{d}t$
• Frequency content at each instant of time
• Have to balance between resolution in time and resolution in frequency

## Correlations in time-series

Correlation
Statistical relationship between two signals
Cross-correlation
Similarity between two time series for different time delays, i.e., lags
Autocorrelation
Cross-correlation of signal with itself. Reveals periodic structures using time-domain techniques (cf. spectral density)

### Cross-correlation

• Detrend the signals first!
• $c(\tau) = \int\limits_{-\infty}^{\infty} s_1(t) s_2(t+\tau)\,\mathrm{d}t$
• Useful, e.g., in signal matching and finding time-delays between signals

### Autocorrelation

• Detrend the signals first!
• $c(\tau) = \int\limits_{-\infty}^{\infty} s(t) s(t+\tau)\,\mathrm{d}t$
• Finds repeating patterns in the signal

## Filtering

$s_\mathrm{filtered}(t) = \mathrm{IFFT}\{ G(f)\, \mathrm{FFT}[s(t), t \to f],\, f \to t\}$
• Filtering = removing certain unwanted components from the signal
• Low-pass filter = allows frequencies lower than $f_c$ to survive. E.g., $G(f) = \Theta(f_c - f)$
• High-pass filter = allows frequencies higher than $f_c$ to survice. E.g., $G(f) = \Theta(f - f_c)$
• Band-pass filter = allows frequencies in a range $f_a,\ldots f_b$ to survive. E.g., $G(f) = \Theta(f - f_a)\Theta(f_b - f)$
• Use Spectral Density for inspecting the data in frequency domain
• Detrend first!

## Low-pass filter

• For removal of, e.g., high-frequency noise

## High-pass filter

• E.g., for inspecting suspected noise

### Lag plot (Poincaré plot)

• Visualization of $\Big( s(t), s(t+\Delta t) \Big)$
• Can reveal non-randomness of data, sometimes even the exact origin

# Testing, debugging, profiling, and optimizing numerical codes

Once your code runs and provides somewhat reasonable-looking results, it's often difficult to gain confidence and evidence that your software does what it's supposed to do. A few tips that will help:

Unit-testing
Provides reasonable trust in the individual modules, classes, and functions (is data read correctly, do intermediate steps work correctly).
Analytically solvable test-cases
Demonstrates that the code provides exact results for analytically solvable cases.
Known asymptotic behavior
E.g., 3D poisson equation should yield a solution with asymptotic behaviour $\sim 1/r$ far from electric charges.
Comparison with other software
Compare results to software that (1) implement the same numerical method and (2) are designed to solve similar problems.

Convergence properties
If your numerical method should approach the exact result as $\mathcal{O}(h^3)$, you should test and demonstrate this.
Scaling of computational cost
Same idea as in the previous, but now you should know, e.g., the CPU and/or memory requirements of your numerical method.

## Debugging

• Use a debugger, it let's you run your code interactively:
• import pdb; pdb.set_trace()
• python3 -m pdb your_script.py
• import your_module
import pdb
pdb.run('your_module.function()')
• Execute a single line including all functions: n
• Execute a single line but step into first function: s
• Fast-forward to line nbr 5: j 5
• Show arguments of current function: a
• Setting a breakpoint: b [filename:]lineno | function
• Continue until next breakpoint: c
• Show currently executing file, function, and line number: w
• Show a few lines of code: l or ll
• Print the value of an expression: p expr

## Narrowing down the location of the bug

• Need to be able to reason expected state of the program for each line
• Use the bisection method to quickly find the problematic section of the code.

## Finding solutions using error messages


$python3 problem1.py File "problem1.py", line 9 ^ SyntaxError: invalid character in identifier  1. Look at the line of code and a few lines before it. Do you see any obvious mistakes? 2. Copy/paste the error message verbatim to Google (the last line). 3. See if there's a stackoverflow thread or a blog post about it. 1. If error arises from a function, check known issues issues on the corresponding project's webpage for possible workarounds. ## Asking for help (and reporting bugs) #### Things to check first • You have debugged the code and narrowed down the issue as much as you can. • You have already used your google-fu skills • You have read the relevant documentation thoroughly (RTFM) • If your problem is related to a more complex setup, extract a minimal working example (of your issue). • Double check FAQs and issue trackers. #### Help request/bug report • Always use the requested method of communication (issue tracker, email, ...) • Explain your issue in a consise manner and include a working example of code that produces the issue. If you try to achieve something, write it as pseudocode to make your meaning clear. • Tell the exact version of your OS, Python interpreter, package etc. ## Benchmarking = assessing software performance • Don't run inside VM • Close all other applications, including AV • Measure the runtime multiple times! • Always compare benchmarking results fromthe exact same machine! $ time qdyn_laser

real	0m33.753s
user	0m31.424s
sys     0m1.100s

• user + sys tells how much CPU time your code has used.

## Profiling

= checking where your code spends its time/memory


$pip3 install pyprof2calltree$ python3 -m cProfile -o output.cprof scripts/qdyn_laser
$pyprof2calltree -k -i qdyn.cprof  ## Coded optimization in HPC = making your software run faster and use less memory • Profile first with a realistic simulation to see where in code time is spent. • Optimize only bottlenecks • Things that take a lot of time: Reading/writing files Do this in a new thread while the main thread(s) can continue with simulations. Avoid data-races. Terminal output Minimize, this is surprisingly slow Memory allocation Allocate large chuncks of memory at once instead of little pieces many times:  a=[]; b=[] for i in range(N): a.append( foo ) b.append( bar ) ab = np.zeros((N,2)) for i in range(N): ab[i,0] = foo ab[i,1] = bar  Fragmented memory  particles = [ ( np.zeros(3), np.zeros(3) ) for i in range(20) ] particles = np.zeros((20, 6))  Transferring data between processes Important in MPI-parallelization ## Optimization = finding extrema of a scalar valued function$f(\mathbf{\overline x})$• Necessary condition for an extrema:$\nabla f(\mathbf{\overline x}) = 0$• Sufficient conditions:$\nabla f(\mathbf{\overline x}) = 0$and Hessian$H(\mathbf{\overline x})$is positive/negative definite. ## Real optimization examples The optimization target (function) may be difficult to evaluate. Closed form is usually not available. Parameter space may be infinite dimensional. • Revenue from a marketing campain. The adjustable variables are, e.g., advertisment media, ad visibility, target group, campain duration, coupons, offers... • Grade from FYS-4096. Adjustable parameters are, e.g., attendance % in lectures and tutorials, how much time is spent on the exercises, how much help asked via Slack... • Maximum electric field enhancement by optimizing shape of metal nanoparticles. ### Bounds Maximize$f(x_0, x_1)$with$x_0\in\mathbb{R}$and$2\leq x_1 \leq 4$.$\Leftrightarrow$maximize$g(y_0, y_1) = f[y_0, 3 + \tanh(y_1)]$with$y_0,y_1\in\mathbb{R}$. Already implemented in most optimization packages. ### Equality constraints Maximize$f(x_0, x_1)$with$x_0, x_1 \in \mathbb{R}$subject to$x_0+x_1 = 4$.$\Leftrightarrow$maximize$g(x_0, x_1, \lambda) = f(x_0, x_1) - \lambda (x_0 + x_1 - 4)$(method of Lagrange multipliers) ### Inequality constraints Maximize$f(x_0, x_1)$with$x_0, x_1 \in \mathbb{R}$subject to$x_0^2-x_1 < 4$.$\Leftrightarrow$maximize$f(x_0, x_1, s)$subject to$x_0^2-x_1-4 = -s^2$, where$s\in\mathbb{R}$.$\Leftrightarrow$maximize$g(x_0, x_1, s, \lambda) = f(x_0, x_1) - \lambda (x_0^2-x_1-4+s^2)$(method of Slack variables) ## Gradient free optimization = numerical optimization of a function$f(x_0, x_1, \dots, x_{N-1})$without access to its gradient. Multiple available algorithms, e.g., ### Example: Coordinate descent A gradient-free optimization algorithm. Let's optimize$f(x, y) = 5x^2-6xy+5y^2$starting from$(x_0,y_0)=(-0.5,-1.0)$. 1. Find a coordinate$x_1$in x-direction for which$f(x_1, y_0) > f(x_0, y_0)$. 2. Find a coordinate$y_1$for which$f(x_1, y_1)>f(x_1, y_0)$. 3. Goto step 1 until converged. ## Gradient descent Idea is simply to use the information of$\nabla f(\mathbf{\overline x})$to maximize/minimize the function. 1. Pick initial condition$\mathbf{\overline x}$2. Compute$\nabla f(\mathbf{\overline x})$3. Set$x \leftarrow x - \gamma \nabla f(\mathbf{\overline x})$for some clever choise of$\gamma$and GOTO step 2. Continue until convergence. ## Optimization tools in Python ### scipy.optimize  import scipy.optimize as sciopt def fun(x): return x[0]**2+2*x[1]**2 res = sciopt.minimize(fun, x0=[0.1,0.2], method='Nelder-Mead') print( res.x ) print( res.fun )  ### nlopt  import nlopt def fun(x, grad): grad[0] = 2*x[0] grad[1] = 4*x[1] return x[0]**2+2*x[1]**2 opt = nlopt.opt(nlopt.LD_MMA, 2) opt.set_min_objective(fun) opt.set_xtol_abs(1e-4) xopt = opt.optimize([0.1, 0.2]) fopt = opt.last_optimum_value()  ### tensorflow  import tensorflow as tf x0 = np.array([0.1, 0.2]) x = tf.Variable(x0, dtype=tf.float64) fun = x[0]**2 + 2*x[1]**2 opt = tf.train.GradientDescentOptimizer(learning_rate=0.01) minimizer = opt.minimize(fun) sess = tf.Session() sess.run( tf.global_variables_initializer() ) sess.run(minimizer) xopt = sess.run(x) while np.linalg.norm(xopt - x0) > 1e-4: x0 = xopt sess.run(minimizer) xopt = sess.run( x ) print(xopt)  ## Curve fitting = minimizing error between data$(\mathbf{\overline x}, \mathbf{\overline y})$and parametrized function$f[\mathbf{\overline c}](\mathbf{\overline x})$. • Typical error functions: • Least squares:$\mathrm{err} = \sum\limits_i \Vert \mathbf{\overline y}_i - f(\mathbf{\overline x}_i) \Vert^2$• Absolute error:$\mathrm{err} = \sum\limits_i \Vert \mathbf{\overline y}_i - f(\mathbf{\overline x}_i) \Vert$• Relative error:$\mathrm{err} = \sum\limits_i \frac{\Vert \mathbf{\overline y}_i - f(\mathbf{\overline x}_i) \Vert}{\Vert \mathbf{\overline y}_i \Vert}\$
Division-by-zero warning!

### scipy.optimize.curve_fit


from scipy.optimize import curve_fit

xdata, ydata = generate_data()

def fun(x, a, b):
return a * np.sin(b*x)

popt, _ = curve_fit(fun, xdata, ydata)


### tensorflow


import tensorflow as tf

xdata_tmp, ydata_tmp = generate_data()

xdata = tf.Variable(xdata_tmp, dtype=tf.float64, trainable=False)
ydata = tf.Variable(xdata_tmp, dtype=tf.float64, trainable=False)

a = tf.Variable(0.2, dtype=tf.float64)
b = tf.Variable(0.1, dtype=tf.float64)

fitted_function = a * tf.sin( xdata * b )

error =  tf.reduce_sum( tf.abs( fitted_function - ydata )**2 )

opt = tf.train.RMSPropOptimizer(learning_rate = 0.01)
minimizer = opt.minimize(error)

sess = tf.Session()
sess.run( tf.global_variables_initializer() )

for i in range(100):
sess.run( minimizer )

aopt, bopt = sess.run([a,b])


## Feedforward neural network

= a many-parameter nonlinear function which is fitted to data

• Data is split to training and validation sets.
• Network is fitted only to training data
• Validation data is used to check whether the network generalizes beyond the training data

Check blackboard for more details.