Overview of week 37

Overview of week 37

Iterative methods, Chapter 6

• Direct solvers such as Gauss elimination and LU decomposition discussed last week.
• Iterative solvers such as Basic iterative solvers, Jacobi, Gauss-Seidel, Successive over-relaxation. These methods are easy to parallelize, as we will se later. Much used in solutions of partial differential equations.
• Other iterative methods such as Krylov subspace methods with Generalized minimum residual (GMRES) and Conjugate gradient etc will not be discussed.

Iterative methods, Jacobi's method

 
$$ \hat{A}{\bf x}={\bf b},$$

 
where \( \hat{A} \) is a matrix and \( {\bf x} \) and \( {\bf b} \) are vectors. The vector \( {\bf x} \) is the unknown.

It is an iterative scheme where we start with a guess for the unknown, and after \( k+1 \) iterations we have

 
$$ {\bf x}^{(k+1)}= \hat{D}^{-1}({\bf b}-(\hat{L}+\hat{U}){\bf x}^{(k)}),$$

 
with \( \hat{A}=\hat{D}+\hat{U}+\hat{L} \) and \( \hat{D} \) being a diagonal matrix, \( \hat{U} \) an upper triangular matrix and \( \hat{L} \) a lower triangular matrix.

If the matrix \( \hat{A} \) is positive definite or diagonally dominant, one can show that this method will always converge to the exact solution.

Iterative methods, Jacobi's method

 
$$ \begin{align} x_1^{(1)} =&(b_1-a_{12}x_2^{(0)} -a_{13}x_3^{(0)} - a_{14}x_4^{(0)})/a_{11} \nonumber \\ x_2^{(1)} =&(b_2-a_{21}x_1^{(0)} - a_{23}x_3^{(0)} - a_{24}x_4^{(0)})/a_{22} \nonumber \\ x_3^{(1)} =&(b_3- a_{31}x_1^{(0)} -a_{32}x_2^{(0)} -a_{34}x_4^{(0)})/a_{33} \nonumber \\ x_4^{(1)}=&(b_4-a_{41}x_1^{(0)} -a_{42}x_2^{(0)} - a_{43}x_3^{(0)})/a_{44}, \nonumber \end{align} $$

 
which after \( k+1 \) iterations reads

 
$$ \begin{align} x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k)} -a_{32}x_2^{(k)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k)} -a_{42}x_2^{(k)} - a_{43}x_3^{(k)})/a_{44}, \nonumber \end{align} $$

Iterative methods, Jacobi's method

 
$$ x_i^{(k+1)}=(b_i-\sum_{j=1, j\ne i}^{n}a_{ij}x_j^{(k)})/a_{ii} $$

 
or in an even more compact form as

 
$$ {\bf x}^{(k+1)}= \hat{D}^{-1}({\bf b}-(\hat{L}+\hat{U}){\bf x}^{(k)}),$$

 
with \( \hat{A}=\hat{D}+\hat{U}+\hat{L} \) and \( \hat{D} \) being a diagonal matrix, \( \hat{U} \) an upper triangular matrix and \( \hat{L} \) a lower triangular matrix.

Iterative methods, Gauss-Seidel's method

 
$$ \begin{align} x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k)} -a_{32}x_2^{(k)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k)} -a_{42}x_2^{(k)} - a_{43}x_3^{(k)})/a_{44}, \nonumber \end{align} $$

 
can be rewritten as

 
$$ \begin{align} x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k+1)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k+1)} -a_{32}x_2^{(k+1)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k+1)} -a_{42}x_2^{(k+1)} - a_{43}x_3^{(k+1)})/a_{44}, \nonumber \end{align} $$

 
which allows us to utilize the preceding solution (forward substitution). This improves normally the convergence behavior and leads to the Gauss-Seidel method!

Iterative methods, Gauss-Seidel's method

 
$$ \begin{align} x_1^{(k+1)} =&(b_1-a_{12}x_2^{(k)} -a_{13}x_3^{(k)} - a_{14}x_4^{(k)})/a_{11} \nonumber \\ x_2^{(k+1)} =&(b_2-a_{21}x_1^{(k+1)} - a_{23}x_3^{(k)} - a_{24}x_4^{(k)})/a_{22} \nonumber \\ x_3^{(k+1)} =&(b_3- a_{31}x_1^{(k+1)} -a_{32}x_2^{(k+1)} -a_{34}x_4^{(k)})/a_{33} \nonumber \\ x_4^{(k+1)}=&(b_4-a_{41}x_1^{(k+1)} -a_{42}x_2^{(k+1)} - a_{43}x_3^{(k+1)})/a_{44}, \nonumber \end{align} $$

 
to the following form

 
$$ x^{(k+1)}_i = \frac{1}{a_{ii}} \left( b_i - \sum_{j > i} a_{ij}x^{(k)}_j - \sum_{j < i}a_{ij}x^{(k+1)}_j \right), \quad i=1,2,\ldots,n. $$

The procedure is generally continued until the changes made by an iteration are below some tolerance.

The convergence properties of the Jacobi method and the Gauss-Seidel method are dependent on the matrix \( \hat{A} \). These methods converge when the matrix is symmetric positive-definite, or is strictly or irreducibly diagonally dominant. Both methods sometimes converge even if these conditions are not satisfied.

Iterative methods, Successive over-relaxation

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

 
where:

 
$$ \hat{A}=\begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}, \qquad \mathbf{x} = \begin{bmatrix} x_{1} \\ x_2 \\ \vdots \\ x_n \end{bmatrix} , \qquad \mathbf{b} = \begin{bmatrix} b_{1} \\ b_2 \\ \vdots \\ b_n \end{bmatrix}. $$

Iterative methods, Successive over-relaxation

 
$$ \hat{A} =\hat{D} + \hat{L} + \hat{U}, $$

 
where

 
$$ D = \begin{bmatrix} a_{11} & 0 & \cdots & 0 \\ 0 & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & a_{nn} \end{bmatrix}, \quad L = \begin{bmatrix} 0 & 0 & \cdots & 0 \\ a_{21} & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\a_{n1} & a_{n2} & \cdots & 0 \end{bmatrix}, \quad U = \begin{bmatrix} 0 & a_{12} & \cdots & a_{1n} \\ 0 & 0 & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \cdots & 0 \end{bmatrix}. $$

 
The system of linear equations may be rewritten as:

 
$$ (D+\omega L) \mathbf{x} = \omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x} $$

 
for a constant \( \omega > 1 \).

Iterative methods, Successive over-relaxation

 
$$ \mathbf{x}^{(k+1)} = (D+\omega L)^{-1} \big(\omega \mathbf{b} - [\omega U + (\omega-1) D ] \mathbf{x}^{(k)}\big). $$

 
However, by taking advantage of the triangular form of \( (D+\omega L) \), the elements of \( x^{(k+1)} \) can be computed sequentially using forward substitution:

 
$$ x^{(k+1)}_i = (1-\omega)x^{(k)}_i + \frac{\omega}{a_{ii}} \left(b_i - \sum_{j>i} a_{ij}x^{(k)}_j - \sum_{j < i} a_{ij}x^{(k+1)}_j \right),\quad i=1,2,\ldots,n. $$

 
The choice of relaxation factor is not necessarily easy, and depends upon the properties of the coefficient matrix. For symmetric, positive-definite matrices it can be proven that \( 0 < \omega < 2 \) will lead to convergence, but we are generally interested in faster convergence rather than just convergence.

Cubic Splines, Chapter 6

A spline function consists of polynomial pieces defined on subintervals. The different subintervals are connected via various continuity relations.

Assume we have at our disposal \( n+1 \) points \( x_0, x_1, \dots x_n \) arranged so that \( x_0 < x_1 < x_2 < \dots x_{n-1} < x_n \) (such points are called knots). A spline function \( s \) of degree \( k \) with \( n+1 \) knots is defined as follows

• On every subinterval \( [x_{i-1},x_i) \) $s$ is a polynomial of degree \( \le k \).
• \( s \) has \( k-1 \) continuous derivatives in the whole interval \( [x_0,x_n] \).

Splines

 
$$ s(x)=\left\{\begin{array}{cc} s_0(x)=a_0x+b_0 & x\in [x_0, x_1) \\ s_1(x)=a_1x+b_1 & x\in [x_1, x_2) \\ \dots & \dots \\ s_{n-1}(x)=a_{n-1}x+b_{n-1} & x\in [x_{n-1}, x_n] \end{array}\right. $$

In this case the polynomial consists of series of straight lines connected to each other at every endpoint. The number of continuous derivatives is then \( k-1=0 \), as expected when we deal with straight lines. Such a polynomial is quite easy to construct given \( n+1 \) points \( x_0, x_1, \dots x_n \) and their corresponding function values.

Splines

 
$$ s_{i-1}(x_i)= y_i = s_i(x_i), $$

 
with \( 1 \le i \le n-1 \). In total we have \( n \) polynomials of the type

 
$$ s_i(x)=a_{i0}+a_{i1}x+a_{i2}x^2+a_{i2}x^3, $$

 
yielding \( 4n \) coefficients to determine.

Splines

 
$$ y_i = s(x_i), $$

 
and

 
$$ s(x_{i+1})= y_{i+1}, $$

 
to be fulfilled. If we also assume that \( s' \) and \( s'' \) are continuous, then

 
$$ s'_{i-1}(x_i)= s'_i(x_i), $$

 
yields \( n-1 \) conditions. Similarly,

 
$$ s''_{i-1}(x_i)= s''_i(x_i), $$

 
results in additional \( n-1 \) conditions. In total we have \( 4n \) coefficients and \( 4n-2 \) equations to determine them, leaving us with \( 2 \) degrees of freedom to be determined.

Splines

 
$$ s''_{i}(x_i)= f_i, $$

 
and

 
$$ s''_{i}(x_{i+1})= f_{i+1}, $$

 
and setting up a straight line between \( f_i \) and \( f_{i+1} \) we have

 
$$ s_i''(x) = \frac{f_i}{x_{i+1}-x_i}(x_{i+1}-x)+ \frac{f_{i+1}}{x_{i+1}-x_i}(x-x_i), $$

 
and integrating twice one obtains

 
$$ s_i(x) = \frac{f_i}{6(x_{i+1}-x_i)}(x_{i+1}-x)^3+ \frac{f_{i+1}}{6(x_{i+1}-x_i)}(x-x_i)^3 +c(x-x_i)+d(x_{i+1}-x). $$

Splines

 
$$ \begin{align} s_i(x) =&\frac{f_i}{6(x_{i+1}-x_i)}(x_{i+1}-x)^3+ \frac{f_{i+1}}{6(x_{i+1}-x_i)}(x-x_i)^3 \nonumber \\ +&(\frac{y_{i+1}}{x_{i+1}-x_i}-\frac{f_{i+1}(x_{i+1}-x_i)}{6}) (x-x_i)+ (\frac{y_{i}}{x_{i+1}-x_i}-\frac{f_{i}(x_{i+1}-x_i)}{6}) (x_{i+1}-x). \end{align} $$

Splines

 
$$ s'_{i-1}(x_i)= s'_i(x_i), $$

 
and set \( x=x_i \). Defining \( h_i=x_{i+1}-x_i \) we obtain finally the following expression

 
$$ h_{i-1}f_{i-1}+2(h_{i}+h_{i-1})f_i+h_if_{i+1}= \frac{6}{h_i}(y_{i+1}-y_i)-\frac{6}{h_{i-1}}(y_{i}-y_{i-1}), $$

 
and introducing the shorthands \( u_i=2(h_{i}+h_{i-1}) \), \( v_i=\frac{6}{h_i}(y_{i+1}-y_i)-\frac{6}{h_{i-1}}(y_{i}-y_{i-1}) \), we can reformulate the problem as a set of linear equations to be solved through e.g., Gaussian elemination

Splines

 
$$ \left[\begin{array}{cccccccc} u_1 & h_1 &0 &\dots & & & & \\ h_1 & u_2 & h_2 &0 &\dots & & & \\ 0 & h_2 & u_3 & h_3 &0 &\dots & & \\ \dots& & \dots &\dots &\dots &\dots &\dots & \\ &\dots & & &0 &h_{n-3} &u_{n-2} &h_{n-2} \\ & && & &0 &h_{n-2} &u_{n-1} \end{array}\right] \left[\begin{array}{c} f_1 \\ f_2 \\ f_3\\ \dots \\ f_{n-2} \\ f_{n-1} \end{array} \right] = \left[\begin{array}{c} v_1 \\ v_2 \\ v_3\\ \dots \\ v_{n-2}\\ v_{n-1} \end{array} \right]. $$

 
Note that this is a set of tridiagonal equations and can be solved through only \( O(n) \) operations.

Splines

spline(double x[], double y[], int n, double yp1,
         double yp2, double y2[])

This function takes as input \( x[0,..,n - 1] \) and \( y[0,..,n - 1] \) containing a tabulation \( y_i = f(x_i) \) with \( x_0 < x_1 < .. < x_{n - 1} \) together with the first derivatives of \( f(x) \) at \( x_0 \) and \( x_{n-1} \), respectively. Then the function returns \( y2[0,..,n-1] \) which contains the second derivatives of \( f(x_i) \) at each point \( x_i \). \( n \) is the number of points. This function provides the cubic spline interpolation for all subintervals and is called only once.

Splines

splint(double x[], double y[], double y2a[], int n,
                    double x, double *y)

which takes as input the tabulated values \( x[0,..,n - 1] \) and \( y[0,..,n - 1] \) and the output y2a[0,..,n - 1] from \( spline \). It returns the value \( y \) corresponding to the point \( x \).

Overview of week 38

Overview of week 38

Eigenvalue Solvers

 
$$ {\bf A}{\bf x^{(\nu)}} = \lambda^{(\nu)}{\bf x^{(\nu)}}, $$

where \( \lambda^{(\nu)} \) are the eigenvalues and \( {\bf x^{(\nu)}} \) the corresponding eigenvectors. Unless otherwise stated, when we use the wording eigenvector we mean the right eigenvector. The left eigenvector is defined as

 
$$ {\bf x^{(\nu)}}_L{\bf A} = \lambda^{(\nu)}{\bf x^{(\nu)}}_L $$

 
The above right eigenvector problem is equivalent to a set of \( n \) equations with \( n \) unknowns \( x_i \).

Eigenvalue Solvers

 
$$ \left( {\bf A}-\lambda^{(\nu)} I \right) {\bf x^{(\nu)}} = 0, $$

with \( I \) being the unity matrix. This equation provides a solution to the problem if and only if the determinant is zero, namely

 
$$ \left| {\bf A}-\lambda^{(\nu)}{\bf I}\right| = 0, $$

which in turn means that the determinant is a polynomial of degree \( n \) in \( \lambda \) and in general we will have \( n \) distinct zeros.

Eigenvalue Solvers

 
$$ P(\lambda) = det(\lambda{\bf I}-{\bf A}), $$

 
or

 
$$ P(\lambda)= \prod_{i=1}^{n}\left(\lambda_i-\lambda\right). $$

 
The set of these roots is called the spectrum and is denoted as \( \lambda({\bf A}) \). If \( \lambda({\bf A})=\left\{\lambda_1,\lambda_2,\dots ,\lambda_n\right\} \) then we have

 
$$ det({\bf A})= \lambda_1\lambda_2\dots\lambda_n, $$

 
and if we define the trace of \( {\bf A} \) as

 
$$ Tr({\bf A})=\sum_{i=1}^n a_{ii}$$

 
then \( Tr({\bf A})=\lambda_1+\lambda_2+\dots+\lambda_n \).

Abel-Ruffini Impossibility Theorem

The content of this theorem is frequently misunderstood. It does not assert that higher-degree polynomial equations are unsolvable. In fact, if the polynomial has real or complex coefficients, and we allow complex solutions, then every polynomial equation has solutions; this is the fundamental theorem of algebra. Although these solutions cannot always be computed exactly with radicals, they can be computed to any desired degree of accuracy using numerical methods such as the Newton-Raphson method or Laguerre method, and in this way they are no different from solutions to polynomial equations of the second, third, or fourth degrees.

The theorem only concerns the form that such a solution must take. The content of the theorem is that the solution of a higher-degree equation cannot in all cases be expressed in terms of the polynomial coefficients with a finite number of operations of addition, subtraction, multiplication, division and root extraction. Some polynomials of arbitrary degree, of which the simplest nontrivial example is the monomial equation \( ax^n = b \), are always solvable with a radical.

Abel-Ruffini Impossibility Theorem

Today, in the modern algebraic context, we say that second, third and fourth degree polynomial equations can always be solved by radicals because the symmetric groups \( S_2, S_3 \) and \( S_4 \) are solvable groups, whereas \( S_n \) is not solvable for \( n \ge 5 \).

Eigenvalue Solvers

 
$$ {\bf D}= \left( \begin{array}{ccccccc} \lambda_1 & 0 & 0 & 0 & \dots &0 & 0 \\ 0 & \lambda_2 & 0 & 0 & \dots &0 &0 \\ 0 & 0 & \lambda_3 & 0 &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &\lambda_{n-1} & \\ 0 & \dots & \dots & \dots &\dots &0 & \lambda_n \end{array} \right). $$

If \( {\bf A} \) is real and symmetric then there exists a real orthogonal matrix \( {\bf S} \) such that

 
$$ {\bf S}^T {\bf A}{\bf S}= \mbox{diag}(\lambda_1,\lambda_2,\dots ,\lambda_n), $$

 
and for \( j=1:n \) we have \( {\bf A}{\bf S}(:,j) = \lambda_j {\bf S}(:,j) \).

Eigenvalue Solvers

We say that a matrix \( {\bf B} \) is a similarity transform of \( {\bf A} \) if

 
$$ {\bf B}= {\bf S}^T {\bf A}{\bf S}, \quad \mbox{where}\quad {\bf S}^T{\bf S}={\bf S}^{-1}{\bf S} ={\bf I}. $$

The importance of a similarity transformation lies in the fact that the resulting matrix has the same eigenvalues, but the eigenvectors are in general different.

Eigenvalue Solvers

 
$$ {\bf Ax}=\lambda{\bf x} \mbox{ and } {\bf B}= {\bf S}^T {\bf A}{\bf S}. $$

We multiply the first equation on the left by \( {\bf S}^T \) and insert \( {\bf S}^{T}{\bf S} ={\bf I} \) between \( {\bf A} \) and \( {\bf x} \). Then we get

 
$$ \begin{equation} ({\bf S^TAS})({\bf S^Tx})=\lambda{\bf S^Tx} , \end{equation} $$

which is the same as

 
$$ {\bf B} \left ( {\bf S^Tx} \right ) = \lambda \left ({\bf S^Tx} \right ). $$

The variable \( \lambda \) is an eigenvalue of \( {\bf B} \) as well, but with eigenvector \( {\bf S^Tx} \).

Eigenvalue Solvers

either apply subsequent similarity transformations (direct method) so that

 
$$ \begin{equation} {\bf S_N^T\dots S_1^TAS_1\dots S_N }={\bf D} , \end{equation} $$

or apply subsequent similarity transformations so that \( \bf A \) becomes tridiagonal (Householder) or upper/lower triangular (\( \bf QR \) method). Thereafter, techniques for obtaining eigenvalues from tridiagonal matrices can be used.

or use so-called power methods

or use iterative methods (Krylov, Lanczos, Arnoldi). These methods are popular for huge matrix problems.

Gaussian Elimination and Tridiagonal matrices, project 1

 
$$ -\frac{u_{i+1} -2u_i +u_{i-i}}{h^2}=f(x_i,u(x_i)), $$

 
with \( i=1,2,\dots, n \). We need to add to this system the two boundary conditions \( u(a) =u_0 \) and \( u(b) = u_{n+1} \). If we define a matrix

 
$$ {\bf A} = \frac{1}{h^2}\left(\begin{array}{cccccc} 2 & -1 & & & & \\ -1 & 2 & -1 & & & \\ & -1 & 2 & -1 & & \\ & \dots & \dots &\dots &\dots & \dots \\ & & &-1 &2& -1 \\ & & & &-1 & 2 \\ \end{array} \right) $$

 
and the corresponding vectors \( {\bf u} = (u_1, u_2, \dots,u_n)^T \) and \( {\bf f}({\bf u}) = f(x_1,x_2,\dots, x_n,u_1, u_2, \dots,u_n)^T \) we can rewrite the differential equation including the boundary conditions as a system of linear equations with a large number of unknowns

 
$$ {\bf A}{\bf u} = {\bf f}({\bf u}). $$

Student project 2, part a-b)

 
$$ -\frac{\hbar^2}{2 m} \left ( \frac{1}{r^2} \frac{d}{dr} r^2 \frac{d}{dr} - \frac{l (l + 1)}{r^2} \right )R(r) + V(r) R(r) = E R(r). $$

 
In our case \( V(r) \) is the harmonic oscillator potential \( (1/2)kr^2 \) with \( k=m\omega^2 \) and \( E \) is the energy of the harmonic oscillator in three dimensions. The oscillator frequency is \( \omega \) and the energies are

 
$$ E_{nl}= \hbar \omega \left(2n+l+\frac{3}{2}\right), $$

 
with \( n=0,1,2,\dots \) and \( l=0,1,2,\dots \).

Since we have made a transformation to spherical coordinates it means that \( r\in [0,\infty) \). The quantum number \( l \) is the orbital momentum of the electron.

Then we substitute \( R(r) = (1/r) u(r) \) and obtain

 
$$ -\frac{\hbar^2}{2 m} \frac{d^2}{dr^2} u(r) + \left ( V(r) + \frac{l (l + 1)}{r^2}\frac{\hbar^2}{2 m} \right ) u(r) = E u(r) . $$

The boundary conditions are \( u(0)=0 \) and \( u(\infty)=0 \).

Student project 2, part a-b)

 
$$ -\frac{\hbar^2}{2 m \alpha^2} \frac{d^2}{d\rho^2} u(\rho) + \left ( V(\rho) + \frac{l (l + 1)}{\rho^2} \frac{\hbar^2}{2 m\alpha^2} \right ) u(\rho) = E u(\rho) . $$

In project 2 we set \( l=0 \). Inserting \( V(\rho) = (1/2) k \alpha^2\rho^2 \) we end up with

 
$$ -\frac{\hbar^2}{2 m \alpha^2} \frac{d^2}{d\rho^2} u(\rho) + \frac{k}{2} \alpha^2\rho^2u(\rho) = E u(\rho) . $$

 
We multiply thereafter with \( 2m\alpha^2/\hbar^2 \) on both sides and obtain

 
$$ -\frac{d^2}{d\rho^2} u(\rho) + \frac{mk}{\hbar^2} \alpha^4\rho^2u(\rho) = \frac{2m\alpha^2}{\hbar^2}E u(\rho) . $$

Student project 2, part a-b)

 
$$ -\frac{d^2}{d\rho^2} u(\rho) + \frac{mk}{\hbar^2} \alpha^4\rho^2u(\rho) = \frac{2m\alpha^2}{\hbar^2}E u(\rho) . $$

 
The constant \( \alpha \) can now be fixed so that

 
$$ \frac{mk}{\hbar^2} \alpha^4 = 1, $$

 
or

 
$$ \alpha = \left(\frac{\hbar^2}{mk}\right)^{1/4}. $$

 
Defining

 
$$ \lambda = \frac{2m\alpha^2}{\hbar^2}E, $$

 
we can rewrite Schr\"odinger's equation as

 
$$ -\frac{d^2}{d\rho^2} u(\rho) + \rho^2u(\rho) = \lambda u(\rho) . $$

 
This is the first equation to solve numerically. In three dimensions the eigenvalues for \( l=0 \) are \( \lambda_0=3,\lambda_1=7,\lambda_2=11,\dots . \)

Student project 2, part a-b)

 
$$ \begin{equation} u''=\frac{u(\rho+h) -2u(\rho) +u(\rho-h)}{h^2} +O(h^2), \tag{1} \end{equation} $$

 
where \( h \) is our step. Next we define minimum and maximum values for the variable \( \rho \), \( \rho_{\mbox{min}}=0 \) and \( \rho_{\mbox{max}} \), respectively. You need to check your results for the energies against different values \( \rho_{\mbox{max}} \), since we cannot set \( \rho_{\mbox{max}}=\infty \).

Student project 2, part a-b)

 
$$ h=\frac{\rho_{\mbox{max}}-\rho_{\mbox{min}} }{n_{\mbox{step}}}. $$

 
Define an arbitrary value of \( \rho \) as

 
$$ \rho_i= \rho_{\mbox{min}} + ih\quad i=0,1,2,\dots , n_{\mbox{step}} $$

 
we can rewrite the Schr\"odinger equation for \( \rho_i \) as

 
$$ -\frac{u(\rho_i+h) -2u(\rho_i) +u(\rho_i-h)}{h^2}+\rho_i^2u(\rho_i) = \lambda u(\rho_i), $$

 
or in a more compact way

 
$$ -\frac{u_{i+1} -2u_i +u_{i-1}}{h^2}+\rho_i^2u_i=-\frac{u_{i+1} -2u_i +u_{i-1} }{h^2}+V_iu_i = \lambda u_i, $$

 
where \( V_i=\rho_i^2 \) is the harmonic oscillator potential.

Student project 2, part a-b)

 
$$ d_i=\frac{2}{h^2}+V_i, $$

 
and the non-diagonal matrix element

 
$$ e_i=-\frac{1}{h^2}. $$

 
In this case the non-diagonal matrix elements are given by a mere constant. All non-diagonal matrix elements are equal. With these definitions the Schr\"odinger equation takes the following form

 
$$ d_iu_i+e_{i-1}u_{i-1}+e_{i+1}u_{i+1} = \lambda u_i, $$

 
where \( u_i \) is unknown. We can write the latter equation as a matrix eigenvalue problem

 
$$ \begin{equation} \left( \begin{array}{ccccccc} d_1 & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & d_2 & e_2 & 0 & \dots &0 &0 \\ 0 & e_2 & d_3 & e_3 &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &d_{n_{\mbox{step}}-2} & e_{n_{\mbox{step}}-1}\\ 0 & \dots & \dots & \dots &\dots &e_{n_{\mbox{step}}-1} & d_{n_{\mbox{step}}-1} \end{array} \right) \left( \begin{array}{c} u_{1} \\ u_{2} \\ \dots\\ \dots\\ \dots\\ u_{n_{\mbox{step}}-1} \end{array} \right)=\lambda \left( \begin{array}{c} u_{1} \\ u_{2} \\ \dots\\ \dots\\ \dots\\ u_{n_{\mbox{step}}-1} \end{array} \right) \tag{2} \end{equation} $$

Student project 2, part a-b)

 
$$ \begin{equation} \left( \begin{array}{ccccccc} \frac{2}{h^2}+V_1 & -\frac{1}{h^2} & 0 & 0 & \dots &0 & 0 \\ -\frac{1}{h^2} & \frac{2}{h^2}+V_2 & -\frac{1}{h^2} & 0 & \dots &0 &0 \\ 0 & -\frac{1}{h^2} & \frac{2}{h^2}+V_3 & -\frac{1}{h^2} &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &\frac{2}{h^2}+V_{n_{\mbox{step}}-2} & -\frac{1}{h^2}\\ 0 & \dots & \dots & \dots &\dots &-\frac{1}{h^2} & \frac{2}{h^2}+V_{n_{\mbox{step}}-1} \end{array} \right) \tag{3} \end{equation} $$

Recall that the solutions are known via the boundary conditions at \( i=n_{\mbox{step}} \) and at the other end point, that is for \( \rho_0 \). The solution is zero in both cases.

Student project 2, part c)

 
$$ -\frac{\hbar^2}{2 m} \frac{d^2}{dr^2} u(r) + \frac{1}{2}k r^2u(r) = E^{(1)} u(r), $$

 
where \( E^{(1)} \) stands for the energy with one electron only. For two electrons with no repulsive Coulomb interaction, we have the following Schr\"odinger equation

 
$$ \left( -\frac{\hbar^2}{2 m} \frac{d^2}{dr_1^2} -\frac{\hbar^2}{2 m} \frac{d^2}{dr_2^2}+ \frac{1}{2}k r_1^2+ \frac{1}{2}k r_2^2\right)u(r_1,r_2) = E^{(2)} u(r_1,r_2) . $$

Student project 2, part c)

With no interaction this can be written out as the product of two single-electron wave functions, that is we have a solution on closed form.

We introduce the relative coordinate \( {\bf r} = {\bf r}_1-{\bf r}_2 \) and the center-of-mass coordinate \( {\bf R} = 1/2({\bf r}_1+{\bf r}_2) \). With these new coordinates, the radial Schr\"odinger equation reads

 
$$ \left( -\frac{\hbar^2}{m} \frac{d^2}{dr^2} -\frac{\hbar^2}{4 m} \frac{d^2}{dR^2}+ \frac{1}{4} k r^2+ kR^2\right)u(r,R) = E^{(2)} u(r,R). $$

Student project 2, part c)

 
$$ E^{(2)}=E_r+E_R. $$

We add then the repulsive Coulomb interaction between two electrons, namely a term

 
$$ V(r_1,r_2) = \frac{\beta e^2}{|{\bf r}_1-{\bf r}_2|}=\frac{\beta e^2}{r}, $$

 
with \( \beta e^2=1.44 \) eVnm.

Student project 2, part c)

 
$$ \left( -\frac{\hbar^2}{m} \frac{d^2}{dr^2}+ \frac{1}{4}k r^2+\frac{\beta e^2}{r}\right)\psi(r) = E_r \psi(r). $$

 
This equation is similar to the one we had previously in parts (a) and (b) and we introduce again a dimensionless variable \( \rho = r/\alpha \). Repeating the same steps, we arrive at

 
$$ -\frac{d^2}{d\rho^2} \psi(\rho) + \frac{mk}{4\hbar^2} \alpha^4\rho^2\psi(\rho)+\frac{m\alpha \beta e^2}{\rho\hbar^2}\psi(\rho) = \frac{m\alpha^2}{\hbar^2}E_r \psi(\rho) . $$

Student project 2, part c)

 
$$ \omega_r^2=\frac{1}{4}\frac{mk}{\hbar^2} \alpha^4, $$

 
and fix the constant \( \alpha \) by requiring

 
$$ \frac{m\alpha \beta e^2}{\hbar^2}=1 $$

 
or

 
$$ \alpha = \frac{\hbar^2}{m\beta e^2}. $$

Student project 2, part c)

 
$$ \lambda = \frac{m\alpha^2}{\hbar^2}E, $$

 
we can rewrite Schr\"odinger's equation as

 
$$ -\frac{d^2}{d\rho^2} \psi(\rho) + \omega_r^2\rho^2\psi(\rho) +\frac{1}{\rho}\psi(\rho) = \lambda \psi(\rho). $$

Student project 2, part c)

Here we will study the cases \( \omega_r = 0.01 \), \( \omega_r = 0.5 \), \( \omega_r =1 \), and \( \omega_r = 5 \) for the ground state only, that is the lowest-lying state.

With no repulsive Coulomb interaction you should get a result which corresponds to the relative energy of a non-interacting system. Make sure your results are stable as functions of \( \rho_{\mbox{max}} \) and the number of steps.

We are only interested in the ground state with \( l=0 \). We omit the center-of-mass energy.

For specific oscillator frequencies, the above equation has closed answers, see the article by M. Taut, Phys. Rev. A 48, 3561 - 3566 (1993). The article can be retrieved from the following web address http://prola.aps.org/abstract/PRA/v48/i5/p3561_1.

Eigenvalue Solvers

The first is the formal method, involving determinants and the \textit{characteristic polynomial}. This proves how many eigenvalues there are, and is the way most of you learned about how to solve the eigenvalue problem, but for matrices of dimensions greater than 2 or 3, it is rather impractical.

The other general approach is to use similarity or unitary tranformations to reduce a matrix to diagonal form. Almost always this is done in two steps: first reduce to for example a \textit{tridiagonal} form, and then to diagonal form. The main algorithms we will discuss in detail, Jacobi's and Householder's (so-called direct method) and Lanczos algorithms (an iterative method), follow this methodology.

Diagonalization methods, direct methods

Diagonalization methods

Important Matrix and vector handling packages

• LINPACK: package for linear equations and least square problems.
• LAPACK:package for solving symmetric, unsymmetric and generalized eigenvalue problems. From LAPACK's website http://www.netlib.org it is possible to download for free all source codes from this library. Both C/C++ and Fortran versions are available.
• BLAS (I, II and III): (Basic Linear Algebra Subprograms) are routines that provide standard building blocks for performing basic vector and matrix operations. Blas I is vector operations, II vector-matrix operations and III matrix-matrix operations. Highly parallelized and efficient codes, all available for download from http://www.netlib.org.

Eigenvalue Solvers, Jacobi

 
$$ {\bf S}= \left( \begin{array}{cccccccc} 1 & 0 & \dots & 0 & 0 & \dots & 0 & 0 \\ 0 & 1 & \dots & 0 & 0 & \dots & 0 & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & 0 & \dots \\ 0 & 0 & \dots & \cos\theta & 0 & \dots & 0 & \sin\theta \\ 0 & 0 & \dots & 0 & 1 & \dots & 0 & 0 \\ \dots & \dots & \dots & \dots & \dots & \dots & 0 & \dots \\ 0 & 0 & \dots & -\sin\theta & 0 & \dots & 1 & \cos\theta \\ 0 & 0 & \dots & 0 & \dots & \dots & 0 & 1 \end{array} \right) $$

with property \( {\bf S^{T}} = {\bf S^{-1}} \). It performs a plane rotation around an angle \( \theta \) in the Euclidean $n-$dimensional space.

Eigenvalue Solvers, Jacobi

 
$$ s_{kk}= s_{ll}=\cos\theta, s_{kl}=-s_{lk}= -\sin\theta, s_{ii}=-s_{ii}=1\quad i\ne k \quad i \ne l, $$

A similarity transformation

 
$$ {\bf B}= {\bf S}^T {\bf A}{\bf S}, $$

results in

 
$$ \begin{align*} b_{ik} &= a_{ik}\cos\theta - a_{il}\sin\theta , i \ne k, i \ne l \\ b_{il} &= a_{il}\cos\theta + a_{ik}\sin\theta , i \ne k, i \ne l \nonumber\\ b_{kk} &= a_{kk}\cos^2\theta - 2a_{kl}\cos\theta \sin\theta +a_{ll}\sin^2\theta\nonumber\\ b_{ll} &= a_{ll}\cos^2\theta +2a_{kl}\cos\theta sin\theta +a_{kk}\sin^2\theta\nonumber\\ b_{kl} &= (a_{kk}-a_{ll})\cos\theta \sin\theta +a_{kl}(\cos^2\theta-\sin^2\theta)\nonumber \end{align*} $$

The angle \( \theta \) is arbitrary. The recipe is to choose \( \theta \) so that all non-diagonal matrix elements \( b_{kl} \) become zero.

Eigenvalue Solvers, Jacobi

 
$$ \mbox{off}({\bf A}) = \sqrt{\sum_{i=1}^n\sum_{j=1,j\ne i}^n a_{ij}^2}. $$

 
To demonstrate the algorithm, we consider the simple \( 2\times 2 \) similarity transformation of the full matrix. The matrix is symmetric, we single out $ 1\le k < l \le n$ and use the abbreviations \( c=\cos\theta \) and \( s=\sin\theta \) to obtain

 
$$ \left( \begin{array}{cc} b_{kk} & 0 \\ 0 & b_{ll} \\\end{array} \right) = \left( \begin{array}{cc} c & -s \\ s &c \\\end{array} \right) \left( \begin{array}{cc} a_{kk} & a_{kl} \\ a_{lk} &a_{ll} \\\end{array} \right) \left( \begin{array}{cc} c & s \\ -s & c \\\end{array} \right). $$

Eigenvalue Solvers, Jacobi

 
$$ a_{kl}(c^2-s^2)+(a_{kk}-a_{ll})cs = b_{kl} = 0. $$

 
If \( a_{kl}=0 \) one sees immediately that \( \cos\theta = 1 \) and \( \sin\theta=0 \).

The Frobenius norm of an orthogonal transformation is always preserved. The Frobenius norm is defined as

 
$$ ||{\bf A}||_F = \sqrt{\sum_{i=1}^n\sum_{j=1}^n |a_{ij}|^2}. $$

 
This means that for our \( 2\times 2 \) case we have

 
$$ 2a_{kl}^2+a_{kk}^2+a_{ll}^2 = b_{kk}^2+b_{ll}^2, $$

 
which leads to

 
$$ \mbox{off}({\bf B})^2 = ||{\bf B}||_F^2-\sum_{i=1}^nb_{ii}^2=\mbox{off}({\bf A})^2-2a_{kl}^2, $$

 
since

 
$$ ||{\bf B}||_F^2-\sum_{i=1}^nb_{ii}^2=||{\bf A}||_F^2-\sum_{i=1}^na_{ii}^2+(a_{kk}^2+a_{ll}^2 -b_{kk}^2-b_{ll}^2). $$

 
This results means that the matrix \( {\bf A} \) moves closer to diagonal form for each transformation.

Eigenvalue Solvers, Jacobi

 
$$\cot 2\theta=\tau = \frac{a_{ll}-a_{kk}}{2a_{kl}}, $$

 
we obtain the quadratic equation (using \( \cot 2\theta=1/2(\cot \theta-\tan\theta) \)

 
$$ t^2+2\tau t-1= 0, $$

 
resulting in

 
$$ t = -\tau \pm \sqrt{1+\tau^2}, $$

 
and \( c \) and \( s \) are easily obtained via

 
$$ c = \frac{1}{\sqrt{1+t^2}}, $$

 
and \( s=tc \). Choosing \( t \) to be the smaller of the roots ensures that \( |\theta| \le \pi/4 \) and has the effect of minimizing the difference between the matrices \( {\bf B} \) and \( {\bf A} \) since

 
$$ ||{\bf B}-{\bf A}||_F^2=4(1-c)\sum_{i=1,i\ne k,l}^n(a_{ik}^2+a_{il}^2) +\frac{2a_{kl}^2}{c^2}. $$

Eigenvalue Solvers, Jacobi algo

Step 2. Setup a while-test where one compares the norm of the newly computed off-diagonal matrix elements

 
$$ \mbox{off}({\bf A}) = \sqrt{\sum_{i=1}^n\sum_{j=1,j\ne i}^n a_{ij}^2} > \epsilon. $$

Step 3. Now choose the matrix elements \( a_{kl} \) so that we have those with largest value, that is \( |a_{kl}|=\mbox{max}_{i\ne j} |a_{ij}| \).

Step 4. Compute thereafter \( \tau = (a_{ll}-a_{kk})/2a_{kl} \), \( \tan\theta \), \( \cos\theta \) and \( \sin\theta \).

Step 5. Compute thereafter the similarity transformation for this set of values \( (k,l) \), obtaining the new matrix \( {\bf B}= {\bf S}(k,l,\theta)^T {\bf A}{\bf S}(k,l,\theta) \).

Step 6. Compute the new norm of the off-diagonal matrix elements and continue till you have satisfied \( \mbox{off}({\bf B}) \le \epsilon \)

The convergence rate of the Jacobi method is however poor, one needs typically \( 3n^2-5n^2 \) rotations and each rotation requires \( 4n \) operations, resulting in a total of \( 12n^3-20n^3 \) operations in order to zero out non-diagonal matrix elements.

Jacobi's method, an example to convice you about the algorithm

 
$$ {\bf B} = \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & c & -s \\ 0 & s & c \end{array} \right)\left( \begin{array}{ccc} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{array} \right) \left( \begin{array}{ccc} 1 & 0 & 0 \\ 0 & c & s \\ 0 & -s & c \end{array} \right). $$

 
We will choose the angle \( \theta \) in order to have \( a_{23}=a_{32}=0 \). We get (symmetric matrix)

 
$$ {\bf B} =\left( \begin{array}{ccc} a_{11} & a_{12}c -a_{13}s& a_{12}s+a_{13}c \\ a_{12}c -a_{13}s & a_{22}c^2+a_{33}s^2 -2a_{23}sc& (a_{22}-a_{33})sc +a_{23}(c^2-s^2) \\ a_{12}s+a_{13}c & (a_{22}-a_{33})sc +a_{23}(c^2-s^2) & a_{22}s^2+a_{33}c^2 +2a_{23}sc \end{array} \right). $$

 
Note that \( a_{11} \) is unchanged! As it should.

Jacobi's method, an example to convice you about the algorithm

 
$$ {\bf B} =\left( \begin{array}{ccc} a_{11} & a_{12}c -a_{13}s& a_{12}s+a_{13}c \\ a_{12}c -a_{13}s & a_{22}c^2+a_{33}s^2 -2a_{23}sc& (a_{22}-a_{33})sc +a_{23}(c^2-s^2) \\ a_{12}s+a_{13}c & (a_{22}-a_{33})sc +a_{23}(c^2-s^2) & a_{22}s^2+a_{33}c^2 +2a_{23}sc \end{array} \right). $$

 
or

 
$$ \begin{align*} b_{11} &= a_{11} \\ b_{12} &= a_{12}\cos\theta - a_{13}\sin\theta , 1 \ne 2, 1 \ne 3 \\ b_{13} &= a_{13}\cos\theta + a_{12}\sin\theta , 1 \ne 2, 1 \ne 3 \nonumber\\ b_{22} &= a_{22}\cos^2\theta - 2a_{23}\cos\theta \sin\theta +a_{33}\sin^2\theta\nonumber\\ b_{33} &= a_{33}\cos^2\theta +2a_{23}\cos\theta \sin\theta +a_{22}\sin^2\theta\nonumber\\ b_{23} &= (a_{22}-a_{33})\cos\theta \sin\theta +a_{23}(\cos^2\theta-\sin^2\theta)\nonumber \end{align*} $$

 
We will fix the angle \( \theta \) so that \( b_{23}=0 \).

Jacobi's method, an example to convice you about the algorithm

 
$$ {\bf B} =\left( \begin{array}{ccc} b_{11} & b_{12}& b_{13} \\ b_{12}& b_{22}& 0 \\ b_{13}& 0& a_{33} \end{array} \right). $$

 
We repeat then assuming that \( b_{12} \) is the largest non-diagonal matrix element and get a new matrix

 
$$ {\bf C} = \left( \begin{array}{ccc} c & -s & 0 \\ s & c & 0 \\ 0 & 0 & 1 \end{array} \right)\left( \begin{array}{ccc} b_{11} & b_{12} & b_{13} \\ b_{12} & b_{22} & 0 \\ b_{13} & 0 & b_{33} \end{array} \right) \left( \begin{array}{ccc} c & s & 0 \\ -s & c & 0 \\ 0 & 0 & 1 \end{array} \right). $$

 
We continue this process till all non-diagonal matrix elements are zero (ideally). You will notice that performing the above operations that the matrix element \( b_{23} \) which was previous zero becomes different from zero. This is one of the problems which slows down the jacobi procedure.

Jacobi's method, an example to convice you about the algorithm

 
$$ \begin{align*} b_{ii} &= a_{ii}, i \ne k, i \ne l \\ b_{ik} &= a_{ik}\cos\theta - a_{il}\sin\theta , i \ne k, i \ne l \\ b_{il} &= a_{il}\cos\theta + a_{ik}\sin\theta , i \ne k, i \ne l \nonumber\\ b_{kk} &= a_{kk}\cos^2\theta - 2a_{kl}\cos\theta \sin\theta +a_{ll}\sin^2\theta\nonumber\\ b_{ll} &= a_{ll}\cos^2\theta +2a_{kl}\cos\theta \sin\theta +a_{kk}\sin^2\theta\nonumber\\ b_{kl} &= (a_{kk}-a_{ll})\cos\theta \sin\theta +a_{kl}(\cos^2\theta-\sin^2\theta)\nonumber \end{align*} $$

 
This is what we will need to code.

Jacobi code example

//  we have defined a matrix A and a matrix R for the eigenvector, both of dim n x n
//  The final matrix R has the eigenvectors in its row elements, it is set to one
//  for the diagonal elements in the beginning, zero else.
....
double tolerance = 1.0E-10;
int iterations = 0;
while ( maxnondiag > tolerance && iterations <= maxiter)
{
   int p, q;
   maxnondiag  = offdiag(A, p, q, n);
   Jacobi_rotate(A, R, p, q, n);
   iterations++;

...

Jacobi code example

//  the offdiag function
double offdiag(double **A, int p, int q, int n);
{
   double max;
   for (int i = 0; i < n; ++i)
   {
       for ( int j = i+1; j < n; ++j)
       {
           double aij = fabs(A[i][j]);
           if ( aij > max)
           {
              max = aij;  p = i; q = j;

   return max;

...

Jacobi code example

void Jacobi_rotate ( double ** A, double ** R, int k, int l, int n )
{
  double s, c;
  if ( A[k][l] != 0.0 ) {
    double t, tau;
    tau = (A[l][l] - A[k][k])/(2*A[k][l]);

    if ( tau >= 0 ) {
      t = 1.0/(tau + sqrt(1.0 + tau*tau));
    } else {
      t = -1.0/(-tau +sqrt(1.0 + tau*tau));


    c = 1/sqrt(1+t*t);
    s = c*t;
  } else {
    c = 1.0;
    s = 0.0;

Jacobi code example

  double a_kk, a_ll, a_ik, a_il, r_ik, r_il;
  a_kk = A[k][k];
  a_ll = A[l][l];
  A[k][k] = c*c*a_kk - 2.0*c*s*A[k][l] + s*s*a_ll;
  A[l][l] = s*s*a_kk + 2.0*c*s*A[k][l] + c*c*a_ll;
  A[k][l] = 0.0;  // hard-coding non-diagonal elements by hand
  A[l][k] = 0.0;  // same here
  for ( int i = 0; i < n; i++ ) {
    if ( i != k && i != l ) {
      a_ik = A[i][k];
      a_il = A[i][l];
      A[i][k] = c*a_ik - s*a_il;
      A[k][i] = A[i][k];
      A[i][l] = c*a_il + s*a_ik;
      A[l][i] = A[i][l];

Jacobi code example

    r_ik = R[i][k];
    r_il = R[i][l];

    R[i][k] = c*r_ik - s*r_il;
    R[i][l] = c*r_il + s*r_ik;

  return;
} // end of function jacobi_rotate

Week 39

Overview of week 39

Eigenvalue Solvers, Householder

 
$$ {\bf S}={\bf S}_1{\bf S}_2\dots{\bf S}_{n-2}, $$

each of which successively transforms one row and one column of \( {\bf A} \) into the required tridiagonal form. Only \( n-2 \) transformations are required, since the last two elements are already in tridiagonal form. In order to determine each \( {\bf S_i} \) let us see what happens after the first multiplication, namely,

 
$$ {\bf S}_1^T{\bf A}{\bf S}_1= \left( \begin{array}{ccccccc} a_{11} & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & a'_{22} &a'_{23} & \dots & \dots &\dots &a'_{2n} \\ 0 & a'_{32} &a'_{33} & \dots & \dots &\dots &a'_{3n} \\ 0 & \dots &\dots & \dots & \dots &\dots & \\ 0 & a'_{n2} &a'_{n3} & \dots & \dots &\dots &a'_{nn} \\ \end{array} \right) $$

where the primed quantities represent a matrix \( {\bf A}' \) of dimension \( n-1 \) which will subsequentely be transformed by \( {\bf S_2} \).

Eigenvalue Solvers, Householder

 
$$ \left ({\bf S}_{1}{\bf S}_{2} \right )^{T} {\bf A}{\bf S}_{1} {\bf S}_{2} = \left( \begin{array}{ccccccc} a_{11} & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & a'_{22} &e_2 & 0 & \dots &\dots &0 \\ 0 & e_2 &a''_{33} & \dots & \dots &\dots &a''_{3n} \\ 0 & \dots &\dots & \dots & \dots &\dots & \\ 0 & 0 &a''_{n3} & \dots & \dots &\dots &a''_{nn} \\ \end{array} \right) $$

{\bf Note that the effective size of the matrix on which we apply the transformation reduces for every new step. In the previous Jacobi method each similarity transformation is in principle performed on the full size of the original matrix.}

Eigenvalue Solvers, Householder

 
$$ a_{11}, a'_{22}, a''_{33}\dots a^{n-1}_{nn}, $$

and off-diagonal matrix elements

 
$$ e_1, e_2,e_3, \dots, e_{n-1}. $$

The resulting matrix reads

 
$$ {\bf S}^{T} {\bf A} {\bf S} = \left( \begin{array}{ccccccc} a_{11} & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & a'_{22} & e_2 & 0 & \dots &0 &0 \\ 0 & e_2 & a''_{33} & e_3 &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &a^{(n-1)}_{n-2} & e_{n-1}\\ 0 & \dots & \dots & \dots &\dots &e_{n-1} & a^{(n-1)}_{n-1} \end{array} \right) . $$

Eigenvalue Solvers, Householder

 
$$ {\bf S_{1}} = \left( \begin{array}{cc} 1 & {\bf 0^{T}} \\ {\bf 0}& {\bf P} \end{array} \right), $$

with \( {\bf 0^{T}} \) being a zero row vector, \( {\bf 0^{T}} = \{0,0,\cdots\} \) of dimension \( (n-1) \). The matrix \( {\bf P} \) is symmetric with dimension (\( (n-1) \times (n-1) \)) satisfying \( {\bf P}^2={\bf I} \) and \( {\bf P}^T={\bf P} \). A possible choice which fullfils the latter two requirements is

 
$$ {\bf P}={\bf I}-2{\bf u}{\bf u}^T, $$

where \( {\bf I} \) is the \( (n-1) \) unity matrix and \( {\bf u} \) is an \( n-1 \) column vector with norm \( {\bf u}^T{\bf u} \) (inner product).

Eigenvalue Solvers, Householder

 
$$ P_{ij}=\delta_{ij}-2u_iu_j, $$

where \( i \) and \( j \) range from \( 1 \) to \( n-1 \). Applying the transformation \( {\bf S}_1 \) results in

 
$$ {\bf S}_1^T{\bf A}{\bf S}_1 = \left( \begin{array}{cc} a_{11} & ({\bf Pv})^T \\ {\bf Pv}& {\bf A}' \end{array} \right) , $$

where \( {\bf v^{T}} = \{a_{21}, a_{31},\cdots, a_{n1}\} \) and {\bf P} must satisfy (\( {\bf Pv})^{T} = \{k, 0, 0,\cdots \} \). Then

 
$$ \begin{equation} {\bf Pv} = {\bf v} -2{\bf u}( {\bf u}^T{\bf v})= k {\bf e}, \tag{4} \end{equation} $$

 
with \( {\bf e^{T}} = \{1,0,0,\dots 0\} \).

Eigenvalue Solvers, Householder

 
$$ ({\bf Pv})^T{\bf Pv} = k^{2} = {\bf v}^T{\bf v}= |v|^2 = \sum_{i=2}^{n}a_{i1}^2, $$

which determines the constant $ k = \pm v$.

Eigenvalue Solvers, Householder

 
$$ {\bf v} - k{\bf e} = 2{\bf u}( {\bf u}^T{\bf v}), $$

 
and taking the scalar product of this equation with itself and obtain

 
$$ \begin{equation} 2( {\bf u}^T{\bf v})^2=(v^2\pm a_{21}v), \tag{5} \end{equation} $$

 
which finally determines

 
$$ {\bf u}=\frac{{\bf v}-k{\bf e}}{2( {\bf u}^T{\bf v})}. $$

 
In solving Eq. (5) great care has to be exercised so as to choose those values which make the right-hand largest in order to avoid loss of numerical precision. The above steps are then repeated for every transformations till we have a tridiagonal matrix suitable for obtaining the eigenvalues.

Eigenvalue Solvers, Householder, brute force

 
$$ {\bf A} = \left( \begin{array}{cccc} d_{1} & e_{1} & 0 & 0 \\ e_{1} & d_{2} & e_{2} & 0 \\ 0 & e_{2} & d_{3} & e_{3} \\ 0 & 0 & e_{3} & d_{4} \end{array} \right). $$

As a first observation, if any of the elements \( e_{i} \) are zero the matrix can be separated into smaller pieces before diagonalization. Specifically, if \( e_{1} = 0 \) then \( d_{1} \) is an eigenvalue.

Eigenvalue Solvers, Householder

 
$$ {\bf S_{1}} = \left( \begin{array}{cccc} \cos \theta & 0 & 0 & \sin \theta\\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \cos \theta & 0 & 0 & \cos \theta \end{array} \right) $$

Eigenvalue Solvers, Householder

 
$$ {\bf S_{1}^{T} A S_{1}} = {\bf A'} = \left( \begin{array}{cccc} d'_{1} & e'_{1} & 0 & 0 \\ e'_{1} & d_{2} & e_{2} & 0 \\ 0 & e_{2} & d_{3} & e'{3} \\ 0 & 0 & e'_{3} & d'_{4} \end{array} \right) $$

produces a matrix where the primed elements in \( {\bf A'} \) have been changed by the transformation whereas the unprimed elements are unchanged. If we now choose \( \theta \) to give the element \( a_{21}^{'} = e^{'}= 0 \) then we have the first eigenvalue \( = a_{11}^{'} = d_{1}^{'} \). (This is actually what you are doing in project 2!!)

Eigenvalue Solvers, Householder's and Francis' algorithm

What we see here is just a special case of the more general procedure developed by Francis in two articles in 1961 and 1962. Using Jacobi's method is not very efficient ether.

The algorithm is based on the so-called {\bf QR} method (or just {\bf QR}-algorithm). It follows from a theorem by Schur which states that any square matrix can be written out in terms of an orthogonal matrix \( \hat{Q} \) and an upper triangular matrix \( \hat{U} \). Historically \( R \) was used instead of \( U \) since the wording right triangular matrix was first used.

Eigenvalue Solvers, Householder's and Francis' algorithm

Schur's theorem

 
$$ \hat{A} = \hat{Q}\hat{U}, $$

 
is used to rewrite any square matrix into a unitary matrix times an upper triangular matrix. We say that a square matrix is similar to a triangular matrix.

Householder's algorithm which we have derived is just a special case of the general Householder algorithm. For a symmetric square matrix we obtain a tridiagonal matrix.

There is a corollary to Schur's theorem which states that every Hermitian matrix is unitarily similar to a diagonal matrix.

Eigenvalue Solvers, Householder's and Francis' algorithm

 
$$ \hat{A}\hat{Q} = \hat{Q}\hat{U}\hat{Q}, $$

 
and multiply from the left with \( \hat{Q}^{-1} \) we get

 
$$ \hat{Q}^{-1}\hat{A}\hat{Q} = \hat{B}=\hat{U}\hat{Q}, $$

 
where the matrix \( \hat{B} \) is a similarity transformation of \( \hat{A} \) and has the same eigenvalues as \( \hat{B} \).

Eigenvalue Solvers, Householder's and Francis' algorithm

 
$$ \hat{A} = \hat{Q}\hat{U}, $$

 
and multiply from the left with \( \hat{Q}^{-1} \) resulting in

 
$$ \hat{Q}^{-1}\hat{A} = \hat{U}. $$

 
Suppose that \( \hat{Q} \) consists of a series of planar Jacobi like rotations acting on sub blocks of \( \hat{A} \) so that all elements below the diagonal are zeroed out

 
$$ \hat{Q}=\hat{R}_{12}\hat{R}_{23}\dots\hat{R}_{n-1,n}. $$

Eigenvalue Solvers, Householder's and Francis' algorithm

 
$$ \hat{R}_{12} = \left( \begin{array}{ccccccccc} c&s &0 &0 &0 & \dots &0 & 0 & 0\\ -s&c &0 &0 &0 & \dots &0 & 0 & 0\\ 0&0 &1 &0 &0 & \dots &0 & 0 & 0\\ \dots&\dots &\dots &\dots &\dots &\dots \\ 0&0 &0 & 0 & 0 & \dots &1 &0 &0 \\ 0&0 &0 & 0 & 0 & \dots &0 &1 &0 \\ 0&0 &0 & 0 & 0 & \dots &0 &0 & 1 \end{array} \right) $$

Eigenvalue Solvers, Householder's and Francis' algorithm

 
$$ \hat{U} = \left( \begin{array}{ccccccccc} x&x &x &0 &0 & \dots &0 & 0 & 0\\ 0&x &x &x &0 & \dots &0 & 0 & 0\\ 0&0 &x &x &x & \dots &0 & 0 & 0\\ \dots&\dots &\dots &\dots &\dots &\dots \\ 0&0 &0 & 0 & 0 & \dots &x &x &x \\ 0&0 &0 & 0 & 0 & \dots &0 &x &x \\ 0&0 &0 & 0 & 0 & \dots &0 &0 & x \end{array} \right) $$

 
which has a second superdiagonal.

Eigenvalue Solvers, Householder's and Francis' algorithm

 
$$ \hat{Q}^{-1}\hat{A}\hat{Q} = \hat{B}, $$

 
from Schur's theorem the matrix \( \hat{B} \) is triangular and the eigenvalues the same as those of \( \hat{A} \) and are given by the diagonal matrix elements of \( \hat{B} \). Why?

Our matrix \( \hat{B}=\hat{U}\hat{Q} \).

Another iterative procedure

The eigenvalues of a matrix can be obtained using the characteristic polynomial

 
$$ P(\lambda) = det(\lambda{\bf I}-{\bf A})= \prod_{i=1}^{n}\left(\lambda_i-\lambda\right), $$

 
which rewritten in matrix form reads

 
$$ P(\lambda)= \left( \begin{array}{ccccccc} d_1-\lambda & e_1 & 0 & 0 & \dots &0 & 0 \\ e_1 & d_2-\lambda & e_2 & 0 & \dots &0 &0 \\ 0 & e_2 & d_3-\lambda & e_3 &0 &\dots & 0\\ \dots & \dots & \dots & \dots &\dots &\dots & \dots\\ 0 & \dots & \dots & \dots &\dots &d_{N_{\mbox{step}}-2}-\lambda & e_{N_{\mbox{step}}-1}\\ 0 & \dots & \dots & \dots &\dots &e_{N_{\mbox{step}}-1} & d_{N_{\mbox{step}}-1}-\lambda \end{array} \right) $$

Eigenvalue Solvers, Householder

 
$$ P_k(\lambda)=(d_k-\lambda)P_{k-1}(\lambda)-e_{k-1}^2P_{k-2}(\lambda). $$

 
Together with the starting values \( P_1(\lambda) \) and \( P_2(\lambda) \) and good root searching methods we arrive at an efficient computational scheme for finding the roots of \( P_n(\lambda) \). However, for large matrices this algorithm is rather inefficient and time-consuming.

Eigenvalue Solvers, Householder functions

C++:

void trd2(double $**$a, int n, double d[], double e[])
void tqli(double d[], double[], int n, double $**$z)

Fortran:

CALL tred2(a, n, d, e)
CALL tqli(d, e, n, z)

Using Lapack to solve linear algebra and eigenvalue problems

Lapack example

#include <iostream>
#define MAX 10
using namespace std;
int main(){
   // Values needed for dgesv
   int n;
   int nrhs = 1;
   double a[MAX][MAX];
   double b[1][MAX];
   int lda = MAX;
   int ldb = MAX;
   int ipiv[MAX];
   int info;
   // Other values
   int i,j;

Lapack example

   // Read the values of the matrix
   cout << "Enter n \n";
   cin >> n;
   cout << "On each line type a row of the matrix A followed by one element of b:\n";
   for(i = 0; i < n; i++){
     cout << "row " << i << " ";
     for(j = 0; j < n; j++)std::cin >> a[j][i];
     cin >> b[0][i];

Lapack example

   // Solve the linear system
   dgesv(n, nrhs, &a[0][0], lda, ipiv, &b[0][0], ldb, &info);
   // Check for success
   if(info == 0)
   {
      // Write the answer
      cout << "The answer is\n";
      for(i = 0; i < n; i++)
        cout << "b[" << i << "]\t" << b[0][i] << "\n";

   else
   {
      // Write an error message
      cerr << "dgesv returned error " << info << "\n";

   return info;

Lanczos' iteration

Lanczos' algorithm generates a sequence of real tridiagonal matrices \( T_k \) of dimension \( k\times k \) with \( k\le n \), with the property that the extremal eigenvalues of \( T_k \) are progressively better estimates of \( \hat{A} \)' extremal eigenvalues.

The method converges to the extremal eigenvalues.

The similarity transformation is

 
$$ \hat{T}= \hat{Q}^{T}\hat{A}\hat{Q}, $$

 
with the first vector \( \hat{Q}\hat{e}_1=\hat{q}_1 \).

Lanczos' iteration

 
$$ \hat{T}= \hat{Q}^{T}\hat{A}\hat{Q}, $$

 
with the first vector \( \hat{Q}\hat{e}_1=\hat{q}_1 \). We can write out the matrix \( \hat{Q} \) in terms of its column vectors

 
$$ \hat{Q}=\left[\hat{q}_1\hat{q}_2\dots\hat{q}_n\right]. $$

Lanczos' iteration

 
$$ \hat{T}= \hat{Q}^{T}\hat{A}\hat{Q}, $$

 
can be written as

 
$$ \hat{T} = \left(\begin{array}{cccccc} \alpha_1& \beta_1 & 0 &\dots & \dots &0 \\ \beta_1 & \alpha_2 & \beta_2 &0 &\dots &0 \\ 0& \beta_2 & \alpha_3 & \beta_3 & \dots &0 \\ \dots& \dots & \dots &\dots &\dots & 0 \\ \dots& & &\beta_{n-2} &\alpha_{n-1}& \beta_{n-1} \\ 0& \dots &\dots &0 &\beta_{n-1} & \alpha_{n} \\ \end{array} \right) $$

Lanczos' iteration

 
$$ \hat{T}= \hat{Q}^{T}\hat{A}\hat{Q}, $$

 
as

 
$$ \hat{Q}\hat{T}= \hat{A}\hat{Q}, $$

 
and if we equate columns (recall from the previous slide)

 
$$ \hat{T} = \left(\begin{array}{cccccc} \alpha_1& \beta_1 & 0 &\dots & \dots &0 \\ \beta_1 & \alpha_2 & \beta_2 &0 &\dots &0 \\ 0& \beta_2 & \alpha_3 & \beta_3 & \dots &0 \\ \dots& \dots & \dots &\dots &\dots & 0 \\ \dots& & &\beta_{n-2} &\alpha_{n-1}& \beta_{n-1} \\ 0& \dots &\dots &0 &\beta_{n-1} & \alpha_{n} \\ \end{array} \right) $$

 
we obtain

 
$$ \hat{A}\hat{q}_k=\beta_{k-1}\hat{q}_{k-1}+\alpha_k\hat{q}_k+\beta_k\hat{q}_{k+1}. $$

Lanczos' iteration

 
$$ \hat{A}\hat{q}_k=\beta_{k-1}\hat{q}_{k-1}+\alpha_k\hat{q}_k+\beta_k\hat{q}_{k+1}, $$

 
with \( \beta_0\hat{q}_0=0 \) for \( k=1:n-1 \). Remember that the vectors \( \hat{q}_k \) are orthornormal and this implies

 
$$ \alpha_k=\hat{q}_k^T\hat{A}\hat{q}_k, $$

 
and these vectors are called Lanczos vectors.

Lanczos' iteration, the algorithm

 
$$ \hat{A}\hat{q}_k=\beta_{k-1}\hat{q}_{k-1}+\alpha_k\hat{q}_k+\beta_k\hat{q}_{k+1}, $$

 
with \( \beta_0\hat{q}_0=0 \) for \( k=1:n-1 \) and

 
$$ \alpha_k=\hat{q}_k^T\hat{A}\hat{q}_k. $$

 
If

 
$$ \hat{r}_k=(\hat{A}-\alpha_k\hat{I})\hat{q}_k-\beta_{k-1}\hat{q}_{k-1}, $$

 
is non-zero, then

 
$$ \hat{q}_{k+1}=\hat{r}_{k}/\beta_k, $$

 
with \( \beta_k=\pm ||\hat{r}_{k}||_2 \).

A simple implementation of the Lanczos algorithm

  r_0 = q_1; beta_0=1; q_0=0; int k = 0;
  while (beta_k != 0)
      q_{k+1} = r_k/beta_k
      k = k+1
      alpha_k = q_k^T A q_k
      r_k = (A-alpha_k I) q_k  -beta_{k-1}q_{k-1}
      beta_k = || r_k||_2
  end while

Differential equations program

• Ordinary differential equations, Runge-Kutta method,chapter 8
• Ordinary differential equations with boundary conditions: one-variable equations to be solved by shooting and Green's function methods, chapter 9
• We can solve such equations by a finite difference scheme as well, turning the equation into an eigenvalue problem. Still one variable. Done in projects 1 and 2.
• If we have more than one variable, we need to solve partial differential equations, see Chapter 10
• Fourier transforms and Fast Fourier transforms if we get time. Project 3 deals with ordinary differential equations (most likely the solar system).

Differential Equations, chapter 8

 
$$ \begin{equation} \frac{dy}{dt}=f(t,y). \end{equation} $$

 
This equation is of first order and \( f \) is an arbitrary function. A second-order equation goes typically like

 
$$ \begin{equation} \frac{d^2y}{dt^2}=f(t,\frac{dy}{dt},y). \end{equation} $$

 
A well-known second-order equation is Newton's second law

 
$$ \begin{equation} m\frac{d^2x}{dt^2}=-kx, \tag{6} \end{equation} $$

 
where \( k \) is the force constant. ODE depend only on one variable

Differential Equations

 
$$ \begin{equation} i\hbar\frac{\partial \psi({\bf x},t)}{\partial t}= -\frac{\hbar^2}{2m}\left( \frac{\partial^2 \psi({\bf r},t)}{\partial x^2} + \frac{\partial^2 \psi({\bf r},t)}{\partial y^2}+ \frac{\partial^2 \psi({\bf r},t)}{\partial z^2}\right) + V({\bf x})\psi({\bf x},t), \end{equation} $$

 
may depend on several variables. In certain cases, like the above equation, the wave function can be factorized in functions of the separate variables, so that the Schr\"odinger equation can be rewritten in terms of sets of ordinary differential equations. These equations are discussed in chapter 10. Involve boundary conditions in addition to initial conditions.

Differential Equations

 
$$ \begin{equation} \frac{dy}{dt}=g^3(t)y(t), \end{equation} $$

 
is an example of a linear equation, while

 
$$ \begin{equation} \frac{dy}{dt}=g^3(t)y(t)-g(t)y^2(t), \end{equation} $$

 
is a non-linear ODE.

Differential Equations

 
$$ \frac{dm}{dr}=4\pi r^{2}\rho (r)/c^2, $$

 
and

 
$$ \frac{dP}{dr}=-\frac{Gm(r)}{r^{2}}\rho (r)/c^2. $$

 
where \( \rho \) is the mass-energy density. The initial conditions are dictated by the mass being zero at the center of the star, i.e., when \( r=0 \), yielding \( m(r=0)=0 \). The other condition is that the pressure vanishes at the surface of the star.

In the solution of the Schr\"odinger equation for a particle in a potential, we may need to apply boundary conditions as well, such as demanding continuity of the wave function and its derivative.

Differential Equations

 
$$ \begin{equation} \frac{dy^{(1)}(t)}{dt}=\frac{dx(t)}{dt}=y^{(2)}(t), \end{equation} $$

 
we can rewrite Newton's second law as two coupled first-order differential equations

 
$$ \begin{equation} m\frac{dy^{(2)}(t)}{dt}=-kx(t)=-ky^{(1)}(t), \tag{7} \end{equation} $$

 
and

 
$$ \begin{equation} \frac{dy^{(1)}(t)}{dt}=y^{(2)}(t). \tag{8} \end{equation} $$

Differential Equations, Finite Difference

 
$$ \begin{equation} y_0=y(t=t_0). \end{equation} $$

 
We are interested in solving a differential equation in a region in space [a,b]. We define a step \( h \) by splitting the interval in \( N \) sub intervals, so that we have

 
$$ \begin{equation} h=\frac{b-a}{N}. \end{equation} $$

 
With this step and the derivative of \( y \) we can construct the next value of the function \( y \) at

 
$$ \begin{equation} y_1=y(t_1=t_0+h), \end{equation} $$

 
and so forth.

Differential Equations

 
$$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h\Delta(t_i,y_i(t_i)) + O(h^{p+1}), \end{equation} $$

 
where \( O(h^{p+1}) \) represents the truncation error. To determine \( \Delta \), we Taylor expand our function \( y \)

 
$$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h(y'(t_i)+\dots +y^{(p)}(t_i)\frac{h^{p-1}}{p!}) + O(h^{p+1}), \tag{9} \end{equation} $$

 
where we will associate the derivatives in the parenthesis with

 
$$ \begin{equation} \Delta(t_i,y_i(t_i))=(y'(t_i)+\dots +y^{(p)}(t_i)\frac{h^{p-1}}{p!}). \tag{10} \end{equation} $$

Differential Equations

 
$$ \begin{equation} y'(t_i)=f(t_i,y_i) \end{equation} $$

 
and if we truncate \( \Delta \) at the first derivative, we have

 
$$ \begin{equation} y_{i+1}=y(t_i) + hf(t_i,y_i) + O(h^2), \tag{11} \end{equation} $$

 
which when complemented with \( t_{i+1}=t_i+h \) forms the algorithm for the well-known Euler method. Note that at every step we make an approximation error of the order of \( O(h^2) \), however the total error is the sum over all steps \( N=(b-a)/h \), yielding thus a global error which goes like \( NO(h^2)\approx O(h) \).

Differential Equations

 
$$ f'_{2c}(x)= \frac{f(x+h)-f(x)}{h}+O(h), $$

 
we can enter into roundoff error problems when we subtract two almost equal numbers \( f(x+h)-f(x)\approx 0 \). Euler's method is not recommended for precision calculation, although it is handy to use in order to get a first view on how a solution may look like. As an example, consider Newton's equation rewritten in Eqs. (7) and (8). We define \( y_0=y^{(1)}(t=0) \) an \( v_0=y^{(2)}(t=0) \). The first steps in Newton's equations are then

 
$$ \begin{equation} y^{(1)}_1=y_0+hv_0+O(h^2) \end{equation} $$

 
and

 
$$ \begin{equation} y^{(2)}_1=v_0-hy_0k/m+O(h^2). \end{equation} $$

Differential Equations

 
$$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+h y^{(2)}_{n+1}+O(h^2) \end{equation} $$

 
and

 
$$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n}+O(h^2). \end{equation} $$

 
The acceleration \( a_n \) is a function of \( a_n(y^{(1)}_{n}, y^{(2)}_{n},t) \) and needs to be evaluated as well. This is the Euler-Cromer method.

Differential Equations

 
$$ \begin{equation} \Delta(t_i,y_i(t_i))=f(t_i)+\frac{h}{2}\frac{df(t_i,y_i)}{dt}+O(h^3). \end{equation} $$

 
The second derivative can be rewritten as

 
$$ \begin{equation} y''=f'=\frac{df}{dt}=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}\frac{\partial y}{\partial t}=\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}f \tag{12} \end{equation} $$

 
and we can rewrite Eq. (9) as

 
$$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) +hf(t_i)+ \frac{h^2}{2}\left(\frac{\partial f}{\partial t}+\frac{\partial f}{\partial y}f\right) + O(h^{3 }), \end{equation} $$

 
which has a local approximation error \( O(h^{3 }) \) and a global error \( O(h^{2}) \).

Differential Equations

 
$$ \begin{equation} y_{i+1}=y(t=t_i+h)=y(t_i) + h(f(t_i,y_i)+\dots f^{(p-1)}(t_i,y_i) \frac{h^{p-1}}{p!}) + O(h^{p+1}). \end{equation} $$

 
These methods, based on higher-order derivatives, are in general not used in numerical computation, since they rely on evaluating derivatives several times. Unless one has analytical expressions for these, the risk of roundoff errors is large.

Differential Equations

 
$$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+\frac{h}{2}\left(y^{(2)}_{n+1}+y^{(2)}_{n}\right)+O(h^2) \end{equation} $$

 
and

 
$$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n}+O(h^2), \end{equation} $$

 
yielding

 
$$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n}+\frac{h^2}{2}a_n+O(h^3) \end{equation} $$

 
implying that the local truncation error in the position is now \( O(h^3) \), whereas Euler's or Euler-Cromer's methods have a local error of \( O(h^2) \).

Differential Equations

One method that avoids this is the so-called half-step method. Here we define

 
$$ \begin{equation} y^{(2)}_{n+1/2}=y^{(2)}_{n-1/2}+h a_{n}+O(h^2), \end{equation} $$

 
and

 
$$ \begin{equation} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n+1/2} +O(h^2). \end{equation} $$

 
Note that this method needs the calculation of \( y^{(2)}_{1/2} \). This is done using e.g., Euler's method

 
$$ \begin{equation} y^{(2)}_{1/2}=y^{(2)}_{0}+h a_{0}+O(h^2). \end{equation} $$

 
As this method is numerically stable, it is often used instead of Euler's method.

Differential Equations

 
$$ \begin{equation} y^{(2)}_{n+1}=y^{(2)}_{n}+h a_{n+1/2}+O(h^2), \tag{13} \end{equation} $$

 
and

 
$$ \begin{equation} \tag{14} y^{(1)}_{n+1}=y^{(1)}_{n}+hy^{(2)}_{n+1/2} +O(h^2). \end{equation} $$

 
The program program2.cpp includes all of the above methods.

Overview of week 40

Overview of week 40

Differential Equations, Runge-Kutta methods

To see this, consider first the following definitions

 
$$ \begin{equation} \frac{dy}{dt}=f(t,y), \end{equation} $$

 
and

 
$$ \begin{equation} y(t)=\int f(t,y) dt, \end{equation} $$

 
and

 
$$ \begin{equation} y_{i+1}=y_i+ \int_{t_i}^{t_{i+1}} f(t,y) dt. \end{equation} $$

Differential Equations, Runge-Kutta methods

 
$$ \begin{equation} \int_{t_i}^{t_{i+1}} f(t,y) dt \approx hf(t_{i+1/2},y_{i+1/2}) +O(h^3). \end{equation} $$

 
This means in turn that we have

 
$$ \begin{equation} y_{i+1}=y_i + hf(t_{i+1/2},y_{i+1/2}) +O(h^3). \end{equation} $$

Differential Equations, Runge-Kutta methods

 
$$ \begin{equation} y_{(i+1/2)}=y_i + \frac{h}{2}\frac{dy}{dt} = y(t_i) + \frac{h}{2}f(t_i,y_i). \end{equation} $$

 
This means that we can define the following algorithm for the second-order Runge-Kutta method, RK2.

 
$$ \begin{equation} k_1=hf(t_i,y_i), \end{equation} $$

 
$$ \begin{equation} k_2=hf(t_{i+1/2},y_i+k_1/2), \end{equation} $$

 
with the final value

 
$$ \begin{equation} y_{i+i}\approx y_i + k_2 +O(h^3). \end{equation} $$

Differential Equations, Runge-Kutta methods

Differential Equations, Runge-Kutta methods

 
$$ \begin{equation} k_1=hf(t_i,y_i), \end{equation} $$

 
$$ \begin{equation} k_2=hf(t_i+h/2,y_i+k_1/2), \end{equation} $$

 
$$ \begin{equation} k_3=hf(t_i+h/2,y_i+k_2/2) \end{equation} $$

 
$$ \begin{equation} k_4=hf(t_i+h,y_i+k_3) \end{equation} $$

 
with the final value

 
$$ \begin{equation} y_{i+1}=y_i +\frac{1}{6}\left( k_1 +2k_2+2k_3+k_4\right). \end{equation} $$

 
Thus, the algorithm consists in first calculating \( k_1 \) with \( t_i \), \( y_1 \) and \( f \) as inputs. Thereafter, we increase the step size by \( h/2 \) and calculate \( k_2 \), then \( k_3 \) and finally \( k_4 \). Global error as \( O(h^4) \).

Simple Example, Block tied to a Wall

 
$$ F=-kx. $$

Simple Example, Block tied to a Wall

 
$$ m\frac{d^2x}{dt^2}=-kx, $$

 
or we could rephrase it as

 
$$ \frac{d^2x}{dt^2}=-\frac{k}{m}x=-\omega_0^2x, \tag{15} $$

 
with the angular frequency \( \omega_0^2=k/m \).

The above differential equation has the advantage that it can be solved analytically with solutions on the form

 
$$ x(t)=Acos(\omega_0t+\nu), $$

 
where \( A \) is the amplitude and \( \nu \) the phase constant. This provides in turn an important test for the numerical solution and the development of a program for more complicated cases which cannot be solved analytically.

Simple Example, Block tied to a Wall

 
$$ \frac{dx(t)}{dt}=v(t), $$

 
and

 
$$ \frac{dv(t)}{dt}=-\omega_0^2x(t). $$

We are now going to solve these equations using the Runge-Kutta method to fourth order discussed previously.

Simple Example, Block tied to a Wall

Since functions like \( cos \) are periodic with a period \( 2\pi \), then the solution \( x(t) \) has also to be periodic. This means that

 
$$ x(t+T)=x(t), $$

 
with \( T \) the period defined as

 
$$ T=\frac{2\pi}{\omega_0}=\frac{2\pi}{\sqrt{k/m}}. $$

Observe that \( T \) depends only on \( k/m \) and not on the amplitude of the solution.

Simple Example, Block tied to a Wall

Suppose we choose the initial conditions

 
$$ x(t=0)=1\quad \mbox{m}\quad v(t=0)=0\quad\mbox{m/s}, $$

 
meaning that block is at rest at \( t=0 \) but with a potential energy

 
$$ E_0=\frac{1}{2}kx(t=0)^2=\frac{1}{2}k. $$

 
The total energy at any time \( t \) has however to be conserved, meaning that our solution has to fulfil the condition

 
$$ E_0=\frac{1}{2}kx(t)^2+\frac{1}{2}mv(t)^2. $$

Simple Example, Block tied to a Wall

• Choose the initial position and speed, with the most common choice \( v(t=0)=0 \) and some fixed value for the position.
• Choose the method you wish to employ in solving the problem.
• Subdivide the time interval \( [t_i,t_f] \) into a grid with step size \( h=\frac{t_f-t_i}{N} \), where \( N \) is the number of mesh points.
• Calculate now the total energy given by \( E_0=\frac{1}{2}kx(t=0)^2=\frac{1}{2}k \).
• The Runge-Kutta method is used to obtain \( x_{i+1} \) and \( v_{i+1} \) starting from the previous values \( x_i \) and \( v_i \)..
• When we have computed \( x(v)_{i+1} \) we upgrade \( t_{i+1}=t_i+h \).
• This iterative process continues till we reach the maximum time \( t_f \).
• The results are checked against the exact solution. Furthermore, one has to check the stability of the numerical solution against the chosen number of mesh points \( N \).

Simple Example, Block tied to a Wall

    y[0] = initial_x;                 // initial position
    y[1] = initial_v;                // initial velocity
    t=0.;                             // initial time
    E0 = 0.5*y[0]*y[0]+0.5*y[1]*y[1];  // the initial total energy
    // now we start solving the differential
    // equations using the RK4 method
    while (t <= tmax){
      derivatives(t, y, dydt);   // initial derivatives
      runge_kutta_4(y, dydt, n, t, h, yout, derivatives);
      for (i = 0; i < n; i++) {
	   y[i] = yout[i];

      t += h;
      output(t, y, E0);   // write to file

Simple Example, Block tied to a Wall

  //   this function sets up the derivatives for this special case
  void derivatives(double t, double *y, double *dydt)
  {
    dydt[0]=y[1];    // derivative of x
    dydt[1]=-y[0]; // derivative of v
  } // end of function derivatives

Runge-Kutta methods, code

void runge_kutta_4(double *y, double *dydx, int n,
                 double x, double h,
          double *yout, void (*derivs)(double, double *, double *))
{
  int i;
  double      xh,hh,h6;
  double *dym, *dyt, *yt;
  //   allocate space for local vectors
  dym = new double [n];
  dyt =  new double [n];
  yt =  new double [n];
  hh = h*0.5;
  h6 = h/6.;
  xh = x+hh;

Runge-Kutta methods, code

  for (i = 0; i < n; i++) {
    yt[i] = y[i]+hh*dydx[i];

  (*derivs)(xh,yt,dyt);     // computation of k2
  for (i = 0; i < n; i++) {
    yt[i] = y[i]+hh*dyt[i];

  (*derivs)(xh,yt,dym); //  computation of k3
  for (i=0; i < n; i++) {
    yt[i] = y[i]+h*dym[i];
    dym[i] += dyt[i];

  (*derivs)(x+h,yt,dyt);    // computation of k4
  //      now we upgrade y in the array yout
  for (i = 0; i < n; i++){
    yout[i] = y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]);

  delete []dym;
  delete [] dyt;
  delete [] yt;
}       //  end of function Runge-kutta 4

The classical pendulum

 
$$ \begin{equation} ml\frac{d^2\theta}{dt^2}+mgsin(\theta)=0, \end{equation} $$

 
with an angular velocity and acceleration given by

 
$$ \begin{equation} v=l\frac{d\theta}{dt}, \end{equation} $$

 
and

 
$$ \begin{equation} a=l\frac{d^2\theta}{dt^2}. \end{equation} $$

More on the Pendulum

 
$$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=0, \tag{16} \end{equation} $$

 
where \( \nu \) is now a positive constant parameterizing the viscosity of the medium in question. In order to maintain the motion against viscosity, it is necessary to add some external driving force. We choose here a periodic driving force. The last equation becomes then

 
$$ \begin{equation} ml\frac{d^2\theta}{dt^2}+\nu\frac{d\theta}{dt} +mgsin(\theta)=Asin(\omega t), \tag{17} \end{equation} $$

 
with \( A \) and \( \omega \) two constants representing the amplitude and the angular frequency respectively. The latter is called the driving frequency.

More on the Pendulum

 
$$ \omega_0=\sqrt{g/l}, $$

 
the so-called natural frequency and the new dimensionless quantities

 
$$ \hat{t}=\omega_0t, $$

 
with the dimensionless driving frequency

 
$$ \hat{\omega}=\frac{\omega}{\omega_0}, $$

 
and introducing the quantity \( Q \), called the quality factor,

 
$$ Q=\frac{mg}{\omega_0\nu}, $$

 
and the dimensionless amplitude

 
$$ \hat{A}=\frac{A}{mg} $$

More on the Pendulum

 
$$ \frac{d^2\theta}{d\hat{t}^2}+\frac{1}{Q}\frac{d\theta}{d\hat{t}} +sin(\theta)=\hat{A}cos(\hat{\omega}\hat{t}). $$

This equation can in turn be recast in terms of two coupled first-order differential equations as follows

 
$$ \frac{d\theta}{d\hat{t}}=\hat{v}, $$

 
and

 
$$ \frac{d\hat{v}}{d\hat{t}}=-\frac{\hat{v}}{Q}-sin(\theta)+\hat{A}cos(\hat{\omega}\hat{t}). $$

These are the equations to be solved. The factor \( Q \) represents the number of oscillations of the undriven system that must occur before its energy is significantly reduced due to the viscous drag. The amplitude \( \hat{A} \) is measured in units of the maximum possible gravitational torque while \( \hat{\omega} \) is the angular frequency of the external torque measured in units of the pendulum's natural frequency.

Classes for ODE methods

Classes for ODE methods

void euler();
void euler_cromer();
void midpoint();
void euler_richardson();
void half_step();
void rk2(); //runge-kutta-second-order
void rk4_step(double,double*,double*,double); // we need it in function rk4() and asc()
void rk4(); //runge-kutta-fourth-order
void asc(); //runge-kutta-fourth-order with adaptive stepsize control

Classes for ODE methods

class pendulum
 {
 private:
   double Q, A_roof, omega_0, omega_roof,g; //
   double y[2];          //for the initial-values of phi and v
   int n;                // how many steps
   double delta_t,delta_t_roof;

 public:
   void derivatives(double,double*,double*);
   void initialise();
   void euler();
   void euler_cromer();
   void midpoint();
   void euler_richardson();
   void half_step();
   void rk2(); //runge-kutta-second-order
   void rk4_step(double,double*,double*,double); // we need it in function rk4() and asc()
   void rk4(); //runge-kutta-fourth-order
   void asc(); //runge-kutta-fourth-order with adaptive stepsize control
 };

Classes for ODE methods

void pendulum::derivatives(double t, double* in, double* out)
{ /* Here we are calculating the derivatives at (dimensionless) time t
     'in' are the values of phi and v, which are used for the calculation
     The results are given to 'out' */

  out[0]=in[1];             //out[0] = (phi)'  = v
  if(Q)
    out[1]=-in[1]/((double)Q)-sin(in[0])+A_roof*cos(omega_roof*t);  //out[1] = (phi)''
  else
    out[1]=-sin(in[0])+A_roof*cos(omega_roof*t);  //out[1] = (phi)''

Classes for ODE methods

int main()
{
  pendulum testcase;
  testcase.initialise();
  testcase.euler();
  testcase.euler_cromer();
  testcase.midpoint();
  testcase.euler_richardson();
  testcase.half_step();
  testcase.rk2();
  testcase.rk4();
  return 0;
}  // end of main function

Classes for ODE methods

MODULE pendulum
   USE CONSTANTS
   IMPLICIT NONE
   REAL(DP), PRIVATE :: Q, A_roof, omega_0, omega_roof,g
   REAL(DP), PRIVATE :: y(2)         ! for the initial-values of phi and v
   INTEGER, PRIVATE ::  n               ! how many steps
   REAL(DP), PRIVATE :: delta_t,delta_t_roof

   CONTAINS
    SUBROUTINE derivatives(..)
    SUBROUTINE initialise(..)
    SUBROUTINE euler(..)
    SUBROUTINE euler_cromer(..)
    SUBROUTINE midpoint(..)
    etc

END MODULE pendulum

The report: how to write a good scienfitic/technical report

The report

The report

The report

The report

The report

The report

Overview of week 41

Overview of week 41

Student project 3: An object oriented example program for project 3

Adaptive methods

Adaptive methods

 
$$ \tilde{x}=x_1+Ch^{M+1}+O(h^{M+2}), $$

 
and

 
$$ \tilde{x}=x_2+2C(h/2)^{M+1}+O(h^{M+2}), $$

 
with \( C \) a constant. Note that we calculate two halves in the last equation. We get then

 
$$ |x_1-x_2| = Ch^{M+1}(1-\frac{1}{2^M}). $$

 
yielding

 
$$ C=\frac{|x_1-x_2|}{(1-2^{-M})h^{M+1}}. $$

 
We rewrite

 
$$ \tilde{x}=x_2+\epsilon+O((h)^{M+2}), $$

 
with

 
$$ \epsilon = \frac{|x_1-x_2|}{2^M-1}. $$

Adaptive methods

 
$$ \tilde{x}=x_2+\epsilon+O((h)^{6}), $$

 
with

 
$$ \epsilon = \frac{|x_1-x_2|}{15}. $$

 
The estimate is one order higher than the original RK4. But this method is normally rather inefficient since it requires a lot of computations. We solve typically the equation three times at each time step. However, we can compare the estimate \( \epsilon \) with some by us given accuracy \( \xi \). We can then ask the question: what is, with a given \( x_j \) and \( t_j \), the largest possible step size \( \tilde{h} \) that leads to a truncation error below \( \xi \)? We want

 
$$ C\tilde{h} \le \xi, $$

 
which leads to

 
$$ \left(\frac{\tilde{h}}{h}\right)^{M+1}\frac{|x_1-x_2|}{(1-2^{-M})}\le \xi, $$

 
meaning that

 
$$ \tilde{h}=h\left(\frac{\xi}{\epsilon}\right)^{1+1/M}. $$

Adaptive methods

 
$$ \tilde{h}=h\left(\frac{\xi}{\epsilon}\right)^{1+1/M}. $$

 
we can design the following algorithm:

• If the two answers are close, keep the approximation to \( h \).
• If \( \epsilon > \xi \) we need to decrease the step size in the next time step.
• If \( \epsilon < \xi \) we need to increase the step size in the next time step. A much used algorithm is the so-called RKF45 which uses a combination of a fourth and fifth order RK methods.

Adaptive methods, RKF45

 
$$ k_1 = h f (t_k , y_k ), $$

 
$$ k_2 = h f (t_k + \frac{1}{4}h, y_k + \frac{1}{4}k_1) , $$

 
$$ k_3 = h f (t_k + \frac{3}{8}h, y_k + \frac{3}{32}k_1 + \frac{9}{32}k_2) , $$

 
$$ k_4 = h f (t_k + \frac{12}{13}h, y_k + \frac{1932}{2197}k_1 + \frac{7200}{2197}k_2+\frac{7296}{2197}k_3), $$

 
$$ k_5 = h f (t_k + h, y_k + \frac{439}{216}k_1 -8k_2+ \frac{3680}{513}k_3+\frac{845}{4104}k_4), $$

 
$$ k_6 = h f (t_k + \frac{1}{2}h, y_k - \frac{8}{27}k_1 + 2k_2-\frac{3544}{2565}k_2+\frac{1859}{4104}k_4-+\frac{11}{40}k_5). $$

Adaptive methods, RKF45

 
$$ y_{k+1} = y_k + \frac{25}{216}k_1+\frac{1408}{2565}k_3 +\frac{2197}{4101}k_4-\frac{1}{5}k_5, $$

 
where the four function values \( k_1 \) , \( k_3 \) , \( k_4 \) , and \( k_5 \) are used. Notice that \( k_2 \) is not used here. A better value for the solution is determined using a Runge-Kutta method of order 5:

 
$$ z_{k+1} = y_k + \frac{16}{135}k_1+\frac{6656}{12825}k_3 +\frac{28561}{56430}k_4-\frac{9}{50}k_5+\frac{2}{55}k_6. $$

The optimal time step \( \alpha h \) is then determined by

 
$$ \alpha = \left( \frac{\xi h}{2|z_{k+1}-y_{k+1}|}\right)^{1/4}, $$

 
with \( \xi \) our defined tolerance.

Partial Differential Equations, chapter 10

 
$$ A(x,y)\frac{\partial^2 U}{\partial x^2}+B(x,y)\frac{\partial^2 U}{\partial x\partial y} +C(x,y)\frac{\partial^2 U}{\partial y^2}=F(x,y,U,\frac{\partial U}{\partial x}, \frac{\partial U}{\partial y}) $$

 
Examples

 
$$ B=C=0, $$

 
give e.g., 1+1-dim diffusion equation

 
$$ A\frac{\partial^2 U}{\partial x^2}=\frac{\partial U}{\partial t} $$

 
and is an example of a parabolic PDE

Partial Differential Equations

 
$$ A\frac{\partial^2 U}{\partial x^2}+C\frac{\partial^2 U}{\partial y^2}=\frac{\partial^2 U}{\partial t^2} $$

 
Poisson's (Laplace's \( \rho =0 \)) equation

 
$$ \nabla^2 U({\bf x})=-4\pi \rho({\bf x}). $$

Heat/Diffusion Equation

 
$$ \frac{\kappa}{C\rho}\nabla^2 T({\bf x},t) =\frac{\partial T({\bf x},t)}{\partial t} $$

 
$$ \frac{\kappa}{C\rho({\bf x},t)}\nabla^2 T({\bf x},t) =\frac{\partial T({\bf x},t)}{\partial t} $$

Explicit Scheme for the Diffusion Equation

 
$$ \nabla^2 u(x,t) =\frac{\partial u(x,t)}{\partial t}, $$

 
or

 
$$ u_{xx} = u_t, $$

 
with initial conditions, i.e., the conditions at \( t=0 \),

 
$$ u(x,0)= g(x) \quad 0 \le x \le L $$

 
with \( L=1 \) the length of the \( x \)-region of interest. The boundary conditions are

 
$$ u(0,t)= a(t) \quad t \ge 0, $$

 
and

 
$$ u(L,t)= b(t) \quad t \ge 0, $$

 
where \( a(t) \) and \( b(t) \) are two functions which depend on time only, while \( g(x) \) depends only on the position \( x \).

Explicit Scheme, Forward Euler

 
$$ u_t\approx \frac{u_{i,j+1}-u_{i,j}}{\Delta t}, $$

 
and

 
$$ u_{xx}\approx \frac{u_{i+i,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2}. $$

 
The one-dimensional diffusion equation can then be rewritten in its discretized version as

 
$$ \frac{u_{i,j+1}-u_{i,j}}{\Delta t}=\frac{u_{i+i,j}-2u_{i,j}+u_{i-1,j}}{\Delta x^2}. $$

 
Defining \( \alpha = \Delta t/\Delta x^2 \) results in the explicit scheme

 
$$ \tag{18} u_{i,j+1}= \alpha u_{i-1,j}+(1-2\alpha)u_{i,j}+\alpha u_{i+1,j}. $$

Explicit Scheme

 
$$ V_{j+1} = AV_{j} $$

 
with

 
$$ A=\left(\begin{array}{cccc}1-2\alpha&\alpha&0& 0\dots\\ \alpha&1-2\alpha&\alpha & 0\dots \\ \dots & \dots & \dots & \dots \\ 0\dots & 0\dots &\alpha& 1-2\alpha\end{array} \right) $$

 
yielding

 
$$ V_{j+1} = AV_{j}=\dots = A^jV_0 $$

 
The explicit scheme, although being rather simple to implement has a very weak stability condition given by

 
$$ \Delta t/\Delta x^2 \le 1/2 $$

Implicit Scheme

 
$$ u_t\approx \frac{u(x_i,t_j)-u(x_i,t_j-k)}{k} $$

 
and

 
$$ u_{xx}\approx \frac{u(x_i+h,t_j)-2u(x_i,t_j)+u(x_i-h,t_j)}{h^2} $$

 
Define \( \alpha = k/h^2 \). Gives

 
$$ u_{i,j-1}= -\alpha u_{i-1,j}+(1-2\alpha)u_{i,j}-\alpha u_{i+1,j} $$

 
Here \( u_{i,j-1} \) is the only unknown quantity.

Have

 
$$ AV_{j} = V_{j-1} $$

 
with

 
$$ A=\left(\begin{array}{cccc}1+2\alpha&-\alpha&0& 0\dots\\ -\alpha&1+2\alpha&-\alpha & 0\dots \\ \dots & \dots & \dots & \dots \\ 0\dots & 0\dots &-\alpha& 1+2\alpha\end{array} \right) $$

 
which gives

 
$$ V_{j} = A^{-1}V_{j-1}=\dots = A^{-j}V_0 $$

 
Need only to invert a matrix

Brute Force Implicit Scheme, inefficient algo

!  now invert the matrix
       CALL matinv( a, ndim, det)
       DO i = 1, m
          DO l=1, ndim
             u(l) = DOT_PRODUCT(a(l,:),v(:))
          ENDDO
          v = u
          t = i*k
          DO  j=1, ndim
              WRITE(6,*) t, j*h, v(j)
          ENDDO
       ENDDO

Brief Summary of the Explicit and the Implicit Methods

 
$$ u_t\approx \frac{u(x,t+\Delta t)-u(x,t)}{\Delta t}=\frac{u(x_i,t_j+\Delta t)-u(x_i,t_j)}{\Delta t} $$

The implicit method can be applied in a brute force way as well as long as the element of the matrix are constants.

 
$$ u_t\approx \frac{u(x,t)-u(x,t-\Delta t)}{\Delta t}=\frac{u(x_i,t_j)-u(x_i,t_j-\Delta t)}{\Delta t} $$

However, it is more efficient to use a linear algebra solver for tridiagonal matrices.

Crank-Nicolson

 
$$ \frac{\theta}{\Delta x^2}\left(u_{i-1,j}-2u_{i,j}+u_{i+1,j}\right)+ \frac{1-\theta}{\Delta x^2}\left(u_{i+1,j-1}-2u_{i,j-1}+u_{i-1,j-1}\right)= \frac{1}{\Delta t}\left(u_{i,j}-u_{i,j-1}\right), $$

 
which for \( \theta=0 \) yields the forward formula for the first derivative and the explicit scheme, while \( \theta=1 \) yields the backward formula and the implicit scheme. These two schemes are called the backward and forward Euler schemes, respectively. For \( \theta = 1/2 \) we obtain a new scheme after its inventors, Crank and Nicolson.

Crank Nicolson

 
$$ -\alpha u_{i-1,j}+\left(2+2\alpha\right)u_{i,j}-\alpha u_{i+1,j}= \alpha u_{i-1,j-1}+\left(2-2\alpha\right)u_{i,j-1}+\alpha u_{i+1,j-1}, $$

 
or in matrix-vector form as

 
$$ \left(2\hat{I}+\alpha\hat{B}\right)V_{j}= \left(2\hat{I}-\alpha\hat{B}\right)V_{j-1}, $$

 
where the vector \( V_{j} \) is the same as defined in the implicit case while the matrix \( \hat{B} \) is

 
$$ \hat{B}=\left(\begin{array}{cccc}2&-1&0& 0\dots\\ -1&2&-1 & 0\dots \\ \dots & \dots & \dots & \dots \\ 0\dots & 0\dots && 2\end{array} \right) $$

Analysis of diffusion equation

 
$$ \begin{align} u(x+\Delta x,t)&=u(x,t)+\frac{\partial u(x,t)}{\partial x} \Delta x+\frac{\partial^2 u(x,t)}{2\partial x^2}\Delta x^2+\mathcal{O}(\Delta x^3), \nonumber\\ u(x-\Delta x,t)&=u(x,t)-\frac{\partial u(x,t)}{\partial x}\Delta x+\frac{\partial^2 u(x,t)}{2\partial x^2} \Delta x^2+\mathcal{O}(\Delta x^3), \nonumber\\ u(x,t+\Delta t)&=u(x,t)+\frac{\partial u(x,t)}{\partial t}\Delta t+ \mathcal{O}(\Delta t^2). \tag{19} \end{align} $$

Analysis of diffusion equation

 
$$ \begin{align} &\left[\frac{\partial u(x,t)}{\partial t}\right]_{\text{approx}} =\frac{\partial u(x,t)}{\partial t}+\mathcal{O}(\Delta t) , \nonumber\\ &\left[\frac{\partial^2 u(x,t)}{\partial x^2}\right]_{\text{approx}} =\frac{\partial^2 u(x,t)}{\partial x^2}+\mathcal{O}(\Delta x^2). \tag{20} \end{align} $$

 
It is easy to convince oneself that the backward Euler method must have the same truncation errors as the forward Euler scheme.

Analysis of diffusion equation

 
$$ \begin{align} u(x+\Delta x, t+\Delta t)&=u(x,t')+\frac{\partial u(x,t')}{\partial x}\Delta x+\frac{\partial u(x,t')}{\partial t} \frac{\Delta t}{2} +\frac{\partial^2 u(x,t')}{2\partial x^2}\Delta x^2+\frac{\partial^2 u(x,t')}{2\partial t^2}\frac{\Delta t^2}{4} +\notag \\ \nonumber &\frac{\partial^2 u(x,t')}{\partial x\partial t}\frac{\Delta t}{2} \Delta x+ \mathcal{O}(\Delta t^3)\\ \nonumber u(x-\Delta x, t+\Delta t)&=u(x,t')-\frac{\partial u(x,t')}{\partial x}\Delta x+\frac{\partial u(x,t')}{\partial t} \frac{\Delta t}{2} +\frac{\partial^2 u(x,t')}{2\partial x^2}\Delta x^2+\frac{\partial^2 u(x,t')}{2\partial t^2}\frac{\Delta t^2}{4} - \notag\\ \nonumber &\frac{\partial^2 u(x,t')}{\partial x\partial t}\frac{\Delta t}{2} \Delta x+ \mathcal{O}(\Delta t^3)\\ u(x+\Delta x,t)&=u(x,t')+\frac{\partial u(x,t')}{\partial x}\Delta x-\frac{\partial u(x,t')}{\partial t} \frac{\Delta t}{2} +\frac{\partial^2 u(x,t')}{2\partial x^2}\Delta x^2+\frac{\partial^2 u(x,t')}{2\partial t^2}\frac{\Delta t^2}{4} -\notag \\ \nonumber &\frac{\partial^2 u(x,t')}{\partial x\partial t}\frac{\Delta t}{2} \Delta x+ \mathcal{O}(\Delta t^3)\\ \nonumber u(x-\Delta x,t)&=u(x,t')-\frac{\partial u(x,t')}{\partial x}\Delta x-\frac{\partial u(x,t')}{\partial t} \frac{\Delta t}{2} +\frac{\partial^2 u(x,t')}{2\partial x^2}\Delta x^2+\frac{\partial^2 u(x,t')}{2\partial t^2}\frac{\Delta t^2}{4} +\notag \\ \nonumber &\frac{\partial^2 u(x,t')}{\partial x\partial t}\frac{\Delta t}{2} \Delta x+ \mathcal{O}(\Delta t^3)\\ \nonumber u(x,t+\Delta t)&=u(x,t')+\frac{\partial u(x,t')}{\partial t}\frac{\Delta_t}{2} +\frac{\partial ^2 u(x,t')}{2\partial t^2}\Delta t^2 + \mathcal{O}(\Delta t^3) \nonumber\\ u(x,t)&=u(x,t')-\frac{\partial u(x,t')}{\partial t}\frac{\Delta t}{2}+\frac{\partial ^2 u(x,t')}{2\partial t^2}\Delta t^2 + \mathcal{O}(\Delta t^3) \tag{21} \end{align} $$

Analysis of diffusion equation

 
$$ \begin{align} &\left[\frac{\partial u(x,t')}{\partial t}\right]_{\text{approx}} =\frac{\partial u(x,t')}{\partial t}+\mathcal{O}(\Delta t^2) , \\ \nonumber &\left[\frac{\partial^2 u(x,t')}{\partial x^2}\right]_{\text{approx}}=\frac{\partial^2 u(x,t')}{\partial x^2}+\mathcal{O}(\Delta x^2). \end{align} $$

Analysis of diffusion equation

Analysis of diffusion equation

 
$$ \nabla^2 u(x,t) =\frac{\partial u(x,t)}{\partial t}, $$

 
with initial conditions

 
$$ u(x,0)= g(x) \quad 0 < x < L. $$

Analysis of diffusion equation

 
$$ u(0,t)= 0 \quad t \ge 0, \quad u(L,t)= 0 \quad t \ge 0, $$

We assume that we have solutions of the form (separation of variable)

 
$$ \begin{equation} u(x,t)=F(x)G(t), \end{equation} $$

 
which inserted in the partial differential equation results in

 
$$ \begin{equation} \frac{F''}{F}=\frac{G'}{G}, \end{equation} $$

 
where the derivative is with respect to \( x \) on the left hand side and with respect to \( t \) on right hand side. This equation should hold for all \( x \) and \( t \). We must require the rhs and lhs to be equal to a constant.

Analysis of diffusion equation

 
$$ \begin{equation} F''+\lambda^2F=0; \quad G'=-\lambda^2G, \end{equation} $$

 
with general solutions

 
$$ \begin{equation} F(x)=A\sin(\lambda x)+B\cos(\lambda x); \quad G(t)=Ce^{-\lambda^2t}. \end{equation} $$

Analysis of diffusion equation

 
$$ \begin{equation} u(x,t)=A_n\sin(n\pi x/L)e^{-n^2\pi^2 t/L^2}. \end{equation} $$

 
But there are infinitely many possible \( n \) values (infinite number of solutions). Moreover, the diffusion equation is linear and because of this we know that a superposition of solutions will also be a solution of the equation. We may therefore write

 
$$ \begin{equation} u(x,t)=\sum_{n=1}^{\infty} A_n \sin(n\pi x/L) e^{-n^2\pi^2 t/L^2}. \end{equation} $$

Analysis of diffusion equation

 
$$ \begin{equation} u(x,0)=g(x)=\sum_{n=1}^{\infty} A_n \sin(n\pi x/L). \end{equation} $$

 
The coefficient \( A_n \) is the Fourier coefficients for the function \( g(x) \). Because of this, \( A_n \) is given by (from the theory on Fourier series)

 
$$ \begin{equation} A_n=\frac{2}{L}\int_0^L g(x)\sin(n\pi x/L) \mbox{d}x. \end{equation} $$

 
Different \( g(x) \) functions will obviously result in different results for \( A_n \).

Overview of week 42

Overview of week 42

Laplace's and Poisson's equations

 
$$ \nabla^2 u({\bf x})=u_{xx}+u_{yy}=0. $$

 
with possible boundary conditions \( u(x,y) = g(x,y) \) on the border. There is no time-dependence. Choosing equally many steps in both directions we have a quadratic or rectangular grid, depending on whether we choose equal steps lengths or not in the \( x \) and the \( y \) directions. Here we set \( \Delta x = \Delta y = h \) and obtain a discretized version

 
$$ u_{xx}\approx \frac{u(x+h,y)-2u(x,y)+u(x-h,y)}{h^2}, $$

 
and

 
$$ u_{yy}\approx \frac{u(x,y+h)-2u(x,y)+u(x,y-h)}{h^2}, $$

Laplace's and Poisson's equations

 
$$ u_{xx}\approx \frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h^2}, $$

 
and

 
$$ u_{yy}\approx \frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h^2}, $$

 
which gives when inserted in Laplace's equation

 
$$ \begin{equation} \tag{22} u_{i,j}= \frac{1}{4}\left[u_{i,j+1}+u_{i,j-1}+u_{i+1,j}+u_{i-1,j}\right]. \end{equation} $$

 
This is our final numerical scheme for solving Laplace's equation. Poisson's equation adds only a minor complication to the above equation since in this case we have

 
$$ u_{xx}+u_{yy}=-\rho({\bf x}), $$

 
and we need only to add a discretized version of \( \rho({\bf x}) \) resulting in

 
$$ \begin{equation} \tag{23} u_{i,j}= \frac{1}{4}\left[u_{i,j+1}+u_{i,j-1}+u_{i+1,j}+u_{i-1,j}\right] +\rho_{i,j}. \end{equation} $$

Solution Approach

Code example for the two-dimensional diff equation/Laplace

int DiffusionJacobi(int N, double dx, double dt,
		      double **A, double **q, double abstol){
  int i,j,k;
  int maxit = 100000;
  double sum;
  double ** Aold = CreateMatrix(N,N);

  double D = dt/(dx*dx);

Code example for the two-dimensional diff equation/Laplace

  for(i=1; i<N-1; i++)
    for(j=1;j<N-1;j++)
      Aold[i][j] = 1.0;
  /* Boundary Conditions -- all zeros */
  for(i=0;i<N;i++){
    A[0][i] = 0.0;
    A[N-1][i] = 0.0;
    A[i][0] = 0.0;
    A[i][N-1] = 0.0;

Code example for the two-dimensional diff equation/Laplace

  for(k=0; k<maxit; k++){
    for(i = 1; i<N-1; i++){
      for(j=1; j<N-1; j++){
	A[i][j] = dt*q[i][j] + Aold[i][j] +
	  D*(Aold[i+1][j] + Aold[i][j+1] - 4.0*Aold[i][j] +
	     Aold[i-1][j] + Aold[i][j-1]);

    sum = 0.0;
    for(i=0;i<N;i++){
      for(j=0;j<N;j++){
	sum += (Aold[i][j]-A[i][j])*(Aold[i][j]-A[i][j]);
	Aold[i][j] = A[i][j];

    if(sqrt(sum)<abstol){DestroyMatrix(Aold,N,N);
      return k;

Two-dimensional wave equation

 
$$ \left\{\begin{array}{cc} \lambda\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right) = \frac{\partial^2u}{\partial t^2}& x,y\in[0,1], t \ge 0 \\ u(x,y,0) = sin(\pi x)sin(2\pi y)& x,y\in (0,1) \\ u = 0 \quad \mbox{boundary} & t \ge 0\\ \partial u/\partial t|_{t=0}=0 & x,y\in (0,1)\\ \end{array}\right. . $$

 
The boundary is defined by \( x=0 \), \( x=1 \), \( y=0 \) and \( y=1 \). Here we set \( \lambda = 1 \).

Two-dimensional wave equation

 
$$ \left\{\begin{array}{cc} t_l=l\Delta t& l \ge 0 \\ x_i=i\Delta x& 0 \le i \le n_x\\ y_j=j\Delta y& 0 \le j \le n_y\end{array}\right. , $$

 
and we will let \( \Delta x=\Delta y = h \) and \( n_x=n_y \) for the sake of simplicity. We have now the following discretized partial derivatives

 
$$ u_{xx}\approx \frac{u_{i+1,j}^l-2u_{i,j}^l+u_{i-1,j}^l}{h^2}, $$

 
and

 
$$ u_{yy}\approx \frac{u_{i,j+1}^l-2u_{i,j}^l+u_{i,j-1}^l}{h^2}, $$

 
and

 
$$ u_{tt}\approx \frac{u_{i,j}^{l+1}-2u_{i,j}^{l}+u_{i,j}^{l-1}}{\Delta t^2}. $$

Two-dimensional wave equation

 
$$ \tag{24} u_{i,j}^{l+1} =2u_{i,j}^{l}-u_{i,j}^{l-1}+\frac{\Delta t^2}{h^2}\left(u_{i+1,j}^l-4u_{i,j}^l+u_{i-1,j}^l+u_{i,j+1}^l+u_{i,j-1}^l\right), $$

 
where again we have an explicit scheme with \( u_{i,j}^{l+1} \) as the only unknown quantity. It is easy to account for different step lengths for \( x \) and \( y \). The partial derivative is treated in much the same way as for the one-dimensional case, except that we now have an additional index due to the extra spatial dimension, viz., we need to compute \( u_{i,j}^{-1} \) through

 
$$ u_{i,j}^{-1}=u_{i,j}^0+\frac{\Delta t}{2h^2}\left(u_{i+1,j}^0-4u_{i,j}^0+u_{i-1,j}^0+u_{i,j+1}^0+u_{i,j-1}^0\right), $$

 
in our setup of the initial conditions.

Code example for the two-dimensional wave equation

  //  After initializations and declaration of variables
  for ( int i = 0; i < n; i++ ) {
    u[i] = new double [n];
    uLast[i] = new double [n];
    uNext[i] = new double [n];
    x[i] = i*h;
    y[i] = x[i];

  // initializing
  for ( int i = 0; i < n; i++ ) {  // setting initial step
    for ( int j = 0; j < n; j++ ) {
      uLast[i][j] = sin(PI*x[i])*sin(2*PI*y[j]);

Code example for the two-dimensional wave equation

  for ( int i = 1; i < (n-1); i++ ) {  // setting first step using the initial derivative
    for ( int j = 1; j < (n-1); j++ ) {
      u[i][j] = uLast[i][j] - ((tStep*tStep)/(2.0*h*h))*
	(4*uLast[i][j] - uLast[i+1][j] - uLast[i-1][j] - uLast[i][j+1] - uLast[i][j-1]);

    u[i][0] = 0;  // setting boundaries once and for all
    u[i][n-1] = 0;
    u[0][i] = 0;
    u[n-1][i] = 0;

    uNext[i][0] = 0;
    uNext[i][n-1] = 0;
    uNext[0][i] = 0;
    uNext[n-1][i] = 0;

Code example for the two-dimensional wave equation

  // iterating in time
  double t = 0.0 + tStep;
  int iter = 0;

  while ( t < tFinal ) {
    iter ++;
    t = t + tStep;

    for ( int i = 1; i < (n-1); i++ ) {  // computing next step
      for ( int j = 1; j < (n-1); j++ ) {
	uNext[i][j] = 2*u[i][j] - uLast[i][j] - ((tStep*tStep)/(h*h))*
	  (4*u[i][j] - u[i+1][j] - u[i-1][j] - u[i][j+1] - u[i][j-1]);

Code example for the two-dimensional wave equation

    for ( int i = 1; i < (n-1); i++ ) {  // shifting results down
      for ( int j = 1; j < (n-1); j++ ) {
	uLast[i][j] = u[i][j];
	u[i][j] = uNext[i][j];

Closed form solution of the wave equation

 
$$ \left\{\begin{array}{cc} c^2(u_{xx}+u_{yy}) = u_{tt}& x,y\in(0,L), t>0 \\ u(x,y,0) = f(x,y) & x,y\in (0,L) \\ u(0,0,t)=u(L,L,t)=0 & t > 0\\ \partial u/\partial t|_{t=0}= g(x,y) & x,y\in (0,L)\\ \end{array}\right. . $$

Closed form solution of the wave equation

 
$$ u(x,y,t) = F(x,y) G(t), $$

 
resulting in the equation

 
$$ FG_{tt}= c^2(F_{xx}G+F_{yy}G), $$

 
or

 
$$ \frac{G_{tt}}{c^2G} = \frac{1}{F}(F_{xx}+F_{yy}) = -\nu^2. $$

 
The lhs and rhs are independent of each other and we obtain two differential equations

 
$$ F_{xx}+F_{yy}+F\nu^2=0, $$

 
and

 
$$ G_{tt} + Gc^2\nu^2 = G_{tt} + G\lambda^2 = 0, $$

 
with \( \lambda = c\nu \).

Closed form solution of the wave equation

 
$$ F(x,y) = H(x)Q(y), $$

 
which results in

 
$$ \frac{1}{H}H_{xx} = -\frac{1}{Q}(Q_{yy}+Q\nu^2)= -\kappa^2. $$

 
Since the lhs and rhs are again independent of each other, we can separate the latter equation into two independent equations, one for \( x \) and one for \( y \), namely

 
$$ H_{xx} + \kappa^2H = 0, $$

 
and

 
$$ Q_{yy} + \rho^2Q = 0, $$

 
with \( \rho^2= \nu^2-\kappa^2 \).

Closed form solution of the wave equation

 
$$ H(x) = A\cos(\kappa x)+B\sin(\kappa x), $$

 
and

 
$$ Q(y) = C\cos(\rho y)+D\sin(\rho y). $$

 
The boundary conditions require that \( F(x,y) = H(x)Q(y) \) are zero at the boundaries, meaning that \( H(0)=H(L)=Q(0)=Q(L)=0 \). This yields the solutions

 
$$ H_m(x) = \sin(\frac{m\pi x}{L}) \quad Q_n(y) = \sin(\frac{n\pi y}{L}), $$

 
or

 
$$ F_{mn}(x,y) = \sin(\frac{m\pi x}{L})\sin(\frac{n\pi y}{L}). $$

Closed form solution of the wave equation

 
$$ G_{mn}(t) = B_{mn}\cos(\lambda_{mn} t)+D_{mn}\sin(\lambda_{mn} t), $$

 
with the general solution of the form

 
$$ u(x,y,t) = \sum_{mn=1}^{\infty} u_{mn}(x,y,t) = \sum_{mn=1}^{\infty}F_{mn}(x,y)G_{mn}(t). $$

Closed form solution of the wave equation

 
$$ B_{mn} = \frac{2}{L}\int_0^L\int_0^L dxdy f(x,y) \sin(\frac{m\pi x}{L})\sin(\frac{n\pi y}{L}), $$

 
and

 
$$ D_{mn} = \frac{2}{L}\int_0^L\int_0^L dxdy g(x,y) \sin(\frac{m\pi x}{L})\sin(\frac{n\pi y}{L}). $$

 
Inserting the particular functional forms of \( f(x,y) \) and \( g(x,y) \) one obtains the final analytic expressions.

Two-dimensional wave equation

 
$$ \Delta t \le \frac{1}{\sqrt{\lambda}}\left(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}\right)^{-1/2} $$

 
where \( \Delta t \), \( \Delta x \) and \( \Delta y \) are the chosen step lengths. In our case \( \Delta x=\Delta y=h \). How do we find this condition? In one dimension we can proceed as we did for the diffusion equation.

Two-dimensional wave equation

 
$$ u(x,y,t) = A \exp{ (i(k_xx+k_yy-\omega t))} $$

 
Then from

 
$$ u_{xx}\approx \frac{u_{i+1,j}^l-2u_{i,j}^l+u_{i-1,j}^l}{\Delta x^2}, $$

 
we get, with \( u_i=\exp{ikx_i} \)

 
$$ u_{xx}\approx \frac{u_i}{\Delta x^2}\left(\exp{ik\Delta x}-2+\exp{(-ik\Delta x)}\right), $$

 
or

 
$$ u_{xx}\approx 2\frac{u_i}{\Delta x^2}\left(cos(k\Delta x)-1\right)= -4\frac{u_i}{\Delta x^2}sin^2(k\Delta x/2) $$

 
We get similar results for \( t \) and \( y \).

Two-dimensional wave equation

 
$$ \lambda\left(\frac{\partial^2u}{\partial x^2}+\frac{\partial^2u}{\partial y^2}\right) = \frac{\partial^2u}{\partial t^2}, $$

 
resulting in

 
$$ \lambda\left(-4\frac{u_{ij}^l}{\Delta x^2}\sin^2{(k_x\Delta x/2)}-4\frac{u_{ij}^l}{\Delta y^2}\sin^2{(k_y\Delta y/2)}\right)=-4\frac{u_{ij}^l}{\Delta t^2}\sin^2{(\omega\Delta t/2)}, $$

 
resulting in

 
$$ \sin{(\omega\Delta t/2)}=\pm \sqrt{\lambda}\Delta t\left(\frac{1}{\Delta x^2}\sin^2{(k_x\Delta x/2)}+\frac{1}{\Delta y^2}\sin^2{(k_y\Delta y/2)}\right)^{1/2}. $$

The squared sine functions can at most be unity. The frequency \( \omega \) is real and our wave is neither damped nor amplified.

Two-dimensional wave equation

 
$$ \sin{(\omega\Delta t/2)}=\pm \sqrt{\lambda}\Delta t\left(\frac{1}{\Delta x^2}\sin^2{(k_x\Delta x/2)}+\frac{1}{\Delta y^2}\sin^2{(k_y\Delta y/2)}\right)^{1/2}. $$

 
The squared sine functions can at most be unity. \( \omega \) is real and our wave is neither damped nor amplified. The numerical \( \omega \) must also be real which is the case when \( \sin{(\omega\Delta t/2)} \) is less than or equal to unity, meaning that

 
$$ \Delta t \le \frac{1}{\sqrt{\lambda}}\left(\frac{1}{\Delta x^2}+\frac{1}{\Delta y^2}\right)^{-1/2}. $$

Two-dimensional wave equation

 
$$ \frac{\partial^2 u}{\partial t^2} = \nabla\cdot (\lambda(x,y) \nabla u). $$

 
If \( \lambda \) is constant, we obtain the standard wave equation discussed in the two previous points. The solution \( u(x,y,t) \) could represent a model for water waves. It represents then the surface elevation from still water. We can model \( \lambda \) as

 
$$ \lambda = gH(x,y), $$

 
with \( g \) being the acceleration of gravity and \( H(x,y) \) is the still water depth.

The function \( H(x,y) \) simulates the water depth using for example measurements of still water depths in say a fjord or the north sea. The boundary conditions are then determined by the coast lines as discussed in point d) below. We have assumed that the vertical motion is negligible and that we deal with long wavelenghts \( \tilde{\lambda} \) compared with the depth of the sea \( H \), that is \( \tilde{\lambda}/H \gg 1 \). We neglect normally Coriolis effects in such calculations.

Two-dimensional wave equation

 
$$ \nabla \cdot (\lambda(x,y) \nabla u)= \frac{\partial }{\partial x}\left(\lambda(x,y)\frac{\partial u}{\partial x}\right)+ \frac{\partial }{\partial y}\left(\lambda(x,y)\frac{\partial u}{\partial y}\right), $$

 
as follows using again a quadratic domain for \( x \) and \( y \):

 
$$ \frac{\partial }{\partial x}\left(\lambda(x,y)\frac{\partial u}{\partial x}\right)\approx \frac{1}{\Delta x} \left(\lambda_{i+1/2,j}\left[\frac{u_{i+1,j}^l-u_{i,j}^l}{\Delta x}\right] -\lambda_{i-1/2,j}\left[\frac{u_{i,j}^l-u_{i-1,j}^l}{\Delta x}\right]\right), $$

 
and

 
$$ \frac{\partial }{\partial y}\left(\lambda(x,y)\frac{\partial u}{\partial y}\right)\approx \frac{1}{\Delta y} \left(\lambda_{i,j+1/2}\left[\frac{u_{i,j+1}^l-u_{i,j}^l}{\Delta y}\right] -\lambda_{i,j-1/2}\left[\frac{u_{i,j}^l-u_{i,j-1}^l}{\Delta y}\right]\right). $$

Two-dimensional wave equation

 
$$ \frac{d}{dx}\left(\lambda(x)\frac{du}{dx}\right)|_{x=x_i} \approx \frac{1}{\Delta x}\left(\lambda\frac{du}{dx}|_{x=x_{i+1/2}}-\lambda\frac{du}{dx}|_{x=x_{i-1/2}}\right), $$

 
where we approximated it at the midpoint by going half a step to the right and half a step to the left. Then we approximate

 
$$ \lambda\frac{du}{dx}|_{x=x_{i+1/2}}\approx \lambda_{x_{i+1/2}}\frac{u_{i+1}-u_i}{\Delta x}, $$

 
and similarly for \( x = x_i-1/2 \).

Overview of week 43

Overview of week 43

Numerical integration and Equal Step Methods

Equal Step Methods

Equal Step Methods

Trapezoidal Rule

double trapezoidal_rule(double a, double b, int n,
                        double (*func)(double))
{
      double trapez_sum;
      double fa, fb, x, step;
      int    j;
      step=(b-a)/((double) n);
      fa=(*func)(a)/2. ;
      fb=(*func)(b)/2. ;
      trapez_sum=0.;
      for (j=1; j <= n-1; j++){
         x=j*step+a;
         trapez_sum+=(*func)(x);

      trapez_sum=(trapez_sum+fb+fa)*step;
      return trapez_sum;
}  // end function for trapezoidal rule

Trapezoidal Rule

double trapezoidal_rule(double a, double b, int n,
                        double (*func)(double))

We call this function simply as something like this

integral =  trapezoidal_rule(a, b, n,  mysuperduperfunction);

Equal Step Methods

Equal Step Methods

Equal Step Methods

 
$$ I=\int_a^bf(x)dx \approx \sum_{i=1}^N\omega_if(x_i), $$

 
where \( \omega \) and \( x \) are the weights and the chosen mesh points, respectively. Simpson's rule gives

 
$$ \omega : \left\{h/3,4h/3,2h/3,4h/3,\dots,4h/3,h/3\right\}, $$

 
for the weights, while the trapezoidal rule resulted in

 
$$ \omega : \left\{h/2,h,h,\dots,h,h/2\right\}. $$

 
In general, an integration formula which is based on a Taylor series using \( N \) points, will integrate exactly a polynomial \( P \) of degree \( N-1 \). That is, the \( N \) weights \( \omega_n \) can be chosen to satisfy \( N \) linear equations

Equal Step Methods, Polynomials and Newton-Cotes

 
$$ p_n(x_j) = y_j\quad j=0,\dots,n $$

 
In the Lagrange representation this interpolation polynomial is given by

 
$$ p_n = \sum_{k=0}^nl_ky_k, $$

 
with the Lagrange factors

 
$$ l_k(x) = \prod_{\begin{array}{c}i=0 \\ i\ne k\end{array}}^n\frac{x-x_i}{x_k-x_i}\quad k=0,\dots,n $$

 
Example: \( n=1 \)

 
$$ p_1(x) = y_0\frac{x-x_1}{x_0-x_1}+y_1\frac{x-x_0}{x_1-x_0}=\frac{y_1-y_0}{x_1-x_0}x-\frac{y_1x_0+y_0x_1}{x_1-x_0}, $$

 
which we recognize as the equation for a straight line.

Equal Step Methods, Polynomials and Newton-Cotes

 
$$ \int_a^bf(x)dx \approx \int_a^bp_n(x)dx = \sum_{k=0}^nw_kf(x_k) $$

 
with

 
$$ w_k = h\frac{(-1)^{n-k}}{k!(n-k)!}\int_0^n\prod_{\begin{array}{c}j=0 \\ j\ne k\end{array}}^n(z-j)dz, $$

 
for \( k=0,\dots,n \).

Equal Step Methods, Polynomials and Newton-Cotes

 
$$ \int_a^bf(x)dx -\frac{b-a}{2}\left[f(a)+f(b)\right]=-\frac{h^3}{12}f^{(2)}(\xi), $$

 
and the global error (composite formula)

 
$$ \int_a^bf(x)dx -T_h(f)=-\frac{b-a}{12}h^2f^{(2)}(\xi). $$

 
For Simpson's rule we have

 
$$ \int_a^bf(x)dx -\frac{b-a}{6}\left[f(a)+4f((a+b)/2)+f(b)\right]=-\frac{h^5}{90}f^{(4)}(\xi), $$

 
and the global error

 
$$ \int_a^bf(x)dx -S_h(f)=-\frac{b-a}{180}h^4f^{(4)}(\xi). $$

 
with \( \xi\in[a,b] \).

Gaussian Quadrature

 
$$ f(x)\approx P_{n}(x), $$

 
with \( n+1 \) mesh points we should be able to integrate exactly the polynomial \( P_{n} \).

Gaussian quadrature methods promise more than this. We can get a better polynomial approximation with order greater than \( n+1 \) to \( f(x) \) and still get away with only \( n+1 \) mesh points. More precisely, we approximate

 
$$ f(x) \approx P_{2n+1}(x), $$

 
and with only \( n+1 \) mesh points these methods promise that

 
$$ \int f(x)dx \approx \int P_{2n+1}(x)dx=\sum_{i=0}^{n} P_{2n+1}(x_i)\omega_i, $$

Gaussian Quadrature

A greater precision for a given amount of numerical work can be achieved if we are willing to give up the requirement of equally spaced integration points. In Gaussian quadrature (hereafter GQ), both the mesh points and the weights are to be determined. The points will not be equally spaced The theory behind GQ is to obtain an arbitrary weight \( \omega \) through the use of so-called orthogonal polynomials. These polynomials are orthogonal in some interval say e.g., [-1,1]. Our points \( x_i \) are chosen in some optimal sense subject only to the constraint that they should lie in this interval. Together with the weights we have then \( 2(n+1) \) (\( n+1 \) the number of points) parameters at our disposal.

Gaussian Quadrature

 
$$ I=\int_a^bf(x)dx =\int_a^bW(x)g(x)dx\approx \sum_{i=0}^n\omega_ig(x_i), $$

 
where \( g \) is smooth and \( W \) is the weight function, which is to be associated with a given orthogonal polynomial.

Gaussian Quadrature

Legendre

 
$$ I=\int_{-1}^{1}f(x)dx $$

 
$$ C(1-x^2)P-m_l^2P+(1-x^2)\frac{d}{dx}\left((1-x^2)\frac{dP}{dx}\right)=0. $$

 
\( C \) is a constant. For \( m_l=0 \) we obtain the Legendre polynomials as solutions, whereas \( m_l \ne 0 \) yields the so-called associated Legendre polynomials. The corresponding polynomials \( P \) are

 
$$ L_k(x)=\frac{1}{2^kk!}\frac{d^k}{dx^k}(x^2-1)^k \quad k=0,1,2,\dots, $$

 
which, up to a factor, are the Legendre polynomials \( L_k \). The latter fulfil the orthorgonality relation

 
$$ \int_{-1}^1L_i(x)L_j(x)dx=\frac{2}{2i+1}\delta_{ij}, $$

 
and the recursion relation

 
$$ (j+1)L_{j+1}(x)+jL_{j-1}(x)-(2j+1)xL_j(x)=0. $$

Laguerre

 
$$ I=\int_0^{\infty}f(x)dx =\int_0^{\infty}x^{\alpha}e^{-x}g(x)dx. $$

 
These polynomials arise from the solution of the differential equation

 
$$ \left(\frac{d^2 }{dx^2}-\frac{d }{dx}+\frac{\lambda}{x}-\frac{l(l+1)}{x^2}\right){\cal L}(x)=0, $$

 
where \( l \) is an integer \( l\ge 0 \) and \( \lambda \) a constant. They fulfil the orthorgonality relation

 
$$ \int_{-\infty}^{\infty}e^{-x}{\cal L}_n(x)^2dx=1, $$

 
and the recursion relation

 
$$ (n+1){\cal L}_{n+1}(x)=(2n+1-x){\cal L}_{n}(x)-n{\cal L}_{n-1}(x). $$

Hermite

 
$$ I=\int_{-\infty}^{\infty}f(x)dx =\int_{-\infty}^{\infty}e^{-x^2}g(x)dx. $$

 
we could use the Hermite polynomials in order to extract weights and mesh points. The Hermite polynomials are the solutions of the following differential equation

 
$$ \frac{d^2H(x)}{dx^2}-2x\frac{dH(x)}{dx}+ (\lambda-1)H(x)=0. \tag{26} $$

 
They fulfil the orthorgonality relation

 
$$ \int_{-\infty}^{\infty}e^{-x^2}H_n(x)^2dx=2^nn!\sqrt{\pi}, $$

 
and the recursion relation

 
$$ H_{n+1}(x)=2xH_{n}(x)-2nH_{n-1}(x). $$

Gaussian Quadrature, general Properties

 
$$ \int_a^bW(x)f(x)dx \approx \sum_{i=0}^n\omega_if(x_i), $$

 
with \( n+1 \) distinct quadrature points (mesh points) is a called a Gaussian quadrature formula if it integrates all polynomials \( p\in P_{2n+1} \) exactly, that is

 
$$ \int_a^bW(x)p(x)dx =\sum_{i=0}^n\omega_ip(x_i), $$

 
It is assumed that \( W(x) \) is continuous and positive and that the integral

 
$$ \int_a^bW(x)dx , $$

 
exists. Note that the replacement of \( f\rightarrow Wg \) is normally a better approximation due to the fact that we may isolate possible singularities of \( W \) and its derivatives at the endpoints of the interval.

Numerical integration: A simple example

 
$$ I= \int_0^{\infty} x\exp{(-x)}\sin{x}=\frac{1}{2}, $$

 
using brute force Trapezoidal rule, Simpson's rule, Gauss-Legendre, Gauss-Laguerre and Gauss-Legendre again but with a smarter mapping.

• Before we start it can be useful to study the integrand.
• How should we pick the integration limits?

Simple integral

 
$$ \int_0^{\infty}f(x)dx \approx \int_0^{\Lambda}f(x)dx $$

     int n;
     double a, b, alf, xx;
     cout << "Read in the number of integration points" << endl;
     cin >> n;
     cout << "Read in integration limits" << endl;
     cin >> a >> b;

Simple integral

//   reserve space in memory for vectors containing the mesh points
//   weights and function values for the use of the gauss-legendre
//   method
     double *x = new double [n];
     double *w = new double [n];
     // Gauss-Laguerre is old-fashioned translation of F77 --> C++
     // arrays start at 1 and end at n
     double *xgl = new double [n+1];
     double *wgl = new double [n+1];

Simple integral

//   set up the mesh points and weights
     gauleg(a, b,x,w, n);
//   set up the mesh points and weights
     alf = 1.0;  //  <---  Note alf
     gauss_laguerre(xgl,wgl, n, alf);
//   evaluate the integral with the Gauss-Legendre method
//   Note that we initialize the sum. Here brute force gauleg
     double int_gauss = 0.;
     for ( int i = 0;  i < n; i++){
        int_gauss+=w[i]*int_function(x[i]);

Simple integral, importance integration/sampling

//   evaluate the integral with the Gauss-Laguerre method
//   Note that we initialize the sum
     double int_gausslag = 0.;
     for ( int i = 1;  i <= n; i++){
       int_gausslag+=wgl[i]*sin(xgl[i]);     }

Simple integral

//   evaluate the integral with the Gauss-Laguerre method
//   Here we change the mesh points with a mapping
//   Need to call gauleg from -1 to + 1
     gauleg(-1.0, 1.0,x,w, n);
     double pi_4 = acos(-1.0)*0.25;
     for ( int i = 0;  i < n; i++){
       xx=pi_4*(x[i]+1.0);
       r[i]= tan(xx);
       s[i]=pi_4/(cos(xx)*cos(xx))*w[i];

     double int_gausslegimproved = 0.;
     for ( int i = 0;  i < n; i++){
       int_gausslegimproved += s[i]*int_function(r[i]);

A six-dimensional integral, project 3 2010

 
$$ \Psi({\bf r}_1,{\bf r}_2) = e^{-\alpha (r_1+r_2)}. $$

 
Note that it is not possible to find a closed form solution to Schr\"odinger's equation for two interacting electrons in the helium atom.

The integral we need to solve is the quantum mechanical expectation value of the correlation energy between two electrons, namely

 
$$ \begin{equation}label{eq:correlationenergy} \langle \frac{1}{|{\bf r}_1-{\bf r}_2|} \rangle = \int_{-\infty}^{\infty} d{\bf r}_1d{\bf r}_2 e^{-2\alpha (r_1+r_2)}\frac{1}{|{\bf r}_1-{\bf r}_2|}=\frac{5\pi^2}{16^2}=0.192765711. \end{equation} $$

 
Note that our wave function is not normalized. There is a normalization factor missing.

Brute force, six-dimensional integral, Gauss-Legendre

     double *x = new double [N];
     double *w = new double [N];
//   set up the mesh points and weights
     gauleg(a,b,x,w, N);

//   evaluate the integral with the Gauss-Legendre method
//   Note that we initialize the sum
     double int_gauss = 0.;
     for (int i=0;i<N;i++){
	     for (int j = 0;j<N;j++){
	     for (int k = 0;k<N;k++){
	     for (int l = 0;l<N;l++){
	     for (int m = 0;m<N;m++){
	     for (int n = 0;n<N;n++){
        int_gauss+=w[i]*w[j]*w[k]*w[l]*w[m]*w[n]
       *int_function(x[i],x[j],x[k],x[l],x[m],x[n]);
     		}}}}}
	}

The six-dimensional integral with Gauss-Legendre

//  this function defines the function to integrate
double int_function(double x1, double y1, double z1, double x2, double y2, double z2)
{
   double alpha = 2.;
// evaluate the different terms of the exponential
   double exp1=-2*alpha*sqrt(x1*x1+y1*y1+z1*z1);
   double exp2=-2*alpha*sqrt(x2*x2+y2*y2+z2*z2);
   double deno=sqrt(pow((x1-x2),2)+pow((y1-y2),2)+pow((z1-z2),2));
  if(deno <pow(10.,-6.)) { return 0;}
  else return exp(exp1+exp2)/deno;
} // end of function to evaluate

Switch to spherical coordinates

 
$$ d{\bf r}_1d{\bf r}_2 = r_1^2dr_1 r_2^2dr_2 dcos(\theta_1)dcos(\theta_2)d\phi_1d\phi_2, $$

 
and

 
$$ \frac{1}{r_{12}}= \frac{1}{\sqrt{r_1^2+r_2^2-2r_1r_2cos(\beta)}} $$

 
with

 
$$ \cos(\beta) = \cos(\theta_1)\cos(\theta_2)+\sin(\theta_1)\sin(\theta_2)\cos(\phi_1-\phi_2)) $$

Switch to spherical coordinates

 
$$ \int_0^{\infty} r_1^2dr_1 \int_0^{\infty}r_2^2dr_2 \int_0^{\pi}dcos(\theta_1)\int_0^{\pi}dcos(\theta_2)\int_0^{2\pi}d\phi_1\int_0^{2\pi}d\phi_2 \times $$

 
$$ \frac{e^{-2\alpha (r_1+r_2)}}{\sqrt{r_1^2+r_2^2-2r_1r_2\cos(\theta_1)\cos(\theta_2)+\sin(\theta_1)\sin(\theta_2)\cos(\phi_1-\phi_2)) }} $$

 
since

 
$$ \frac{1}{r_{12}}= \frac{1}{\sqrt{r_1^2+r_2^2-2r_1r_2cos(\beta)}} $$

 
with

 
$$ \cos(\beta) = \cos(\theta_1)\cos(\theta_2)+\sin(\theta_1)\sin(\theta_2)\cos(\phi_1-\phi_2)) $$

Adaptive methods

The above methods are all based on a defined step length, normally provided by the user, dividing the integration domain with a fixed number of subintervals. This is rather simple to implement may be inefficient, in particular if the integrand varies considerably in certain areas of the integration domain. In these areas the number of fixed integration points may not be adequate. In other regions, the integrand may vary slowly and fewer integration points may be needed.

Adaptive methods

Adaptive methods

Step 1. We compute our first approximation by computing the integral for the full domain. We label this as \( I^{(0)} \). It is obtained by calling our previously discussed function trapezoidal_rule as

I0 = trapezoidal_rule(a, b, n, function);

Step 2. We split the integration in two, with \( c= (a+b)/2 \). We compute then the two integrals \( I^{(1L)} \) and \( I^{(1R)} \)

I1L = trapezoidal_rule(a, c, n, function);

and

I1R = trapezoidal_rule(c, b, n, function);

With a given defined tolerance, being a small number provided by us, we estimate the difference \( |I^{(1L)}+I^{(1R)}-I^{(0)}| < \mbox{tolerance} \). If this test is satisfied, our first approximation is satisfactory. If not, we can set up a recursive procedure where the integral is split into subsequent subintervals until our tolerance is satisfied.

Adaptive methods

//     Simple recursive function that implements the
//     adaptive integration using the trapezoidal rule
const int maxrecursions = 50;
const double tolerance = 1.0E-10;
//  Takes as input the integration  limits, number of points, function to integrate
//  and the number of steps
void adaptive_integration(double a, double b, double *Integral, int n, int steps, double (*func)(double))
     if ( steps > maxrecursions){
        cout << 'Too many recursive steps, the function varies too much' << endl;
        break;

Adaptive methods

     double c = (a+b)*0.5;
     // the whole integral
     double I0 = trapezoidal_rule(a, b,n, func);
     //  the left half
     double I1L = trapezoidal_rule(a, c,n, func);
     //  the right half
     double I1R = trapezoidal_rule(c, b,n, func);
     if (fabs(I1L+I1R-I0) < tolerance )  integral = I0;
     else
     {
        adaptive_integration(a, c, integral, int n, ++steps, func)
        adaptive_integration(c, b, integral, int n, ++steps, func)

// end function Adaptive integration

The variables {\bf Integral} and {\bf steps} should be initialized to zero by the function that calls the adaptive procedure.

Overview of week 44

Overview of week 44

Going Parallel: Divide et impera - Divide and Conquer (accordingly after Julius Caesar)

{\bf Task parallelism}: the work of a global problem can be divided into a number of independent tasks, which rarely need to synchronize. Monte Carlo simulation or integrations are examples of this. It is almost embarrassingly trivial to parallelize Monte Carlo and numerical integration codes.

We will use MPI=Message Passing Interface. MPI is a message-passing library where all the routines have corresponding C/C++-binding

   MPI_Command_name

and Fortran-binding (routine names are in uppercase, but can also be in lower case)

   MPI_COMMAND_NAME

What is Message Passing Interface (MPI)? Yet another library!

MPI is a specification, not a particular implementation. MPI programs should be able to run on all possible machines and run all MPI implementetations without change.

An MPI computation is a collection of processes communicating with messages.

See chapter 5.5 of lecture notes for more details.

MPI

• independent of hardware;
• not a language or compiler specification;
• not a specific implementation or product. A message passing standard for portability and ease-of-use. Designed for high performance. Insert communication and synchronization functions where necessary.

Demands from the HPC community

Development of the capacity for single-processor computers can hardly keep up with the pace of scientific computing:

• processor speed
• memory size/speed Solution: parallel computing!

The basic ideas of parallel computing

• Pursuit of shorter computation time and larger simulation size gives rise to parallel computing.
• Multiple processors are involved to solve a global problem.
• The essence is to divide the entire computation evenly among collaborative processors. Divide and conquer.

A rough classification of hardware models

• SIMD (single-instruction-multiple-data) machines incorporate the idea of parallel processing, which use a large number of process- ing units to execute the same instruction on different data.
• Modern parallel computers are so-called MIMD (multiple-instruction- multiple-data) machines and can execute different instruction streams in parallel on different data.

Shared memory and distributed memory

• One way of categorizing modern parallel computers is to look at the memory configuration.
• In shared memory systems the CPUs share the same address space. Any CPU can access any data in the global memory.
• In distributed memory systems each CPU has its own memory. The CPUs are connected by some network and may exchange messages.

Different parallel programming paradigms

• {\bf Task parallelism:} the work of a global problem can be divided into a number of independent tasks, which rarely need to synchronize. Monte Carlo simulation is one example. Integration is another. However this paradigm is of limited use.
• {\bf Data parallelism:} use of multiple threads (e.g. one thread per processor) to dissect loops over arrays etc. This paradigm requires a single memory address space. Communication and synchronization between processors are often hidden, thus easy to program. However, the user surrenders much control to a specialized compiler. Examples of data parallelism are compiler-based parallelization and OpenMP directives.

Today's situation of parallel computing

• Distributed memory is the dominant hardware configuration. There is a large diversity in these machines, from MPP (massively parallel processing) systems to clusters of off-the-shelf PCs, which are very cost-effective.
• Message-passing is a mature programming paradigm and widely accepted. It often provides an efficient match to the hardware. It is primarily used for the distributed memory systems, but can also be used on shared memory systems. In these lectures we consider only message-passing for writing parallel programs.

Overhead present in parallel computing

• {\bf Uneven load balance}: not all the processors can perform useful work at all time.
• {\bf Overhead of synchronization.}
• {\bf Overhead of communication}.
• {Extra computation due to parallelization}. Due to the above overhead and that certain part of a sequential algorithm cannot be parallelized we may not achieve an optimal parallelization.

Parallelizing a sequential algorithm

• Identify the part(s) of a sequential algorithm that can be executed in parallel. This is the difficult part,
• Distribute the global work and data among \( P \) processors.

Process and processor

• We refer to process as a logical unit which executes its own code, in an MIMD style.
• The processor is a physical device on which one or several processes are executed.
• The MPI standard uses the concept process consistently throughout its documentation.

Bindings to MPI routines

   MPI_Command_name

and Fortran-binding (routine names are in uppercase, but can also be in lower case)

   MPI_COMMAND_NAME

The discussion in these slides focuses on the C++ binding.

The most important/used MPI routines

• MPI_ Init - initiate an MPI computation
• MPI_Finalize - terminate the MPI computation and clean up
• MPI_Comm_size - how many processes participate in a given MPI communicator?
• MPI_Comm_rank - which one am I? (A number between 0 and size-1.)
• MPI_Reduce(Allreduce) - Collect data from all nodes and either sum them up in one or all (Allreduce). Useful for numerical integration
• MPI_Send - send a message to a particular process within an MPI communicator
• MPI_Recv - receive a message from a particular process within an MPI communicator

Note the MPI_COMM_WORLD declaration

• A group of MPI processes with a name (context).
• Any process is identified by its rank. The rank is only meaningful within a particular communicator.
• By default communicator MPI_COMM_WORLD contains all the MPI processes.
• Mechanism to identify subset of processes.
• Promotes modular design of parallel libraries.

The first MPI C/C++ program

using namespace std;
#include <mpi.h>
#include <iostream>
int main (int nargs, char* args[])
{
int numprocs, my_rank;
//   MPI initializations
MPI_Init (&nargs, &args);
MPI_Comm_size (MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank (MPI_COMM_WORLD, &my_rank);
cout << "Hello world, I have  rank " << my_rank << " out of "
     << numprocs << endl;
//  End MPI
MPI_Finalize ();

The Fortran program

PROGRAM hello
INCLUDE "mpif.h"
INTEGER:: size, my_rank, ierr

CALL  MPI_INIT(ierr)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierr)
CALL MPI_COMM_RANK(MPI_COMM_WORLD, my_rank, ierr)
WRITE(*,*)"Hello world, I've rank ",my_rank," out of ",size
CALL MPI_FINALIZE(ierr)

END PROGRAM hello

Note 1

Ordered output with MPI_Barrier

int main (int nargs, char* args[])
{
 int numprocs, my_rank, i;
 MPI_Init (&nargs, &args);
 MPI_Comm_size (MPI_COMM_WORLD, &numprocs);
 MPI_Comm_rank (MPI_COMM_WORLD, &my_rank);
 for (i = 0; i < numprocs; i++) {}
 MPI_Barrier (MPI_COMM_WORLD);
 if (i == my_rank) {
 cout << "Hello world, I have  rank " << my_rank <<
        " out of " << numprocs << endl;}
      MPI_Finalize ();

Note 2

Strategies

• Develop codes locally, run with some few processes and test your codes. Do benchmarking, timing and so forth on local nodes, for example your laptop. You can install MPICH2 on your laptop (most new laptos come with dual cores). You can test with one node at the lab.
• When you are convinced that your codes run correctly, you start your production runs on available supercomputers. At UiO we have a new machine among the top 100, Abel.

How do I run MPI on the machines at the lab (MPICH2)

Step 1. Compile with mpicxx or mpic++

Step 2. Set up collaboration between processes and run

mpd --ncpus=4 & # run code with mpiexec -n 4 ./nameofprog

Here we declare that we will use 4 processes via the -ncpus option and via -n 4 when running.

Step 3. End with

mpdallexit

Can I do it on my own PC/laptop?

• go to http://www.mcs.anl.gov/research/projects/mpich2/
• follow the instructions and install it on your own PC/laptop It works on Windows, Mac and Linux. For Linux (Ubuntu, Linux Mint), go to synaptic package manager and search for MPI. You will get several options.

Ordered output with MPI_Recv and MPI_Send

.....
int numprocs, my_rank, flag;
MPI_Status status;
MPI_Init (&nargs, &args);
MPI_Comm_size (MPI_COMM_WORLD, &numprocs);
MPI_Comm_rank (MPI_COMM_WORLD, &my_rank);
if (my_rank > 0)
MPI_Recv (&flag, 1, MPI_INT, my_rank-1, 100,
           MPI_COMM_WORLD, &status);
cout << "Hello world, I have  rank " << my_rank << " out of "
<< numprocs << endl;
if (my_rank < numprocs-1)
MPI_Send (&my_rank, 1, MPI_INT, my_rank+1,
          100, MPI_COMM_WORLD);
MPI_Finalize ();

Note 3

int MPI_Send(void *buf, int count,
             MPI_Datatype datatype,
             int dest, int tag, MPI_Comm comm)}

This single command allows the passing of any kind of variable, even a large array, to any group of tasks. The variable {\bf buf} is the variable we wish to send while {\bf count} is the number of variables we are passing. If we are passing only a single value, this should be 1. If we transfer an array, it is the overall size of the array. For example, if we want to send a 10 by 10 array, count would be \( 10\times 10=100 \) since we are actually passing 100 values.

Note 4

int MPI_Recv( void *buf, int count, MPI_Datatype datatype,
            int source,
            int tag, MPI_Comm comm, MPI_Status *status )

The arguments that are different from those in MPI_SEND are {\bf buf} which is the name of the variable where you will be storing the received data, {\bf source} which replaces the destination in the send command. This is the return ID of the sender.

Finally, we have used MPI_Status status where one can check if the receive was completed.

The output of this code is the same as the previous example, but now process 0 sends a message to process 1, which forwards it further to process 2, and so forth.

Armed with this wisdom, performed all hello world greetings, we are now ready for serious work.

Integrating \( \pi \)

Integration algos

 
$$ I=\int_a^bf(x) dx=h\left(f(a)/2 + f(a+h) +f(a+2h)+ \dots +f(b-h)+ f_{b}/2\right). $$

Another very simple approach is the so-called midpoint or rectangle method. In this case the integration area is split in a given number of rectangles with length \( h \) and heigth given by the mid-point value of the function. This gives the following simple rule for approximating an integral

 
$$ I=\int_a^bf(x) dx \approx h\sum_{i=1}^N f(x_{i-1/2}), $$

 
where \( f(x_{i-1/2}) \) is the midpoint value of \( f \) for a given rectangle. This is used in program5.cpp.

MPI_reduce and MPI_Allreduce

MPI_reduce

MPI_reduce( void *senddata, void* resultdata, int count,
     MPI_Datatype datatype, MPI_Op, int root, MPI_Comm comm)

The two variables \( senddata \) and \( resultdata \) are obvious, besides the fact that one sends the address of the variable or the first element of an array. If they are arrays they need to have the same size. The variable \( count \) represents the total dimensionality, 1 in case of just one variable, while MPI_Datatype defines the type of variable which is sent and received.

The new feature is MPI_Op. It defines the type of operation we want to do. In our case, since we are summing the rectangle contributions from every process we define MPI_Op = MPI_SUM. If we have an array or matrix we can search for the largest og smallest element by sending either MPI_MAX or MPI_MIN. If we want the location as well (which array element) we simply transfer MPI_MAXLOC or MPI_MINOC. If we want the product we write MPI_PROD.

MPI_Allreduce is defined as

MPI_Alreduce( void *senddata, void* resultdata, int count,
          MPI_Datatype datatype, MPI_Op, MPI_Comm comm)}.

Dissection of example program6.cpp

//    Trapezoidal rule and numerical integration usign MPI, example program6.cpp
using namespace std;
#include <mpi.h>
#include <iostream>

//     Here we define various functions called by the main program

double int_function(double );
double trapezoidal_rule(double , double , int , double (*)(double));

//   Main function begins here
int main (int nargs, char* args[])
{
  int n, local_n, numprocs, my_rank;
  double a, b, h, local_a, local_b, total_sum, local_sum;
  double  time_start, time_end, total_time;

Dissection of example program6.cpp

  //  MPI initializations
  MPI_Init (&nargs, &args);
  MPI_Comm_size (MPI_COMM_WORLD, &numprocs);
  MPI_Comm_rank (MPI_COMM_WORLD, &my_rank);
  time_start = MPI_Wtime();
  //  Fixed values for a, b and n
  a = 0.0 ; b = 1.0;  n = 1000;
  h = (b-a)/n;    // h is the same for all processes
  local_n = n/numprocs;
  // make sure n > numprocs, else integer division gives zero
  // Length of each process' interval of
  // integration = local_n*h.
  local_a = a + my_rank*local_n*h;
  local_b = local_a + local_n*h;

Dissection of example program6.cpp

  total_sum = 0.0;
  local_sum = trapezoidal_rule(local_a, local_b, local_n,
                               &int_function);
  MPI_Reduce(&local_sum, &total_sum, 1, MPI_DOUBLE,
              MPI_SUM, 0, MPI_COMM_WORLD);
  time_end = MPI_Wtime();
  total_time = time_end-time_start;
  if ( my_rank == 0) {
    cout << "Trapezoidal rule = " <<  total_sum << endl;
    cout << "Time = " <<  total_time
         << " on number of processors: "  << numprocs  << endl;

  // End MPI
  MPI_Finalize ();
  return 0;
}  // end of main program

Dissection of example program6.cpp

//  this function defines the function to integrate
double int_function(double x)
{
  double value = 4./(1.+x*x);
  return value;
} // end of function to evaluate

Dissection of example program6.cpp

//  this function defines the trapezoidal rule
double trapezoidal_rule(double a, double b, int n,
                         double (*func)(double))
{
  double trapez_sum;
  double fa, fb, x, step;
  int    j;
  step=(b-a)/((double) n);
  fa=(*func)(a)/2. ;
  fb=(*func)(b)/2. ;
  trapez_sum=0.;
  for (j=1; j <= n-1; j++){
    x=j*step+a;
    trapez_sum+=(*func)(x);

  trapez_sum=(trapez_sum+fb+fa)*step;
  return trapez_sum;
}  // end trapezoidal_rule

Plan for Monte Carlo Lectures, chapters 11-14 in Lecture notes

• This week: intro, MC integration and probability distribution functions (PDFs)
• Next week: More on integration, PDFs, MC integration and random walks.
• Third week: random walks and statistical physics.
• Fourth week: Statistical physics.
• Fifth week: Most likely quantum Monte Carlo Approximately from this week till the end of November.

Monte Carlo Keywords

• Be able to generate random variables following a given probability distribution function PDF
• Find a probability distribution function (PDF).
• Sampling rule for accepting a move
• Compute standard deviation and other expectation values
• Techniques for improving errors Enhances algorithmic thinking!

Probability Distribution Functions PDF

Expectation Values

 
$$ E[x^k]= \langle x^k\rangle=\frac{1}{N}\sum_{i=1}^{N}x_i^kp(x_i), $$

 
provided that the sums (or integrals) $ \sum_{i=1}^{N}p(x_i)$ converge absolutely (viz , $ \sum_{i=1}^{N}|p(x_i)|$ converges)

Continuous PDF

 
$$ E[x^k]=\langle x^k\rangle=\int_a^b x^kp(x)dx, $$

Function \( f(x) \)

 
$$ E[f^k]=\langle f^k\rangle=\int_a^b f^kp(x)dx, $$

Variance

 
$$ \sigma^2_f=E[f^2]-(E[f])^2=\langle f^2\rangle-\langle f\rangle^2 $$

Important PDFs

 
$$ p(x)=\frac{1}{b-a}\Theta(x-a)\Theta(b-x), $$

 
which gives for \( a=0,b=1 \) $p(x)=1$ for \( x\in [0,1] \) and zero else.

exponential distribution

 
$$ p(x)=\alpha e^{-\alpha x}, $$

 
with probability different from zero in \( [0,\infty] \)

normal distribution (Gaussian)

 
$$ p(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp{\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)} $$

 
with probability different from zero in \( [-\infty,\infty] \)

All random number generators use the uniform distribution for \( x\in[0,1] \).

Why Monte Carlo?

 
$$ \langle H \rangle = $$

 
$$ \frac{\int d{\bf R}_1d{\bf R}_2\dots d{\bf R}_N \Psi^{\ast}({\bf R_1},{\bf R}_2,\dots,{\bf R}_N) H({\bf R_1},{\bf R}_2,\dots,{\bf R}_N) \Psi({\bf R_1},{\bf R}_2,\dots,{\bf R}_N)} {\int d{\bf R}_1d{\bf R}_2\dots d{\bf R}_N \Psi^{\ast}({\bf R_1},{\bf R}_2,\dots,{\bf R}_N) \Psi({\bf R_1},{\bf R}_2,\dots,{\bf R}_N)}, $$

an in general intractable problem.

This integral is actually the starting point in a Variational Monte Carlo calculation. {\bf Gaussian quadrature: Forget it!} given 10 particles and 10 mesh points for each degree of freedom and an ideal 1 petaflops machine (all operations take the same time), how long will it ta ke to compute the above integral? Lifetime of the universe $T\approx 4.7 \times 10^{17}$s.

More on dimensionality

 
$$ \hat{H}\Psi({\bf r}_1,..,{\bf r}_A,\alpha_1,..,\alpha_A)=E\Psi({\bf r}_1,..,{ \bf r}_A,\alpha_1,..,\alpha_A) $$

 
where

 
$$ {\bf r}_1,..,{\bf r}_A, $$

 
are the coordinates and

 
$$ \alpha_1,..,\alpha_A, $$

 
are sets of relevant quantum numbers such as spin and isospin for a system of \( A \) nucleons (\( A=N+Z \), \( N \) being the number of neutrons and \( Z \) the number of protons).

Even more on dimensionality

 
$$ 2^A\times \left(\begin{array}{c} A\\ Z\end{array}\right) $$

 
coupled second-order differential equations in \( 3A \) dimensions.

For a nucleus like $^{10}$Be this number is {\bf 215040}. This is a truely challenging many-body problem.

But what do we gain by Monte Carlo Integration?

 
$$ I=\int_0^1 f(x)dx\approx \frac{1}{N}\sum_{i=1}^Nf(x_i), $$

 
can be rewritten using the concept of the average of the function \( f \) for a given PDF \( p(x) \) as

 
$$ E[f]=\langle f \rangle = \frac{1}{N}\sum_{i=1}^Nf(x_i)p(x_i), $$

 
and identify \( p(x) \) with the uniform distribution, viz $ p(x)=1$ when \( x\in [0,1] \) and zero for all other values of \( x \). The integral is then the average of \( f \) over the interval \( x \in [0,1] \)

 
$$ I=\int_0^1 f(x)dx\approx E[f]=\langle f \rangle. $$

But what do we gain by Monte Carlo Integration?

 
$$ \sigma^2_f=\frac{1}{N}\sum_{i=1}^N(f(x_i)-\langle f\rangle)^2p(x_i), $$

 
and inserting the uniform distribution this yields

 
$$ \sigma^2_f=\frac{1}{N}\sum_{i=1}^Nf(x_i)^2- \left(\frac{1}{N}\sum_{i=1}^Nf(x_i)\right)^2, $$

 
or

 
$$ \sigma^2_f=E[f^2]-(E[f])^2=\left(\langle f^2\rangle - \langle f \rangle^2\right). $$

 
which is nothing but a measure of the extent to which \( f \) deviates from its average over the region of integration. The standard deviation is defined as the square root of the variance.

But what do we gain by Monte Carlo Integration?

 
$$ \langle I \rangle_M=\frac{1}{M}\sum_{l=1}^{M}\langle f\rangle_l. $$

 
If we can consider the probability of correlated events to be zero, we can rewrite the variance of these series of measurements as (equating \( M=N \))

 
$$ \sigma^2_N\approx \frac{1}{N}\left(\langle f^2\rangle - \langle f \rangle^2\right)=\frac{\sigma^2_f}{N}. $$

We note that the standard deviation is proportional with the inverse square root of the number of measurements

 
$$ \sigma_N \sim \frac{1}{\sqrt{N}}. $$

 
The aim in Monte Carlo calculations is to have \( \sigma_N \) as small as possible after \( N \) samples. The results from one sample represents, since we are using concepts from statistics, a 'measurement'.

But what do we gain by Monte Carlo Integration?

• We saw that the trapezoidal rule carries a truncation error \( O(h^2) \), with \( h \) the step length.

• Quadrature rules such as Newton-Cotes have a truncation error which goes like \( \sim O(h^k) \), with \( k \ge 1 \). Recalling that the step size is defined as \( h=(b-a)/N \), we have an error which goes like \( \sim N^{-k} \).
• Monte Carlo integration is more efficient in higher dimensions. Assume that our integration volume is a hypercube with side \( L \) and dimension \( d \). This cube contains hence \( N=(L/h)^d \) points and therefore the error in the result scales as \( N^{-k/d} \) for the traditional methods.
• The error in the Monte carlo integration is however independent of \( d \) and scales as \( \sigma\sim 1/\sqrt{N} \), always!

* Comparing this with traditional methods, shows that Monte Carlo integration is more efficient than an order-k algorithm when \( d > 2k \)

Some simple examples: Example 1: Particles in a Box

After some time the system reaches its equilibrium state with equally many particles in both halves, \( N/2 \). Instead of determining complicated initial conditions for a system of \( N \) particles, we model the system by a simple statistical model. In order to simulate this system, which may consist of \( N \gg 1 \) particles, we assume that all particles in the left half have equal probabilities of going to the right half.

Particles in a Box

• Choose the number of particles \( N \).
• Make a loop over time, where the maximum time should be larger than the number of particles \( N \).
• For every time step \( \Delta t \) there is a probability \( n_l/N \) for a move to the right. Compare this probability with a random number \( x \).
• If $ x \le n_l/N$, decrease the number of particles in the left half by one, i.e., \( n_l=n_l-1 \). Else, move a particle from the right half to the left, i.e., \( n_l=n_l+1 \). {\bf This is our sampling rule}
• Increase the time by one unit (the external loop). In this case, a Monte Carlo sample corresponds to one time unit \( \Delta t \).

Particles in a Box

  // setup of initial conditions
  nleft = initial_n_particles;
  max_time = 10*initial_n_particles;
  idum = -1;
  // sampling over number of particles
  for( time=0; time <= max_time; time++){
    random_n = ((int) initial_n_particles*ran0(&idum));
    if ( random_n <= nleft){
      nleft -= 1;

    else{
      nleft += 1;

    ofile << setiosflags(ios::showpoint | ios::uppercase);
    ofile << setw(15) << time;
    ofile << setw(15) << nleft << endl;

Example 2: Radioactive Decay

 
$$ dN(t)=-\omega N(t)dt, $$

 
whose solution is

 
$$ N(t)=N(0)e^{-\omega t}, $$

 
where we have defined the mean lifetime \( \tau \) of \( X \) as

 
$$ \tau =\frac{1}{\omega}. $$

If a nucleus \( X \) decays to a daugther nucleus \( Y \) which also can decay, we get the following coupled equations

 
$$ \frac{dN_X(t)}{dt}=-\omega_XN_X(t), $$

 
and

 
$$ \frac{dN_Y(t)}{dt}=-\omega_YN_Y(t)+\omega_XN_X(t). $$

Radioactive Decay

 
$$ \frac{\Delta N(t)}{N(t)\Delta t}= -\lambda $$

 
\( \lambda \) is inversely proportional to the lifetime

• Choose the number of particles \( N(t=0)=N_0 \).
• Make a loop over the number of time steps, with maximum time bigger than the number of particles \( N_0 \)
• At every time step there is a probability \( \lambda \) for decay. Compare this probability with a random number \( x \).
• If $ x \le \lambda$, reduce the number of particles with one i.e., \( N=N-1 \). If not, keep the same number of particles till the next time step. {\bf This is our sampling rule}
• Increase by one the time step (the external loop)

Radioactive Decay

  idum=-1;  // initialise random number generator
  // loop over monte carlo cycles
  // One monte carlo loop is one sample
  for (cycles = 1; cycles <= number_cycles; cycles++){
    n_unstable = initial_n_particles;
    //  accumulate the number of particles per time step per trial
    ncumulative[0] += initial_n_particles;
    // loop over each time step
    for (time=1; time <= max_time; time++){
      // for each time step, we check each particle
      particle_limit = n_unstable;
      for ( np = 1; np <=  particle_limit; np++) {
        if( ran0(&idum) <= decay_probability) {
          n_unstable=n_unstable-1;

      }  // end of loop over particles
      ncumulative[time] += n_unstable;
    }  // end of loop over time steps
  }    // end of loop over MC trials
}   // end mc_sampling function

Monte Carlo Keywords

• Be able to generate random variables following a given probability distribution function PDF
• Find a probability distribution function (PDF).
• Sampling rule for accepting a move
• Compute standard deviation and other expectation values
• Techniques for improving errors Enhances algorithmic thinking!

Why Monte Carlo integration?

 
$$ \langle H \rangle = $$

 
$$ \frac{\int d{\bf R}_1d{\bf R}_2\dots d{\bf R}_N \Psi^{\ast}({\bf R_1},{\bf R}_2,\dots,{\bf R}_N) H({\bf R_1},{\bf R}_2,\dots,{\bf R}_N) \Psi({\bf R_1},{\bf R}_2,\dots,{\bf R}_N)} {\int d{\bf R}_1d{\bf R}_2\dots d{\bf R}_N \Psi^{\ast}({\bf R_1},{\bf R}_2,\dots,{\bf R}_N) \Psi({\bf R_1},{\bf R}_2,\dots,{\bf R}_N)}, $$

an in general intractable problem.

This integral is actually the starting point in a Variational Monte Carlo calculation. {\bf Gaussian quadrature: Forget it!} given 10 particles and 10 mesh points for each degree of freedom and an ideal 1 petaflops machine (all operations take the same time), how long will it ta ke to compute the above integral? Lifetime of the universe $T\approx 4.7 \times 10^{17}$s.

But what do we gain by Monte Carlo Integration?

 
$$ I=\int_0^1 f(x)dx\approx \frac{1}{N}\sum_{i=1}^Nf(x_i), $$

 
can be rewritten using the concept of the average of the function \( f \) for a given PDF \( p(x) \) as

 
$$ E[f]=\langle f \rangle = \frac{1}{N}\sum_{i=1}^Nf(x_i)p(x_i), $$

 
and identify \( p(x) \) with the uniform distribution, viz $ p(x)=1$ when \( x\in [0,1] \) and zero for all other values of \( x \). The integral is then the average of \( f \) over the interval \( x \in [0,1] \)

 
$$ I=\int_0^1 f(x)dx\approx E[f]=\langle f \rangle. $$

But what do we gain by Monte Carlo Integration?

 
$$ \sigma^2_f=\frac{1}{N}\sum_{i=1}^N(f(x_i)-\langle f\rangle)^2p(x_i), $$

 
and inserting the uniform distribution this yields

 
$$ \sigma^2_f=\frac{1}{N}\sum_{i=1}^Nf(x_i)^2- \left(\frac{1}{N}\sum_{i=1}^Nf(x_i)\right)^2, $$

 
or

 
$$ \sigma^2_f=E[f^2]-(E[f])^2=\left(\langle f^2\rangle - \langle f \rangle^2\right). $$

 
which is nothing but a measure of the extent to which \( f \) deviates from its average over the region of integration. The standard deviation is defined as the square root of the variance.

But what do we gain by Monte Carlo Integration?

 
$$ \langle I \rangle_M=\frac{1}{M}\sum_{l=1}^{M}\langle f\rangle_l. $$

 
If we can consider the probability of correlated events to be zero, we can rewrite the variance of these series of measurements as (equating \( M=N \))

 
$$ \sigma^2_N\approx \frac{1}{N}\left(\langle f^2\rangle - \langle f \rangle^2\right)=\frac{\sigma^2_f}{N}. $$

We note that the standard deviation is proportional with the inverse square root of the number of measurements

 
$$ \sigma_N \sim \frac{1}{\sqrt{N}}. $$

 
The aim in Monte Carlo calculations is to have \( \sigma_N \) as small as possible after \( N \) samples. The results from one sample represents, since we are using concepts from statistics, a 'measurement'.

But what do we gain by Monte Carlo Integration?

• We saw that the trapezoidal rule carries a truncation error \( O(h^2) \), with \( h \) the step length.

• Quadrature rules such as Newton-Cotes have a truncation error which goes like \( \sim O(h^k) \), with \( k \ge 1 \). Recalling that the step size is defined as \( h=(b-a)/N \), we have an error which goes like \( \sim N^{-k} \).
• Monte Carlo integration is more efficient in higher dimensions. Assume that our integration volume is a hypercube with side \( L \) and dimension \( d \). This cube contains hence \( N=(L/h)^d \) points and therefore the error in the result scales as \( N^{-k/d} \) for the traditional methods.
• The error in the Monte carlo integration is however independent of \( d \) and scales as \( \sigma\sim 1/\sqrt{N} \), always!

* Comparing this with traditional methods, shows that Monte Carlo integration is more efficient than an order-k algorithm when \( d>2k \)

Example: Acceptance-Rejection Method

 
$$ p(x) \le M \quad x\in [a,b]. $$

 
Then we generate a random number \( x \) from the uniform distribution for \( x\in [a,b] \) and a corresponding number \( s \) for the uniform distribution between \( [0,M] \). If

 
$$ p(x) \ge s, $$

 
we accept the new value of \( x \), else we generate again two new random numbers \( x \) and \( s \) and perform the test in the latter equation again.

Acceptance-Rejection Method

 
$$ I=\int_0^3\exp{(x)}dx. $$

 
Obviously to derive it analytically is much easier, however the integrand could pose some more difficult challenges. The aim here is simply to show how to implent the acceptance-rejection algorithm. The integral is the area below the curve \( f(x)=\exp{(x)} \). If we uniformly fill the rectangle spanned by \( x\in [0,3] \) and \( y\in [0,\exp{(3)}] \), the fraction below the curve obatained from a uniform distribution, and multiplied by the area of the rectangle, should approximate the chosen integral. It is rather easy to implement this numerically, as shown in the following code.

Simple Plot of the Accept-Reject Method

Acceptance-Rejection Method

//   Loop over Monte Carlo trials n
     integral =0.;
     for ( int i = 1;  i <= n; i++){
//   Finds a random value for x in the interval [0,3]
          x = 3*ran0(&idum);
//   Finds y-value between [0,exp(3)]
          y = exp(3.0)*ran0(&idum);
//   if the value of y at exp(x) is below the curve, we accept
//   THIS IS OUR SAMPLING RULE
          if ( y  < exp(x)) s = s+ 1.0;
//   The integral is area enclosed below the line f(x)=exp(x)

//  Then we multiply with the area of the rectangle and
//  divide by the number of cycles
    Integral = 3.*exp(3.)*s/n

Monte Carlo Integration

 
$$ I=\int_0^1 f(x)dx\approx \frac{1}{N}\sum_{i=1}^Nf(x_i), $$

 
$$ I=\int_0^1 f(x)dx\approx E[f]=\langle f \rangle. $$

 
$$ \sigma^2_f=\frac{1}{N}\sum_{i=1}^Nf(x_i)^2- \left(\frac{1}{N}\sum_{i=1}^Nf(x_i)\right)^2, $$

 
or

 
$$ \sigma^2_f=E[f^2]-(E[f])^2=\left(\langle f^2\rangle - \langle f \rangle^2\right). $$

Brute Force Algorithm for Monte Carlo Integration

• Choose the number of Monte Carlo samples \( N \).
• Make a loop over \( N \) and for every step generate a random number \( x_i \) in the interval \( x_i\in [0,1] \) by calling a random number generator.
• Use this number to compute \( f(x_i) \).
• Find the contribution to the variance and the mean value for every loop contribution.
• After \( N \) samplings, compute the final mean value and the standard deviation

Brute Force Integration

// crude mc function to calculate pi
     int i, n;
     long idum;
     double crude_mc, x, sum_sigma, fx, variance;
     cout << "Read in the number of Monte-Carlo samples" << endl;
     cin >> n;
     crude_mc = sum_sigma=0. ; idum=-1 ;
//    evaluate the integral with the a crude Monte-Carlo method
      for ( i = 1;  i <= n; i++){
           x=ran0(&idum);
           fx=func(x);
           crude_mc += fx;
           sum_sigma += fx*fx;

      crude_mc = crude_mc/((double) n );
      sum_sigma = sum_sigma/((double) n );
      variance=sum_sigma-crude_mc*crude_mc;

Or: another Brute Force Integration

// crude mc function to calculate pi
int main()
{
  const int n = 1000000;
  double x, fx, pi, invers_period, pi2;
  int i;
  invers_period = 1./RAND_MAX;
  srand(time(NULL));
  pi = pi2 = 0.;
  for (i=0; i<n;i++)
    {
      x = double(rand())*invers_period;
      //   This is our sampling rule, all points accepted
      fx = 4./(1+x*x);
      pi += fx;
      pi2 += fx*fx;

  pi /= n;  pi2 = pi2/n - pi*pi;
  cout << "pi=" << pi << " sigma^2=" << pi2 << endl;
  return 0;

Brute Force Integration

     long idum;
     idum=-1 ;
     .....
     x=ran0(&idum);
     ....

or

  ...
  invers_period = 1./RAND_MAX;
  srand(time(NULL));
  ...
  x = double(rand())*invers_period;

Algorithm for Monte Carlo Integration, Results

Transformation of Variables

 
$$ p(x)dx=\left\{\begin{array}{cc} dx & 0 \le x \le 1\\ 0 & else\end{array}\right. $$

 
with \( p(x)=1 \) and satisfying

 
$$ \int_{-\infty}^{\infty}p(x)dx=1. $$

 
All random number generators provided in the program library generate numbers in this domain.

When we attempt a transformation to a new variable \( x\rightarrow y \) we have to conserve the probability

 
$$ p(y)dy=p(x)dx, $$

 
which for the uniform distribution implies

 
$$ p(y)dy=dx. $$

Transformation of Variables

 
$$ x(y)=\int_0^y p(y')dy', $$

 
which is nothing but the cumulative distribution of \( p(y) \), i.e.,

 
$$ x(y)=P(y)=\int_0^y p(y')dy'. $$

This is an important result which has consequences for eventual improvements over the brute force Monte Carlo.

Example 1, a general Uniform Distribution

 
$$ p(y)dy=\left\{\begin{array}{cc} \frac{dy}{b-a} & a \le y \le b\\ 0 & else\end{array}\right. $$

 
If we wish to relate this distribution to the one in the interval \( x \in [0,1] \) we have

 
$$ p(y)dy=\frac{dy}{b-a}=dx, $$

 
and integrating we obtain the cumulative function

 
$$ x(y)=\int_a^y \frac{dy'}{b-a}, $$

 
yielding

 
$$ y=a+(b-a)x, $$

 
a well-known result!

Example 2, from Uniform to Exponential

 
$$ p(y)=e^{-y}, $$

 
which is the exponential distribution, important for the analysis of e.g., radioactive decay. Again, \( p(x) \) is given by the uniform distribution with \( x \in [0,1] \), and with the assumption that the probability is conserved we have

 
$$ p(y)dy=e^{-y}dy=dx, $$

 
which yields after integration

 
$$ x(y)=P(y)=\int_0^y \exp{(-y')}dy'=1-\exp{(-y)}, $$

 
or

 
$$ y(x)=-ln(1-x). $$

 
This gives us the new random variable \( y \) in the domain \( y \in [0,\infty) \) determined through the random variable \( x \in [0,1] \) generated by our favorite random generator.

Example 2, from Uniform to Exponential

 
$$ I=\int_0^{\infty}F(y)dy=\int_0^{\infty}\exp{(-y)}G(y)dy $$

 
which we rewrite as

 
$$ \int_0^{\infty}\exp{(-y)}G(y)dy= \int_0^{\infty}\frac{dx}{dy}G(y)dy\approx \frac{1}{N}\sum_{i=1}^NG(y(x_i)), $$

 
where \( x_i \) is a random number in the interval [0,1].

Note that in practical implementations, our random number generators for the uniform distribution never return exactly 0 or 1, but we we may come very close. We should thus in principle set \( x\in (0,1) \).

Example 2, from Uniform to Exponential

.....
idum=-1;
x=ran0(&idum);
y=-log(1.-x);
.....

Example 3

 
$$ p(y)dy=\frac{dy}{(a+by)^n}, $$

 
with \( n > 1 \). It is normalizable, positive definite, analytically integrable and the integral is invertible, allowing thereby the expression of a new variable in terms of the old one. The integral

 
$$ \int_0^{\infty} \frac{dy}{(a+by)^n}=\frac{1}{(n-1)ba^{n-1}}, $$

 
gives

 
$$ p(y)dy=\frac{(n-1)ba^{n-1}}{(a+by)^n}dy, $$

 
which in turn gives the cumulative function

 
$$ x(y)=P(y)=\int_0^y \frac{(n-1)ba^{n-1}}{(a+bx)^n}dy'=, $$

 
resulting in

 
$$ y=\frac{a}{b}\left((1-x)^{-1/(n-1)}-1\right). $$

Example 4, from Uniform to Normal

 
$$ g(x,y)=\exp{(-(x^2+y^2)/2)}dxdy. $$

 
it is rather difficult to find an inverse since the cumulative distribution is given by the error function \( erf(x) \).

If we however switch to polar coordinates, we have for \( x \) and \( y \)

 
$$ r=\left(x^2+y^2\right)^{1/2} \quad \theta =tan^{-1}\frac{x}{y}, $$

 
resulting in

 
$$ g(r,\theta)=r\exp{(-r^2/2)}drd\theta, $$

 
where the angle \( \theta \) could be given by a uniform distribution in the region \( [0,2\pi] \). Following example 1 above, this implies simply multiplying random numbers \( x\in [0,1] \) by \( 2\pi \).

Example 4, from Uniform to Normal

 
$$ u=\frac{1}{2}r^2, $$

 
and define a PDF

 
$$ \exp{(-u)}du, $$

 
with \( u\in [0,\infty) \). Using the results from example 2, we have that

 
$$ u=-ln(1-x'), $$

 
where \( x' \) is a random number generated for \( x'\in [0,1] \). With

 
$$ x=r\cos(\theta)=\sqrt{2u}\cos(\theta), $$

 
and

 
$$ y=r\sin(\theta)=\sqrt{2u}\sin(\theta), $$

 
we can obtain new random numbers \( x,y \) through

 
$$ x=\sqrt{-2ln(1-x')}\cos(\theta), $$

 
and

 
$$ y=\sqrt{-2ln(1-x')}\sin(\theta), $$

 
with \( x'\in [0,1] \) and \( \theta \in 2\pi [0,1] \).

Example 4, from Uniform to Normal

.....
idum=-1;
radius=sqrt(-2*ln(1.-ran0(idum)));
theta=2*pi*ran0(idum);
x=radius*cos(theta);
y=radius*sin(theta);
.....

Box-Mueller Method for Normal Deviates

// random numbers with gaussian distribution
double gaussian_deviate(long * idum)
{
  static int iset = 0;
  static double gset;
  double fac, rsq, v1, v2;
  if ( idum < 0) iset =0;
  if (iset == 0) {
    do {
      v1 = 2.*ran0(idum) -1.0;
      v2 = 2.*ran0(idum) -1.0;
      rsq = v1*v1+v2*v2;
    } while (rsq >= 1.0 || rsq == 0.);
    fac = sqrt(-2.*log(rsq)/rsq);
    gset = v1*fac;
    iset = 1;
    return v2*fac;
  } else {
    iset =0;
    return gset;

Importance Sampling, chapter 11.4

Let us assume that \( p(y) \) is a PDF whose behavior resembles that of a function \( F \) defined in a certain interval \( [a,b] \). The normalization condition is

 
$$ \int_a^bp(y)dy=1. $$

 
We can rewrite our integral as

 
$$ \begin{equation} I=\int_a^b F(y) dy =\int_a^b p(y)\frac{F(y)}{p(y)} dy. \tag{27} \end{equation} $$

Importance Sampling

 
$$ x(y)=\int_a^y p(y')dy', $$

 
where we used

 
$$ p(x)dx=dx=p(y)dy. $$

 
If we can invert \( x(y) \), we find \( y(x) \) as well.

Importance Sampling

 
$$ I=\int_a^b p(y)\frac{F(y)}{p(y)} dy=\int_a^b\frac{F(y(x))}{p(y(x))} dx, $$

 
meaning that a Monte Carlo evalutaion of the above integral gives

 
$$ \int_a^b\frac{F(y(x))}{p(y(x))} dx= \frac{1}{N}\sum_{i=1}^N\frac{F(y(x_i))}{p(y(x_i))}. $$

 
The advantage of such a change of variables in case \( p(y) \) follows closely \( F \) is that the integrand becomes smooth and we can sample over relevant values for the integrand. It is however not trivial to find such a function \( p \). The conditions on \( p \) which allow us to perform these transformations are

• \( p \) is normalizable and positive definite,
• it is analytically integrable and
• the integral is invertible, allowing us thereby to express a new variable in terms of the old one.

Importance Sampling

Use the uniform distribution to find the random variable \( y \) in the interval [0,1]. \( p(x) \) is a user provided PDF.

Evaluate thereafter

 
$$ I=\int_a^b F(x) dx =\int_a^b p(x)\frac{F(x)}{p(x)} dx, $$

 
by rewriting

 
$$ \int_a^b p(x)\frac{F(x)}{p(x)} dx = \int_a^b\frac{F(x(y))}{p(x(y))} dy, $$

 
since

 
$$ \frac{dy}{dx}=p(x). $$

Perform then a Monte Carlo sampling for

 
$$ \int_a^b\frac{F(x(y))}{p(x(y))} dy,\approx \frac{1}{N}\sum_{i=1}^N\frac{F(x(y_i))}{p(x(y_i))}, $$

 
with \( y_i\in [0,1] \),

Evaluate the variance

Demonstration of Importance Sampling

 
$$ I=\int_0^1 F(x) dx = \int_0^1 \frac{1}{1+x^2} dx = \frac{\pi}{4}. $$

We choose the following PDF (which follows closely the function to integrate)

 
$$ p(x)=\frac{1}{3}\left(4-2x\right) \quad \int_0^1p(x)dx=1, $$

 
resulting

 
$$ \frac{F(0)}{p(0)}=\frac{F(1)}{p(1)}=\frac{3}{4}. $$

 
Check that it fullfils the requirements of a PDF. We perform then the change of variables (via the Cumulative function)

 
$$ y(x)=\int_0^x p(x')dx'=\frac{1}{3}x\left(4-x\right), $$

 
or

 
$$ x=2-\left(4-3y\right)^{1/2} $$

 
We have that when \( y=0 \) then \( x=0 \) and when \( y=1 \) we have \( x=1 \).

Simple Code

//   evaluate the integral with importance sampling
     for ( int i = 1;  i <= n; i++){
       x = ran0(&idum);  // random numbers in [0,1]
       y = 2 - sqrt(4-3*x);  // new random numbers
       fy=3*func(y)/(4-2*y); // weighted function
       int_mc += fy;
       sum_sigma += fy*fy;

     int_mc = int_mc/((double) n );
     sum_sigma = sum_sigma/((double) n );
     variance=(sum_sigma-int_mc*int_mc);

Test Runs and Comparison with Brute Force for \( \pi=3.14159 \)

\begin{tabular}{rllll}\hline $N$&$I_{cr}$ &$\sigma_{cr}$ &$I_{is}$ &$\sigma_{is}$\\\hline 10000 & 3.13395E+00 & 4.22881E-01 & 3.14163E+00 & 6.49921E-03 \\ 100000 & 3.14195E+00& 4.11195E-01 & 3.14163E+00 & 6.36837E-03 \\ 1000000& 3.14003E+00 & 4.14114E-01 & 3.14128E+00& 6.39217E-03\\ 10000000 &3.14213E+00 & 4.13838E-01 & 3.14160E+00 & 6.40784E-03 \\ \hline \end{tabular}

However, it is unfair to study one-dimensional integrals with MC methods!

Multidimensional Integrals

 
$$ I=\int_0^1dx_1\int_0^1dx_2\dots \int_0^1dx_d g(x_1,\dots,x_d), $$

 
with \( x_i \) defined in the interval \( [a_i,b_i] \) we would typically need a transformation of variables of the form

 
$$ x_i=a_i+(b_i-a_i)t_i, $$

 
if we were to use the uniform distribution on the interval \( [0,1] \). In this case, we need a Jacobi determinant

 
$$ \prod_{i=1}^d (b_i-a_i), $$

 
and to convert the function \( g(x_1,\dots,x_d) \) to

 
$$ g(x_1,\dots,x_d)\rightarrow g(a_1+(b_1-a_1)t_1,\dots,a_d+(b_d-a_d)t_d). $$

Example: 6-dimensional Integral

 
$$ \int_{-\infty}^{\infty}{\bf dxdy}g({\bf x, y}), $$

 
where

 
$$ g({\bf x, y})=\exp{(-{\bf x}^2-{\bf y}^2-({\bf x-y})^2/2)}, $$

 
with \( d=6 \).

Example: 6-dimensional Integral

 
$$ \frac{1}{\sqrt{\pi}}\exp{(-x^2)}, $$

 
and through

 
$$ \pi^3\int\prod_{i=1}^6\left( \frac{1}{\sqrt{\pi}}\exp{(-x_i^2)}\right) \exp{(-({\bf x-y})^2/2)}dx_1.\dots dx_6, $$

 
we can rewrite our integral as

 
$$ \int f(x_1,\dots,x_d)F(x_1,\dots,x_d)\prod_{i=1}^6dx_i, $$

 
where \( f \) is the gaussian distribution.

Brute Force I

.....
//   evaluate the integral without importance sampling
//   Loop over Monte Carlo Cycles
     for ( int i = 1;  i <= n; i++){
//   x[] contains the random numbers for all dimensions
       for (int j = 0; j< 6; j++) {
           x[j]=-length+2*length*ran0(&idum);

       fx=brute_force_MC(x);
       int_mc += fx;
       sum_sigma += fx*fx;

     int_mc = int_mc/((double) n );
     sum_sigma = sum_sigma/((double) n );
     variance=sum_sigma-int_mc*int_mc;
......

Brute Force II

double  brute_force_MC(double *x)
{
   double a = 1.; double b = 0.5;
// evaluate the different terms of the exponential
   double xx=x[0]*x[0]+x[1]*x[1]+x[2]*x[2];
   double yy=x[3]*x[3]+x[4]*x[4]+x[5]*x[5];
   double xy=pow((x[0]-x[3]),2)+pow((x[1]-x[4]),2)+pow((x[2]-x[5]),2);
   return exp(-a*xx-a*yy-b*xy);

Importance Sampling I

..........
//   evaluate the integral with importance sampling
     for ( int i = 1;  i <= n; i++){
//   x[] contains the random numbers for all dimensions
       for (int j = 0; j < 6; j++) {
	 x[j] = gaussian_deviate(&idum)*sqrt2;

       fx=gaussian_MC(x);
       int_mc += fx;
       sum_sigma += fx*fx;

     int_mc = int_mc/((double) n );
     sum_sigma = sum_sigma/((double) n );
     variance=sum_sigma-int_mc*int_mc;
.............

Importance Sampling II

// this function defines the integrand to integrate

double  gaussian_MC(double *x)
{
   double a = 0.5;
// evaluate the different terms of the exponential
   double xy=pow((x[0]-x[3]),2)+pow((x[1]-x[4]),2)+pow((x[2]-x[5]),2);
   return exp(-a*xy);
} // end function for the integrand

Test Runs for six-dimensional Integral

\begin{tabular}{rll}\hline $N$&$I_{cr}$&$I_{gd}$\\\hline 10000 & 1.15247E+01& 1.09128E+01 \\ 100000 & 1.29650E+01 & 1.09522E+01 \\ 1000000& 1.18226E+01 & 1.09673E+01 \\ 10000000&1.04925E+01 & 1.09612E+01 \\ \hline \end{tabular}

Overview of week 45

Overview of week 45

Monte Carlo Keywords

• Be able to generate random variables following a given probability distribution function PDF
• Find a probability distribution function (PDF).
• Sampling rule for accepting a move
• Compute standard deviation and other expectation values
• Techniques for improving errors Enhances algorithmic thinking!

Example, six-dimensional integral (project 3 2012)

 
$$ {\bf r}_i = x_i {\bf e}_x + y_i {\bf e}_y +z_i {\bf e}_z , $$

 
as

 
$$ \psi_{1s}({\bf r}_i) = e^{-\alpha r_i}, $$

 
where \( \alpha \) is a parameter and

 
$$ r_i = \sqrt{x_i^2+y_i^2+z_i^2}. $$

 
We will fix \( \alpha=2 \), which should correspond to the charge of the helium atom \( Z=2 \).

Switch to spherical coordinates

 
$$ d{\bf r}_1d{\bf r}_2 = r_1^2dr_1 r_2^2dr_2 dcos(\theta_1)dcos(\theta_2)d\phi_1d\phi_2, $$

 
and

 
$$ \frac{1}{r_{12}}= \frac{1}{\sqrt{r_1^2+r_2^2-2r_1r_2cos(\beta)}} $$

 
with

 
$$ \cos(\beta) = \cos(\theta_1)\cos(\theta_2)+\sin(\theta_1)\sin(\theta_2)\cos(\phi_1-\phi_2)) $$

How do I do importance sampling in spherical coordinates

\( \theta_{1,2} \in [0,\pi] \), use mapping \( \theta_{1,2}=\pi*ran \) with \( ran \in [0,1] \) a uniform distribution point.

\( \phi_{1,2} \in [0,2\pi] \), use mapping \( \phi_{1,2}=2\pi*ran \) with \( ran \in [0,1] \) a uniform distribution point.

Be careful with the integrand

 
$$ \frac{\exp{(-4(r_1+r_2))}r_1^2dr_1 r_2^2dr_2 dcos(\theta_1)dcos(\theta_2)d\phi_1d\phi_2}{\sqrt{r_1^2+r_2^2-2r_1r_2cos(\beta)}} $$

Example, six-dimensional integral (project 3 2012)

 
$$ \Psi({\bf r}_1,{\bf r}_2) = e^{-\alpha (r_1+r_2)}. $$

 
Note that it is not possible to find an analytic solution to Schr\"odinger's equation for two interacting electrons in the helium atom.

The integral we need to solve is the quantum mechanical expectation value of the correlation energy between two electrons, namely

 
$$ \begin{equation}label{eq:correlationenergy} \langle \frac{1}{|{\bf r}_1-{\bf r}_2|} \rangle = \int_{-\infty}^{\infty} d{\bf r}_1d{\bf r}_2 e^{-2\alpha (r_1+r_2)}\frac{1}{|{\bf r}_1-{\bf r}_2|}=\frac{5\pi^2}{16^2}=0.192765711. \end{equation} $$

 
Note that our wave function is not normalized. There is a normalization factor missing, but for this project we don't need to worry about that.

Example, six-dimensional integral (project 3 2012)

• a-b) Use Gaussian quadrature and compute the integral by integrating

for each variable \( x_1,y_1,z_1,x_2,y_2,z_2 \) from \( -\infty \) to \( \infty \). How many mesh points do you need before the results converges at the level of the fourth leading digit? Hint: the single-particle wave function \( e^{-\alpha r_i} \) is more or less zero at \( r_i \approx ? \). You can therefore replace the integration limits \( -\infty \) and \( \infty \) with \( -\Lambda \) and \( \Lambda \), respectively. You need to check that this approximation is satisfactory.

• c) Compute the same integral but now with brute force Monte Carlo

and compare your results with those from the previous point. Discuss the differences. With bruce force we mean that you should use the uniform distribution.

• d) Improve your brute force Monte Carlo calculation by using importance sampling. Hint: use the exponential distribution.

Does the variance decrease? Does the CPU time used compared with the brute force Monte Carlo decrease in order to achieve the same accuracy? Comment your results.

Example, six-dimensional integral, Gauss-Legendre

     double *x = new double [N];
     double *w = new double [N];
//   set up the mesh points and weights
     gauleg(a,b,x,w, N);

//   evaluate the integral with the Gauss-Legendre method
//   Note that we initialize the sum
     double int_gauss = 0.;
     for (int i=0;i<N;i++){
	     for (int j = 0;j<N;j++){
	     for (int k = 0;k<N;k++){
	     for (int l = 0;l<N;l++){
	     for (int m = 0;m<N;m++){
	     for (int n = 0;n<N;n++){
        int_gauss+=w[i]*w[j]*w[k]*w[l]*w[m]*w[n]
       *int_function(x[i],x[j],x[k],x[l],x[m],x[n]);
     		}}}}}
	}

Example, six-dimensional integral, Gauss-Legendre

//  this function defines the function to integrate
double int_function(double x1, double y1, double z1,
                    double x2, double y2, double z2)
{
   double alpha = 2.;
// evaluate the different terms of the exponential
   double exp1=-2*alpha*sqrt(x1*x1+y1*y1+z1*z1);
   double exp2=-2*alpha*sqrt(x2*x2+y2*y2+z2*z2);
   double deno=sqrt(pow((x1-x2),2)+pow((y1-y2),2)+pow((z1-z2),2));
  if(deno <pow(10.,-6.)) { return 0;}
  else return exp(exp1+exp2)/deno;
} // end of function to evaluate

Example, six-dimensional integral, brute force MC

     double int_mc = 0.;  double variance = 0.;
     double sum_sigma= 0. ; long idum=-1 ;
     double length=1.5; // we fix the max size of the box to L=3
     double volume=pow((2*length),6.);

//   evaluate the integral with importance sampling
     for ( int i = 1;  i <= n; i++){
//   x[] contains the random numbers for all dimensions
       for (int j = 0; j< 6; j++) {
           // Maps U[0,1] to U[-L,L]
           x[j]=-length+2*length*ran0(&idum);

       fx=brute_force_MC(x);
       int_mc += fx;
       sum_sigma += fx*fx;

     int_mc = int_mc/((double) n );
     sum_sigma = sum_sigma/((double) n );
     variance=sum_sigma-int_mc*int_mc;
     ....

Example, six-dimensional integral, brute force MC

double brute_force_MC(double *x)
{
   double alpha = 2.;
// evaluate the different terms of the exponential
   double exp1=-2*alpha*sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);
   double exp2=-2*alpha*sqrt(x[3]*x[3]+x[4]*x[4]+x[5]*x[5]);
   double deno=sqrt(pow((x[0]-x[3]),2)
          +pow((x[1]-x[4]),2)+pow((x[2]-x[5]),2));
   double value=exp(exp1+exp2)/deno;
	return value;
} // end function for the integrand