FYS-4096 Computational Physics

Introduction to the Art of Numerical Experiments

By Janne Solanpää

Table of contents

Week 1: Basics of scientific computing
Course contents
Grading
Linux 101
Python refresher
Crash course on Git
Automating tasks
Numerical experiments workflow
Software design in Computational Physics
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
Word of warning about chaos
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 Google for help?
How to ask for help and report a bug?
Benchmarking your code
Profiling your code
Code optimization in HPC
Week 14: Optimization
Optimization
Gradient free optimization
Gradient based optimization
Curve fitting
Feedforward neural networks

Course contents

Core content (grade 1) Complementary knowledge (grade 3) Specialist knowledge (grade 5)
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
Core content (grade 1) Complementary knowledge (grade 3) Specialist knowledge (grade 5)
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
Core content (grade 1) Complementary knowledge (grade 3) Specialist knowledge (grade 5)
Correlations and spectral analysis Signal decomposition Behind the noise
Machine learning Neural networks Pitfalls and catastrophic failures

Communication channels

Course books

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?
YES! You will explain your solutions to your peers.

Project work

How many?
2
When are topics released?
February 19th
April 16th
Deadlines?
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 %

Grades

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

What is Computational Physics?

Computational Physics is...

Computational Physics is...

Computational Physics is...

Experimental physics ... on a computer

Linux 101

Why Linux?

Installation instructions for this course (1/2)

  1. Go to OnTheHub and login with your TUT credentials.
  2. Select tab VMware
  3. Select product corresponding to your operating system:
    Windows ↔ VMware Workstation 14
    macOS ↔ VMware Fusion 10
    Linux ↔ VMware Workstation 14
  4. Add to cart, order (free), download, and install

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.

Desktop environment

Shell 101

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

Live demo

Python refresher

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
  • Place a README file at root of your git repo
  • 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

Designing software for Scientific Computing

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
  • Document your code
  • Pick a coding style and stick to it
  • Do not comment/uncomment sections of the code to change the simulation
  • Automatic testing

Indentation with tabs vs spaces

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, ...)

Non-encapsulated design example

Encapsulated design example

Note: Bad design for other reasons (vectorization).

Python and Scientific Computing

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

Matplotlib

  • 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 calculus

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., \begin{equation} f'(x) = \frac{f(x+h) - f(x-h)}{2h}+\mathcal{O}(h^2) \end{equation}

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).

Python, FFT

Discrete Fourier transformations part 2

To obtain the original samples, we just need to do an inverse discrete fourier transformation: \begin{equation} 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] \end{equation}

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: \begin{equation} f(x) = \sum\limits_{k=-\infty}^\infty F_k \exp\left( -\frac{2\pi \mathrm{i}}{L} k x\right) \end{equation}

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}

More after linear algebra

Numerical Integrals

Trapezoidal rule

\begin{equation} \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] \end{equation}

Course book

Ch 5.9 Trapezoid rule

Simpson's rule

\begin{equation} \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} \end{equation}

Course book

Ch 5.10 Simpson's rule

Improper integrals

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

Change of variables, e.g., $x \to \arctan(x)$ yields: \begin{equation} \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 \end{equation} 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
    • If not, VEGAS bad

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

Course book ch. 7.2: Bisection method

Newton's method

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

Course book ch 7.3: Newton's method

Visualization

Matplotlib

  • 2D and simple 3D plots

Mayavi

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!

Example

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: \begin{equation} \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] \end{equation}

Numerical linear algebra with Python

Vectors

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


# 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

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


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

\begin{equation} \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] \end{equation} \begin{equation} \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] \end{equation}


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 ]

Element-wise product (Hadamart product)

\begin{equation} \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] \end{equation}


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

\begin{equation} \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] \end{equation}

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

\begin{equation} \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] \end{equation} \begin{equation} \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] \end{equation}


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

\begin{equation} \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] \end{equation} \begin{equation} \Vert \mathrm{\overline v} \Vert = \sqrt{ \sum\limits_{i=0}^2 \vert v_i \vert^2 } = \sqrt{1+9+25+1} = \sqrt{36} = 6 \end{equation} \begin{equation} \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} \end{equation} \begin{equation} \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 \end{equation}


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

\begin{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], \end{equation} 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., \begin{equation} \mathbf{A}^{-1} \mathbf{\overline v}. \end{equation}

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

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$

Sources:
R. Mäkinen's Lecture notes
rosettacode.org
LAPACK documentation

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

Dense and sparse matrices

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 \begin{equation} \mathbf{\overline u}(t) = \color{LimeGreen}{\exp(t \mathbf A)}\mathbf{\overline u}_0 \end{equation}
  • 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 \begin{equation} \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} \end{equation}

  • 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 \begin{equation} \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 \end{equation}
  • Tempting to look for an approximate solution $\mathbf{\overline w}_\mathrm{approx}$ in a subspace \begin{equation} \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\} \end{equation} 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})$: \begin{equation} \mathbf{\overline v}_0,\mathbf{\overline v}_1,\ldots,\mathbf{\overline v}_{m-1} \end{equation} 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 \begin{equation} \mathbf V_m = \left[ \mathbf{\overline v}_0, \ldots, \mathbf{\overline v}_{m-1}\right] \end{equation}
  • We also get an upper Hessenberg matrix \begin{equation} \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] \end{equation}

What we get from $m$ Arnoldi iterations?

  • By inspecting the Arnoldi scheme, we get the relation \begin{equation} \mathbf H_m = \mathbf V_m^\mathrm{T}\mathbf A\mathbf V_m. \end{equation} $\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})$ \begin{equation} \mathbf A \mathbf V_m = \mathbf V_m \mathbf H_m + H_{m,m-1}\mathbf{\overline v}_{m}\otimes \mathbf{\overline e}_{m} \end{equation}

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

The standard approximation is \begin{equation} \exp(t \mathbf A)\mathbf{\overline v} \approx \beta \mathbf V_m \exp(t \mathbf H_m) \mathbf{\overline e}_0, \end{equation} 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: \begin{equation} \exp(t \mathbf A)\mathbf{\overline v} \approx \beta \mathbf V_{m+1} \exp(t \mathbf{ \overline{H}}_{m}) \mathbf{\overline e}_0, \end{equation} where \begin{equation} \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], \end{equation} 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 \begin{equation} \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, \end{equation} 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})$

Other Krylov subspace techniques include, e.g.,

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})$ \begin{equation} \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}) \end{equation}
  • 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 \begin{equation} \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) \end{equation}
    • Possible outcomes of the measurement can be continuous or discrete, but they are always given by the eigenvalues $\lambda$ of the operator in question: \begin{equation} \hat{O}\psi_\lambda(\mathbf{\overline r}) = \lambda \psi_\lambda(\mathbf{\overline r}) \end{equation}
    • If the system is in a quantum state described by the wavefunction $\psi(\mathbf{\overline r})$, the expectation value of the measurement is \begin{equation} \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}) \end{equation}
  • 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 \begin{equation} P(\lambda) = \vert \langle \psi_\lambda \vert \psi \rangle \vert^2 \end{equation}

More material on quantum mechanics

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

Time-evolution in QM

  • Time-evolution of a quantum system follows the time-dependent Schrödinger equation (TDSE): \begin{equation} \mathrm{i}\frac{\partial}{\partial t}\psi(\mathbf{\overline r}, t) = \hat{H}(t) \psi(\mathbf{\overline r}, t) \end{equation} with the initial condition \begin{equation} \psi(\mathbf{\overline r}, t=0) = \phi(\mathbf{\overline r}) \end{equation}
  • $\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})$ , \begin{equation} \hat{H}_0 \phi(\mathbf{\overline r}) = E_\mathrm{gs}\phi(\mathbf{\overline r}) \end{equation}
  • If $\hat{H}(t) \equiv \hat{H}_0$, i.e., it has no time-dependence, the time-evolved wavefunction is given by \begin{equation} \psi(\mathbf{\overline r}, t) = e^{-\mathrm{i} t \hat{H_0} } \psi(\mathbf{\overline r}, t=0) \end{equation}
  • If $\hat{H}(t)$ depends on time, the time-evolved wavefunction is given by \begin{equation} \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), \end{equation} 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 \begin{equation} W = \mathbf{\overline r}\cdot \mathbf{\overline E}(t), \end{equation} 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: \begin{equation} H(t) = -\frac{1}{2} \frac{\partial^2}{\partial x^2} -\frac{1}{\sqrt{x^2+1}} + x \epsilon(t) \end{equation}

Discretizing PDEs with finite differences

  • Discretize = continuous $\to$ discrete
  • Let's first discretize 1D TISE \begin{equation} \left[ -\frac{1}{2}\frac{\mathrm{d}^2}{\mathrm{d}x^2} + V(x) \right] \psi(x) = E \psi(x), \end{equation} where $x\in[a,b]\subset \mathbb R$. We often use Dirichlet boundary conditions \begin{equation} \psi(a) = 0 \text{ and } \psi(b) = 0. \end{equation}
  • 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 \begin{equation} [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\} \end{equation}
  • Denote the values of the wavefunction at these coordinate by \begin{equation} \{\psi_a, \psi_0, \psi_1, \ldots, \psi_k = \psi[x_a+(k+1)\Delta x],\ldots,\psi_{N-1},\psi_N=\psi_b\} \end{equation}

  • Approximate derivatives with finite differences: \begin{equation} \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} \end{equation}
  • The TISE must hold for all points in our coordinate grid, and thus \begin{equation} -\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 \end{equation}
  • Remember that $\psi_a = \psi_b = 0$, and we get \begin{equation} \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] \end{equation}

  • The previous can be written more clearly as \begin{equation} \mathbf H \overline \psi = E \overline \psi, \end{equation} 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 \begin{equation} \frac{\partial}{\partial t}\psi(x,t) = -\mathrm{i} \hat{H}(t) \psi(x, t) \end{equation} with the initial condition \begin{equation} \psi(x, t=0) = \phi(x) \end{equation}
  • 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 \begin{equation} \frac{\partial}{\partial t}\overline \psi(t) = - \mathrm{i} \mathbf{H}(t) \overline \psi(t) \end{equation} with the initial condition \begin{equation} \overline \psi(t=0) = \overline \phi \end{equation} 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 \begin{equation} \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), \end{equation} 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"): \begin{equation} \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. \end{equation}
    • Multiplying the wavefunction by a mask function between every time-step: \begin{equation} \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. \end{equation}

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: Server room

Cluster autopsy: Blade

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
    

Skipping password-prompt

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
-bash: python3: command not found

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 --cpus-per-task=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

Tips and tricks

  • Run simulations in $WRKDIR (optimized for simulations)
  • Modify .bashrc for nicer command line prompt
  • Automate multiple simulations with bash or python-scripting:
    #!/usr/bin/env python3
    
    import subprocess
    
    job_script_template="""#!/bin/bash
    #SBATCH -J {jobname}
    #SBATCH -e out/{jobname}
    #SBATCH -o out/{jobname}
    #SBATCH -t 0:10:00
    #SBATCH -p serial
    
    ml python-env/intelpython3.5
    
    srun myscript --option1 {option1} --savefile data_{option1}
    """
    
    option1_values = [1, 2, 3.5, 50, 55]
    
    for option1 in option1_values:
        args = {'option1' : option1}
        with open("job.sub", "w") as file:
            job_script = job_script_template.format(**args)
            file.write(job_script)
    
        subprocess.run(["sbatch", "job.sub"])
    
    

Partial differential equations, part 2

What is a PDE?

  • An equation which contains partial derivatives of an unknown function $u(q_0,q_1,\ldots,q_{N-1})$ that depends on multiple variables
E.g., \begin{equation} \frac{\partial u}{\partial x}+u(x,y) + \frac{\partial^2 u}{\partial x \partial y} = \cos(x)+\sin(y) \end{equation} for an unknown function $u(x,y)$.

Course book ch. 19

Types of PDEs

Boundary value problems

  • We have a PDE for an unknown function $u(\mathbf{\overline r})$ in some domain $\Omega\subset \mathbb R^n$
  • The value of the unknown function or its derivative is known at the boundary of the domain, $\partial \Omega$ = boundary condition
  • Example:
    A CPU (copper slab) sitting on top of a circuit board and an aluminum heat sink attached. We turn on the CPU which generates heat with approximately 45 W power. What's the temperature distribution of the system?
    • Governing equation: \[ - \nabla \cdot \left[ k(\mathbf{\overline r}) \nabla T(\mathbf {\overline r}) \right] = Q(\mathbf {\overline r}), \] where $k$ is the thermal conductivity, $T$ the temperature, and $Q$ the volumetric heat generation rate.
    • Boundary conditions: \[ -k \nabla T \cdot \mathbf{\overline n} = h (T_{\mathrm{ext}}-T), \] where $h$ is the heat transfer coefficient for the boundary in question, and $T_\mathrm{ext}$ the external air temperature (assumed constant, i.e., we have some serious air flow).

Computed solution for the heatsink

Possible boundary conditions

Let's denote our unknown function by $u$. The boundary conditions are given by the physics we model.

  • Dirichlet boundary conditions \[ u(\mathbf{\overline r}) = D_{\mathrm{known}}(\mathbf{\overline r})\,\forall \mathbf{\overline r} \in \partial \Omega \]
  • Neumann boundary conditions \[ \nabla u(\mathbf{\overline r}) \cdot n(\mathbf{\overline r}) = N_\mathrm{known}(\mathbf{\overline r}), \] where $n$ is the boundary normal.
  • Robin boundary conditions \[ \alpha u(\mathbf{\overline r}) + \beta \nabla u(\mathbf{\overline r}) \cdot n(\mathbf{\overline r}) = R_{\mathrm{known}}(\mathbf{\overline r}), \] where $\alpha$ and $\beta$ are some constants.
  • Cauchy boundary conditions
    $ u(\mathbf{\overline r}) = D_{\mathrm{known}}(\mathbf{\overline r}) $ and $ \nabla u(\mathbf{\overline r}) \cdot n(\mathbf{\overline r}) = N_\mathrm{known}(\mathbf{\overline r}) $

Combining boundary conditions

  • We can impose different boundary conditions on different parts of the boundary.
  • Depending on the PDE in question, some (combinations of) boundary conditions may render the PDE unsolvable.
  • The function/derivative value on the boundary (e.g., $D_\mathrm{known}$) is often called Dirichlet/Neumann/Robin/Cauchy data.

Eigenvalue problems

  • Already encountered these with time-independent Schrödinger equation
  • Problems such as \[ \left[\nabla^2 + \mathbf{\overline g}\cdot \nabla\right] u = \lambda u, \] where $u$ needs to satisfy some boundary conditions as in boundary value problems.
  • Often used for finding out resonance modes in systems (e.g., acoustics).

Example: Eigenmodes of sound pressure in the lecture hall

  • Given by the Helmholtz equation for the sound pressure $p$ inside the lecture hall (the domain $\Omega$), \[ \nabla^2 p = -k^2 p \] with Neumann boundary conditions on all surfaces $\nabla p \cdot \mathbf{\overline n} = 0$ (hard wall boundaries).

Initial value problems

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) \]

Weak formulation of PDEs

Weak formulation 101

  • PDEs are terribly difficult to solve
    $\rightarrow$ let's do linear algebra instead
  • Let's first work through the Poisson equation \[ \nabla^2 u(\mathbf{\overline r}) = \rho(\mathbf{\overline r})\,\, \forall \mathbf{\overline r}\in \Omega \] with boundary conditions $ u(\mathbf{\overline r}) = D(\mathbf{\overline r})\,\, \forall \mathbf{\overline r}\in \Gamma_D \subset \partial \Omega$ and $ \nabla u(\mathbf{\overline r}) \cdot \mathbf{\overline n}(\mathbf{\overline r}) = N(\mathbf{\overline r})\,\,\forall\mathbf{\overline r} \in \Gamma_N \subset \partial \Omega $

  1. $ \nabla^2 u(\mathbf{\overline r}) = \rho(\mathbf{\overline r}) $
  2. Multiply by a test function $v(\mathbf{\overline r})$ and integrate over the domain \[ \int\limits_\Omega v(\mathbf{\overline r}) \nabla^2 u(\mathbf{\overline r})\,\mathrm{d}V = \int\limits_\Omega v(\mathbf{\overline r}) \rho(\mathbf{\overline r})\,\mathrm{d}V \]
  3. Decrease the order of differentiation by using Green's first identity \[ -\int\limits_\Omega \nabla v \cdot \nabla u \,\mathrm{d}V + \int\limits_{\partial \Omega} v \nabla u \cdot \mathrm{d}\mathbf{\overline S} = \int\limits_\Omega v \rho\,\mathrm{d}V \]
  4. Impose Neumann boundary conditions \[ -\int\limits_\Omega \nabla v \cdot \nabla u \,\mathrm{d}V + \int\limits_{\Gamma_N} v N \mathrm{d}S + \int\limits_{\Gamma_D} v \nabla u \cdot \mathrm{d}\mathbf{\overline S} = \int\limits_\Omega v \rho\,\mathrm{d}V \]
  5. Select appropriate function spaces:
    Test function space: $v \in \hat{V}=\left\{v \in H^1(\Omega)\,:\,v\vert_{\Gamma_D}= 0\right\}$
    Trial function space: $u \in V = \left\{ u \in H^1(\Omega)\, \,:u\vert_{\Gamma_D} = D\right\}$
  6. Profit: \[ -\int\limits_\Omega \nabla v \cdot \nabla u \,\mathrm{d}V + \int\limits_{\Gamma_N} v N \mathrm{d}S = \int\limits_\Omega v \rho\,\mathrm{d}V \]

Example continues

Finally we have a weak form of the original Poisson equation:

Find $u \in V$ such that \[ -\int\limits_\Omega \nabla v \cdot \nabla u \,\mathrm{d}V + \int\limits_{\Gamma_N} v N \mathrm{d}S = \int\limits_\Omega v \rho\,\mathrm{d}V \] for all $v\in \hat{V}$.

Numerical implementation of this goes as follows:

  1. Select finite-dimensional subspaces $\hat{V_h}\subset \hat{V}$ and $V_h \subset V$.
  2. Express our approximate solution as \[ u_h = \sum\limits_{j=0}^{N-1}x_j \phi_j \] where $\left(\{ \phi_j \}_{j=0}^{N}\right)$ forms a basis for $V_h$.
  3. Multiply the weak form equation by all basis vectors of the test function space $\hat{V}_h$, $\psi_i$.
  4. Profit: \[ \mathbf A \mathbf{\overline x} = \mathbf{\overline b}, \] where $\mathbf A_{ij} = -\int\limits_\Omega \nabla \psi_i \cdot \nabla \phi_j \,\mathrm{d}V$ and $b_i = \int\limits_\Omega \psi_i \rho \,\mathrm{d}V - \int\limits_{\Gamma_N} \psi_i N\,\mathrm{d}S$

Rules of thumb for boundary conditions and function spaces

Dirichlet boundary condition $u(\mathbf{\overline r}) = D(\mathbf{\overline r})\,\forall \mathbf{\overline r}\in \partial\Omega_D$
Test functions $v\in \hat{V}=\left\{ v \in H^1(\Omega) : v\vert_{\partial \Omega_D}=0 \right\}$
Trial functions $u\in V=\left\{ v \in H^1(\Omega) : v\vert_{\partial \Omega_D} = D \right\}$
Note: This trial space is not a vector space if $D\neq 0$. A better formulation would be to use 'Dirichlet lift ansatz', i.e., $u = D_+ + u_+$ where $D_+ = D$ on the boundary and zero inside the domain, and $u_+ \in V_{+}=\left\{ v \in H^1(\Omega) \right\}$, which is a vector space.
Neumann boundary condition $\frac{\partial u(\mathbf{\overline r})}{\partial \eta} = N(\mathbf{\overline r}) \,\forall\mathbf{\overline r}\in \partial \Omega_N$
Test functions $v\in \hat{V}=\left\{ v \in H^1(\Omega) \right\}$
Trial functions $u\in V=\left\{ v \in H^1(\Omega) \right\}$
Neumann BCs can be imposed directly in the surface integral term(s): $\int\limits_{\partial \Omega_N} v \frac{\partial u}{\partial \eta}\;\mathrm{d}S = \int\limits_{\partial \Omega_N} v N\;\mathrm{d}S$
Note: For boundary value problems with only Neumann BCs and only derivatives of the function in the equation (i.e., not the function itself), we need to add extra restrictions to $V$ in order to have a unique solution. E.g.,$ V=\left\{ u \in H^1(\Omega) : \int\limits_{\Omega} u\,\mathrm{d}V =0 \right\}$
Robin boundary condition $\alpha u + \beta \frac{\partial u(\mathbf{\overline r})}{\partial\eta} = R(\mathbf{\overline r}) \,\forall\mathbf{\overline r}\in \partial \Omega_R$
Test functions $v\in\hat{V}=\left\{ v \in H^1(\Omega)\right\}$
Trial functions $u\in V=\left\{ v \in H^1(\Omega)\right\}$
Robin BCs can be imposed directly in the surface integral term(s): $\int\limits_{\partial \Omega_R} v \frac{\partial u}{\partial \eta}\;\mathrm{d}S = \frac{1}{\beta}\int\limits_{\partial \Omega_R} v \left[R-\alpha u \right]\,\mathrm{d}S$

Weak formulation, cont.

The aim in weak formulation is to transform the PDE into a more convenient functional equation \[ L(u; v) = b(v), \] where $L$ is hopefully:

  • Symmetric wrt. $u$ and $v$.
  • Bilinear
  • Or at lest semilinear (linear wrt. the test functions $v$)

and $b$ is linear in $v$.

For more advanced examples on obtaining the weak form, check

Finite Element Method (FEM)

Recap from last week

We have the weak form of our PDE, e.g., the Poisson equation:

Find $u \in V$ such that \[ -\int\limits_\Omega \nabla v \cdot \nabla u \,\mathrm{d}V + \int\limits_{\Gamma_N} v N \mathrm{d}S = \int\limits_\Omega v \rho\,\mathrm{d}V \] for all $v\in \hat{V}$.

Numerical implementation of this goes as follows:

  1. Select finite-dimensional subspaces $\hat{V_h}\subset \hat{V}$ and $V_h \subset V$.
  2. Express our approximate solution as \[ u_h = \sum\limits_{j=0}^{N-1}x_j \phi_j \] where $\left(\{ \phi_j \}_{j=0}^{N}\right)$ forms a basis for $V_h$.
  3. Multiply the weak form equation by all basis vectors of the test function space $\hat{V}_h$, $\psi_i$.
  4. Profit: \[ \mathbf A \mathbf{\overline x} = \mathbf{\overline b}, \] where $\mathbf A_{ij} = -\int\limits_\Omega \nabla \psi_i \cdot \nabla \phi_j \,\mathrm{d}V$ and $b_i = \int\limits_\Omega \psi_i \rho \,\mathrm{d}V - \int\limits_{\Gamma_N} \psi_i N\,\mathrm{d}S$

The power of FEM

We've expressed our PDE in the weak form \[ \mathbf A \mathbf{\overline x} = \mathbf{\overline b}, \] where $\mathbf A_{ij} = -\int\limits_\Omega \nabla \psi_i \cdot \nabla \phi_j \,\mathrm{d}V$ and $b_i = \int\limits_\Omega \psi_i \rho \,\mathrm{d}V - \int\limits_{\Gamma_N} \psi_i N\,\mathrm{d}S$
$\Rightarrow$ choose $\phi_j \in V_h$ and $\psi_i \in \hat{V}_h$ so that $\mathbf A_{ij}$ is super sparse and very easy to compute.

Meshing

  • Meshing = approximating the domain $\Omega$ as a union of finite number of elementary elements \[ \Omega \approx \Omega_H = \bigcup\limits_{i}\Omega_i \]
  • Above, the elementary element of the mesh is a triangle

Different types of elements

Triangular (2D)

Quadrilateral (2D)

Tetrahedral (3D)

Mesh generation

  • Typical workflow:
    1. Define domain geometry ($\Omega$)
    2. Mesh geometry with automatic tools ($\Omega_h$)
    3. Refine mesh where appropriate
    4. Run a FEM simulation
    5. Refine mesh where error seems highest, goto step 4.
  • Many FEM software suites include basic domain generation tools, if not, use any 3D modeling software, e.g., Autocad.
  • Most FEM suites include basic meshing tools. If need for more advanced tools, use, e.g., Gmsh.
  • Often one needs to convert between file formats for meshes

Function space approximation

  • We've now meshed the domain with a mesh $\Omega_h = \bigcup\limits_{i=0}^{N-1} \Omega_i$
  • Select our approximate function spaces $V_h, \hat{V}_h$ by choosing a set of overlapping basis functions $u_j:\;\mathrm{span}(\{u_j\}) = V_h$ and $v_i:\;\mathrm{span}(\{v_i\}) = \hat{V}_h$ compact support: \[ v_i(\mathbf{\overline r}) = 0\;\forall \mathbf{\overline r} \not\in \mathrm{supp} (v_i) \subset \Omega_h \] and \[ u_j(\mathbf{\overline r}) = 0 \,\forall \mathbf{\overline r} \not\in \mathrm{supp}( u_j ) \subset \Omega_h \]
  • With these kinds of basis functions we have, e.g., $-\int\limits_\Omega \nabla v_i \cdot \nabla u_j \,\mathrm{d}V = 0 \;$ if $\mathrm{supp}(v_i) \cap \mathrm{supp}(u_j) = \emptyset$

Example basis: Hat functions in 1D

Example basis: Hat functions in 2D

FEM with FEniCS

FEniCS is a popular open-source (LGPLv3) computing platform for solving partial differential equations (PDEs). FEniCS enables users to quickly translate scientific models into efficient finite element code. With the high-level Python and C++ interfaces to FEniCS, it is easy to get started, but FEniCS offers also powerful capabilities for more experienced programmers.

Program/workflow structure:
  1. Create your domain and mesh it (you can use some other tools for this)
  2. Import mesh to FEniCS
  3. Define your unconstrained function space (i.e., no constraints due to BCs, but can empose, e.g., $\nabla\cdot \mathbf A = 0$)
    • Often you want to use the same unconstrained function space for both test and trial functions
  4. Define Dirichlet boundary conditions if any
  5. Write down the weak form equation (including Neumann and/or Robin BCs as surface integrals)
  6. Ask FEniCS to solve the weak form equation
  7. Save your solution to file
  8. Visualize using, e.g., matplotlib or paraview

See also the FEniCS tutorials

FEniCS example: Poisson equation

\[ \nabla^2 u(x,y) = \rho(x,y)\,\forall (x,y) \in [-1,1]\times [-1,1] = \Omega, \] where $\rho(x,y)=\exp[-(x^2+y^2)]$ and BCs are
  • $u(x=-1,y) = 0$
  • $u(x=1,y) = 3$
  • $\partial_y u(x, y=-1) = 1$
  • $\partial_y u(x,y=1) = 3u(x, y=1)$
The weak form is
Find $u \in V=\left\{ u\in H^1(\Omega): u(x=-1,y) = 0 \land u(x=1,y) = 3\right\}$ (trial function space) such that \[ -\int\limits_{\Omega} \nabla v \cdot \nabla u \,\mathrm{d}A + 3\!\int\limits_{y=1}v u \,\mathrm{d}s -\!\!\!\! \int\limits_{y=-1} v \,\mathrm{d}s = \int\limits_{\Omega} v \rho \,\mathrm{d}A \] for all $v \in \hat{V} = \left\{v\in H^1(\Omega) : v(x=\pm 1, y)=0\right\}$ (test function space).

FENICS EXAMPLE: Poisson equation cont.

Find $u \in V=\left\{ u\in H^1(\Omega): u(x=-1,y) = 0 \land u(x=1,y) = 3\right\}$ (trial function space) such that \[ -\int\limits_{\Omega} \nabla v \cdot \nabla u \,\mathrm{d}A + 3\!\int\limits_{y=1}v u \,\mathrm{d}s -\!\!\!\! \int\limits_{y=-1} v \,\mathrm{d}s = \int\limits_{\Omega} v \rho \,\mathrm{d}A \] for all $v \in \hat{V} = \left\{v\in H^1(\Omega) : v(x=\pm 1, y)=0\right\}$ (test function space).

FEniCS example: Eigenmodes of acoustic waves in a hard cylinder with one soft boundary

\[ \nabla^2 u(\mathbf{\overline r}) = k^2 u(\mathbf{\overline r})\,\forall \mathbf{\overline r} \in \text{Cylinder of radius 10cm, height 30cm}, \] BCs are
  • $u=0$ for one end of the cylinder (soft boundary)
  • $\frac{\partial u}{\partial \eta} = 0$ on other boundaries (hard boundaries)
The weak form is
Find $u_\lambda \in V=\left\{ u\in H^1(\Omega): u(\partial \Omega_{\text{soft boundary}}) = 0 \right\}, \lambda \in \mathbb{R}$ such that \[ -\int\limits_{\Omega} \nabla v \cdot \nabla u \,\mathrm{d}A = \lambda \int\limits_{\Omega} v u \,\mathrm{d}A \] for all $v \in \hat{V} = \left\{v\in H^1(\Omega) : v(\partial \Omega_{\text{soft boundary}}) = 0\right\}$ (test function space).

FENICS Example: Eigenmodes

Find $u_\lambda \in V=\left\{ u\in H^1(\Omega): u(\partial \Omega_{\text{soft boundary}}) = 0 \right\}, \lambda \in \mathbb{R}$ such that \[ -\int\limits_{\Omega} \nabla v \cdot \nabla u \,\mathrm{d}A = \lambda \int\limits_{\Omega} v u \,\mathrm{d}A \] for all $v \in \hat{V} = \left\{v\in H^1(\Omega) : v(\partial \Omega_{\text{soft boundary}}) = 0\right\}$ (test function space).
We use a premeshed geometry here. It's generated with Gmsh and converted to FEniCS's fileformat with
$ dolfin-convert mesh.msh mesh.xml

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{equation} \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} \end{equation} 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{equation} \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} \end{equation}
  • 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

Example: Electrons confined by 2D harmonic potential
wiggling as a Wigner molecule

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{equation} \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} \end{equation} 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

Chua's strange attractor

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. \begin{equation} \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. \end{equation} Add them up!
  2. Choose your time-stepping method
    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{equation} \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} \end{equation}

Note: Implicit Euler becomes unstable very early!

  1. Choose your time-stepping method
    1. Implicit Euler
    2. Pretty much every other solver:
      \begin{equation} \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} \end{equation}
      \begin{equation} \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} \end{equation} 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

Stationary

Non-stationary

Trend

Seasonal / cyclic component

Irregular component

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
      • Quadratic
      • Cubic
      • Power law
      • Exponential
      • Piecewise linear

Curve fitting 101

Detrended using a power-law fit

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. This is your IMF
  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.

Original data vs imfs


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

Data detrended with EEMD

Comparison of power-law and EEMD detrending results

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

Verifying your code

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.

Verifying your code

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.

FF network with tensorflow

Next courses...