Numerical Methods for Physicists

by

Volker Hohmann
Institute of Physics
University of Oldenburg, Germany

 Dr. V. Hohmann, University of Oldenburg

http://www.physik.uni-oldenburg.de
Version 1.0, April 2004

1
Contents

1 INTRODUCTION 5

1.1 Definition: Computer physics 5

2 MATLAB – FIRST STEPS 6

2.1 Using Matlab 6

2.2 Help 7

2.3 Variables 7

2.4 Operators 8

2.5 Loops and conditions 9

2.6 Generating matrices/vectors 10

2.7 Functions 10

2.8 Visualization 11

2.9 Input and output 12

2.10 Writing individual commands (“M-Files“) 12

2.11 Exercises 14

3 REPRESENTATION OF NUMBERS AND NUMERICAL ERRORS 15

3.1 Error measures 15

3.2 Representation of numbers and roundoff errors 15

3.2.1 Integers 16
3.2.2 Floating point numbers 17

3.3 Range errors 18

3.4 Truncation errors 18

3.5 Error propagation in arithmetic operations 18

3.6 Error propagation in iterated algorithms 19

3.7 Exercises 21

4 NUMERICAL DIFFERENTIATION AND INTEGRATION 22

4.1 Differentiation 22
4.1.1 Right-hand formula (”naive“ ansatz) 22
4.1.2 Centered formula 23
4.1.3 Richardson extrapolation 24

4.2 Integration 26
4.2.1 Trapezoidal rule 26
4.2.2 Recursive trapezoidal rule 27
2
 4.2.3 Error analysis and Romberg integration 27
4.2.4 Gaussian integration 28

4.3 Exercises 29

5 ORDINARY DIFFERENTIAL EQUATIONS (ODE) 30

5.1 Mathematical formulation of the initial value problem 30

5.2 Simple methods 31

5.2.1 Euler’s method 31
5.2.2 Euler-Cromer method 32
5.2.3 Midpoint method 32
5.2.4 Verlet method 33

5.3 Global and local truncation error 33

5.4 Runge-Kutta method of the 4th order 33

5.4.1 Adaptive step size selection 35

5.5 Application of various methods to solve initial value problems 36

5.5.1 Projectile motion 36
5.5.2 Physical pendulum 37
5.5.3 Planetary motion 39
5.5.4 Lorenz model 40

5.6 Boundary value problems 40

5.7 Exercises 41

6 SOLVING SYSTEMS OF EQUATIONS 43

6.1 Linear systems of equations 43

6.1.1 Elementary operations 44
6.1.2 Gaussian elimination with backsubstitution 44
6.1.3 Gauss-Jordan method 46
6.1.4 Singular value decomposition 47

6.2 Nonlinear systems of equations 50

6.2.1 Newton’s method 50

6.3 Exercises 51

7 MODELING OF DATA 52

7.1 General mathematical formulation of the fitting problem 52

7.2 Establishing the merit function: Maximum-Likelihood Methods 53

7.2.1 Least-Squares method 53
7.2.2 Chi-Square method 54
7.2.3 Robust Fit methods 54

7.3 Minimization of functions 55

7.3.1 ”General Least-Squares“ method 55
7.3.2 Methods for nonlinear fitting: The ”Simplex“ method 57

7.4 Error estimation and confidence limits 59

7.4.1 Chi-square fitting: Confidence regions of parameters 59
7.4.2 The “Bootstrap“ method for estimating the confidence regions 60
3
7.5 Goodness of fit 60

7.6 Summary of methods 62

7.7 Exercises 63

8 ANALOG-DIGITAL TRANSFORM AND DISCRETE FOURIER TRANSFORM 65

8.1 Analog-Digital Transform (A/D Transform) 65

8.2 Diskrete Fourier Transform (DFT) 66

8.2.1 Real-valued time functions 68
8.2.2 Intensity spectrum 68

8.3 Fast Fourier Transform (FFT) 69

8.4 Filtering and convolution theorem 70

8.4.1 Convolution and cyclic convolution 71
8.4.2 Deconvolution 72

8.5 Exercises 73

9 PARTIAL DIFFERENTIAL EQUATIONS 74

9.1 Classification of partial differential equations 74

9.1.1 Parabolic equations 74
9.1.2 Hyperbolic equations 74
9.1.3 Elliptic equations 74

9.2 Initial value problems 75

9.2.1 Discretization 75
9.2.2 Heat equation: FTCS scheme 76
9.2.3 Time-dependent Schrödinger equation: Implicit Scheme (Crank-Nicholson) 76
9.2.4 Advection equation: Lax-Wendroff scheme 78
9.2.5 von-Neumann stability analysis 81

9.3 Boundary value problems: Poisson and Laplace equations 83

9.3.1 Laplace equation: Relaxation methods (Jacobi method) 83
9.3.2 Poisson equation 85

9.4 Exercises 89

4
1 Introduction
This lecture does not deal with the physics of computers but is rather about how to solve physical problems by
means of computers. Nowadays, computers play an important part as tools in experimental physics as well as in
theoretical physics. Special fields of use are:

1. Data logging
2. Systems control
3. Data processing (calculation of functions, visualization)
4. Transformation of theoretical models on the computer

Data logging and systems control cannot be done by hand and eye any longer and the aid of fast computers is
required for most experiments. Data processing enables large quantities of data to be analyzed and replaces
generations of ”human computers“, who were forced to calculate functional values from tables and to draw
diagrams by hand. Using computers for modeling is certainly the most complex way of employing computers in
physics, which is always necessary when analytical (mathematical) solutions do not exist. The three-body
problem, the nonlinear (physical) pendulum as well as solving the Schrödinger equation for complex atoms are
such examples. In these cases numerical methods are applied on a computer which seek the (approximate)
solution by calculating numbers according to simple recurrent rules (algorithms)1. The validity of such
numerical solutions and the applied numerical methods must always be examined very closely and it is advisable
to be suspicious about the results obtained in this way. This can be summarized in one sentence: It is insight not
numbers that we are interested in.

The objective of the lecture is to understand the fundamentals of these applications, especially points 3 and 4, by
means of simple examples and to acquire practical skills which can be applied or extended to concrete problems
later on. Even if the tasks are sometimes very complex: In the end it is always a matter of ’merely’ calculating
numbers according to simple rules.

1.1 Definition: Computer physics

The term ”computer physics” is to be defined on a mathematical basis:

1. Mathematics (e.g. analysis)

• is concerned with relatively abstract terms and uses symbols rather than numbers (numbers being also
symbols which are, however, described as an abstract quantity, e.g. group of natural numbers)
• Solutions are always analytical, i.e. determined with arbitrary precision in principle; exceptions from this
rule are dealt with theoretically (e.g. singularities)
• If analytical solutions are not known, proofs of existence or non-existence are investigated

2. Numerical mathematics
• deals with the calculation of numbers according to simple recurrent rules (algorithms)
• deals with questions of principle regarding the computability of numbers (convergence orders)
• Validation of algorithms on the basis of error estimates

3. Numerical mathematics on a computer

• Numerical mathematics + (unknown) loss of precision (e.g. by finite computational accuracy)

4. Computer physics
• Application of 3. to physical problems
• Additional possibility of validating algorithms on the basis of physical boundary conditions

1
Of course, the computer also helps to find analytical solutions; however, the respective computer algebra
programs (e.g. Maple or Mathematica) are based on the knowledge of the programmers, i.e. they do not find
solutions which are unknown in principle to the programmers.

5
2 Matlab – first steps
Matlab is a special program for numerical mathematics and is used throughout this course. It is easy to use and
allows us to rapidly enter the world of Numerics. Matlab is a very good complement to the known programs of
symbolic (analytical) mathematics (Maple, Mathematica) and is applied when no consistent analytical solutions
can be found for a problem. The fundamental data structure of Matlab are two-dimensional matrices (hence the
name: MATrix LABoratory). Matlab can carry out mathematical calculations in a syntax similar to that of C and
FORTRAN, respectively. It comprises all fundamental mathematical functions as well as special functions
(Bessel etc.). The essential algorithms known from numerical mathematics (solution of linear systems of
equations, eigenvalues/eigenvectors, differential equations etc.) have been implemented. Additionally, the
program includes options for graphically representing functions and data sets. It can generate two-dimensional
xy plots, plots in polar coordinates as well as two-dimensional representations. So-called toolboxes with further
functions are available for special applications (image processing, signal processing), which are not contained in
the standard package.

The following table shows the advantages and disadvantages of Matlab compared to programming languages
such as C or Fortran:

Advantages Disadvantages
fast program development possibly slower execution as compared to optimized
C programs
very good graphic options
many important algorithms implemented
portable programming (Linux, Windows)

Matlab version 6.x including a toolbox for signal processing is available in the CIP room of the Physics
Department and other public computer rooms on the campus. Since Matlab is licensed, it may not be transferred
to other computers! A login for the Windows network of the computing center of the university is required for
using Matlab (for details see: http://www.physik.uni-oldenburg.de/docs/cip/start.html)

In the following the fundamental structure of Matlab as well as the most important commands are explained.
Further below you will find some exercises, which will help you to practice Matlab.

2.1 Using Matlab

Upon selection of Matlab there appears a command line editor into which commands can be written, which are
executed by the program when the Return key (↵) has been pressed. Command lines executed before can be
recalled using the arrow keys (↑↓). Writing an initial and pressing the arrow keys produces only those previous
commands starting with that letter (or group of letters). The command line can be edited as known from text
processing programs (marking symbols/words with the mouse, replacing and inserting symbols/words etc.).

Matlab searches for the commands written in the command line in the directories stated in the path. With the
command “addpath” the path can be extended, if, for example, your personal directory with your own
commands is to be added to the path.

The command diary file name enables Matlab to remember all commands as well as the “answers” given by
Matlab in the stated file. Thus, the sequence of commands and results can be traced back afterwards.

Matlab stores the names and values of all variables generated during program execution in the so-called
“workspace“. With the commands save and load the workspace can be stored into a file and loaded
(reconstructed again later on to continue the task. With the command clear the variables in the workspace can
be deleted. The commands can be restricted to variable names by attaching variable names. For example, the
command save test.dat x y would only store the variables x and y in the file test.dat. The commands who and
whos list the variables available in the workspace, whos yielding more detailed information.

The graphs appear in additional windows on the screen. With the command figure a new graphic window is
generated, so that several graphs can be produced.

6
2.2 Help

In order to get help for a specific Matlab command write “help command”. Help offers a list of possible subjects
without further arguments. If the name of a command is not known, it is recommended to use the command
lookfor keyword. The keyword needs to be an English word or phrase. Its proper choice is crucial for the
succesful application of the lookfor-command.

Access to the entire documentation is obtained by the command helpdesk. The documentation then appears in
the HTML format. The manual ”Getting started“ which includes a detailed introduction to Matlab is of special
interest to beginners.

Furthermore, “Examples and Demos“ can be started from the help menu. There are examples arranged according
to subjects. Generally, the Matlab commands to be used are shown, such that the examples can be used as a
template for individual applications.

2.3 Variables

Names for variables may be assigned freely. Capitals and lower case letters are to be considered. Variables are
introduced by assignment and do not have to be declared. The type of variable (scalar, vector, matrix) is
determined by Matlab according to the assignment. Examples:

a = pi % the value Pi=3.1415... is assigned to a

b = [1 2] % the (row) vector (1, 2) is assigned to b
c = [1 ; 2] % the (column) vector (1, 2) is assigned to c
d = [[1 2 3];[4 5 6i]] % a 2X3 matrix is assigned to d

(Hints: The text following the symbol % is considered a comment.

The variable pi is predefined.
The variable i = − 1 is predefined for defining complex values).

The value of a variable can be displayed at any time by writing the name of the variable without any additions in
the command line. Examples:

a % input of variable
a = 3.14 % answer by MATLAB (scalar)

b % input of variable
b=1 2 % answer by MATLAB (line vector)

c % input of variable
c=1 % answer by MATLAB (column vector)
2

d % input of variable
d=123 % answer by MATLAB (2X3 matrix)
4 5 6i

7
Elements or sectors of matrices and vectors can also be selected and assigned to variables, respectively.
Examples:

d(1,3) % select the first tow, third column

d=3 % answer by MATLAB

d(1,2:3) % select the first, row, elements 2 to 3

d=23 % answer by MATLAB

z = d(2,1:3) % select row 2, columns 1 to 3 to z

z = 4 5 6i % answer by MATLAB

d(2,2) = 3 % assign the value 3 to the matrix d (row 2, column 2)

d=123 % answer by MATLAB (2X3 matrix)
4 3 6i

d(1,:) % address the first row

=123 % answer by MATLAB

d(1,1:2:3) % address the first row, every second element)

=1 3 % answer by MATLAB

d(1,[3 1]) % address the first row, third and first elements
=3 1 % answer by MATLAB

(Hints: Each operation causes Matlab to document its work on the screen. This can be stopped by typing a
semicolon (;) after a command. Indexing of elements of vectors/matrices always starts with index 1.) The term
‘end’ is a place holder for the highest index (last entry in a vector). E.g., a(1,end) always indexes the last column
in row 1.

2.4 Operators

The apostrophe or inverted comma operator ‘ (Attention: not to be confused with double quotation marks!)
performs a complex conjugation:

y = c’; % transform c into a row vector

y % show result
y=12 % answer by Matlab

y = d’; % calculate complex conjugated matrix

y % show result
y=1 4 % answer by MATLAB (3X2 matrix)
2 5
3 -6i

However, the operator .‘ performs a transposition:

y = d.’; % calculate transposed matrix

y % show result
y=1 4 % answer by MATLAB (3X2 matrix)
2 5
3 6i

The known operators +,-,*,/ function in Matlab just as in a pocket calculator. Applying these operators to
matrices and vectors, however, may lead to confusions, since Matlab analyses the variables to the left and right
of the operators and selects the operators 'intuitively', i.e. it selects scalar or matrix operations depending on the
dimension of the variables. Examples:

y = a * a; % multiply two scalars

y % show result

8
 y = 9.8696 % answer by Matlab

y = b * c; % multiply row by column vector

y % result is a scalar product
y=5

y = b .* c' % multiply two row vectors component by component

y=14 % result is a row vector again

y = d .^ d % exponentiating component by component

y=d^3 % exponentiating matrix

Similarly matrices can be multiplied by each other or by vectors, added etc. If an operation is to be done
component by component, a dot has always to be put before the operator. In any operation Matlab is very exact
about the dimension of matrices and vectors and gives the error message ''matrix dimensions must agree'' in case
of discrepancies.

2.5 Loops and conditions

Repetitive operations are performed in loops. The following loop calculates a row vector and a column vector:

for i=1:5 % start of loop

p(i) = i^2; % entries of a row vector
q(i,1) = i^3; % entries of a column vector
end % end of loop

Several loops can be interleaved. The index increment can be selected by the colon operator:

for i=1:2:5 % start of loop over uneven indices

q(i) = i; % entries for uneven indices
q(i+1) = -i; % entries for even indices
end % end of loop

When the command break is given within a loop, the loop is aborted. There are also loops which are carried out
as long as a certain condition is fulfilled:

x=100;
while( x > 1) % start of loop
x = x/2;
end % end of loop

(Hint: Loops are relatively slow in Matlab. If ever possible, they should be replaced by vector operations).

The following operators are available for establishing relations:

Operator Function
== even
~= uneven
> greater
>= greater than or
equal to
< smaller
<= smaller than or
equal to
& logical “and“
| logical “or“
~ logical “not“

9
The conditional execution of commands is done with so-called if-constructions.

if( x == 0 | x == 1 ) % start of if-construction

status = 1;
elseif( x < 0 & x ~= -1 ) % Attention: elseif must be written in one word
status = -1;
else
status = 0;
end % end of if-construction

2.6 Generating matrices/vectors

The colon operator (:) serves to address parts of matrices/vectors (cf. section Variables), but also to generate
matrices/vectors. Examples:

x = [1:10]; % generate row vector with integers between 1 and 10

x % show result
x = 1 2 3 4 5 6 7 8 9 10 % answer by Matlab

x = [1:3 ; 1 3 4]; % generate (2X3) matrix

x % show result
x=123 % answer by Matlab
134

A step size can also be selected:

x = [0:pi/4:pi]; % generate row vector with numbers from 0 to Pi,

% step size Pi/4.
x % show result
x = 0 0.7854 1.5708 2.3562 3.1416 % answer by Matlab

The following functions are available for generating matrices/vectors:

zeros - generates matrices/vectors and fills them up with zeros

ones - generates matricen/vectors and fills them up with ones
rand - generates matrices/vectors and fills them up with evenly distributed random
numbers in the interval (0,1)
randn - generates matrices/vectors and fills them up with Gaussian random numbers
(mean 0, variance 1);

The required numbers of lines and columns each are given as arguments. Example:

x = randn(2,3); % generate (2X3) matrix with random numbers (Gaussian)

2.7 Functions

Matlab offers all possible mathematical and numerical functions. All functions are applied to the entire
vector/matrix. For example, the following command generates a vector containing the values of the sinusoidal
function between 0 and 2 Pi with a resolution of Pi/10:

x = sin([0:pi/10:2*pi]);

The following table shows fundamental and special mathematical functions:

Function Description
sqrt(x) square root
sin(x) sine
asin(x) arc sine
sinh(x) hyperbolic sine

10
asinh(x) hyperbolic arc sine
cos(x) cosine
tan(x) tangent
exp(x) exponential function
expm(X) matrix exponential function
log(x) natural logarithm
logm(x) matrix logarithm
log10(x) decimal logarithm
rem(i,n) modulo function
round(x) rounding towards nearest integer
floor(x) rounding towards minus infinity
ceil(x) rounding towards plus infinity
inv(X) calculate the inverse of a matrix
fft(x) fast Fourier transform
abs(x) absolute value
real(x) real part
imag(x) imaginary part
angle(x) phase angle
sum sum of vector elements
prod product of vector elements
max search of maximum
min search of minimum
mean mean
std standard deviation
cov (co-)variance
median median
bessel Bessel functions
erf Gaussian error function
gamma gamma function
size calculate quantity of a matrix/vector
hist histogram
sort sorting algorithm

2.8 Visualization

Various visualization options are available to represent data graphically. The simplest one are xy-plots which
represent data points (x,y) in a cartesian coordinate system:

x = [0:pi/10:2*pi];
y = sin(x);
plot(x,y);

These commands generate a period of the sinusoidal function as a graph. The following table shows further
fundamental commands available for visualization:

Function Description
plot(x,y) xy plot, linear axes
loglog(x,y) xy plot, double logarithmic representation
semilogx(x,y) xy plot, logarithmic x-axis
semilogy(x,y) xy plot, logarithmic y-axis
title(‘text’) set image title
xlabel(‘text’) set lable of x-axis
ylabel(‘text’) set lable of y-axis
axis set value range of axes
polar(theta,rho) polar plot
contour(z) 2D representation of contour of matrix z
image(z) 2D color representation of matrix z

11
imagesc(z) like image, but with scaling of data
pcolor(z) similar to image(z)
colormap set color scale in case of color representation
plot3(x,y,z) 3D representation of vectors
contour3(z) 3D representation of contour of matrix z
mesh(z) 3D mesh plot of matrix z
surf(z) 3D mesh plot of matrix z with color
representation

In order to practice 2D and 3D representations a text matrix can be generated with the command peaks. For
example, the command imagesc(peaks) plots a 2-D color representation of the matrix generated by peaks.

2.9 Input and output

The command input reads an input from the keyboard:

x = input(‘write value of x: ‘)

This command displays the text on the screen and waits for an input from the keyboard. The input is assigned to
the variable x.

The command disp displays text and values of variables on the screen:

disp(‘the value of x is: ‘)

disp(x)

These commands display the text as well as the value of variables x.

The command fprintf is available for a formatted printout as known from the programming language C.
Additionally, the commands fopen, fread a well as fwrite and fclose are available for reading or writing in files
(Hint: not to be confused with the commands load and save, which store the workspace in a special Matlab data
format).

2.10 Writing individual commands (“M-Files“)

Individual programs and functions are easy to write in Matlab. For this purpose a so-called script is written,
which is a text file comprising Matlab commands. Matlab interprets the file name as the name of the command;
the extension “m” of the file is fixed. Therefore, the files are called M-Files. How such M-Files have to be
arranged is explained by an example in the following. Further information can be obtained by the Matlab Help
function. Write help script or help function.

(Hint: The directory in which the files are stored must be in the search path (command path, addpath)
or the files must be in the actual directory (command cd) so that Matlab can find the individual
M-Files.)

A script is produced by writing the commands into the text file line by line. A usual text editor (e.g. notepad) can
be used for this, but also the special Matlab text editor (call the editor by selecting Open or New in the menu
File). The script is started by writing the file name in the command line. Matlab then executes all commands and
keeps the variables generated upon execution of the commands within the workspace.

(Hint: When an M-File is changed, Matlab perhaps does not read it in anew. The old version is executed
(flaw in the network installation of Matlab, e.g. in the CIP room). The command clear all forces
the program to read in a new version which, however, produces the side-effect that all variables
in the workspace are deleted.

The script can be transformed into a function, which behaves like a special Matlab function. For this certain
formalities are to be considered, which are explained by an example in the following:

12
function [mean,stdev] = stat(x)
%
% [mean,stdev] = stat(x)
%
% This function calculates the mean and standard deviation of a data set
%
% Transfer parameter:
% x – data set as vector
%
% Return values:
% mean - mean
% stdev - standard deviation
%

n = length(x);
mean = sum(x) / n;
stdev = sqrt(sum((x - mean).^2)/n);

The first line begins with the keyword function and declares the name of the function (here: stat), the arguments
in parentheses (here: x) as well as the return values (here: mean and stdev). The following list of comments is
displayed by Matlab, if the help function is called for this command (help stat). Then the actual function starts
which can consist of Matlab commands. It is important that the return values are calculated from the arguments.
For this purpose subsidiary variables can be generated, too (here: n). The function controls its own workspace,
i.e. the function only knows the values of the arguments and only the return values are returned to the calling
program. All subsidiary variables are local and are deleted upon termination of the function. Now the function
can be called from the command line as follows:

noise = randn(1,100); % generate random numbers

[mnoise,snoise] = stat(noise); % calculate statistics
disp(‘mean and standard deviation are:’);
mnoise
snoise

13
2.11 Exercises

Remark: Save the solutions of all exercises in Script files ("M-Files").

I. Generating data

a. Generate a time vector from 0ms to 10 ms at a sampling frequency of 10 kHz.

b. Generate a sine with a frequency of 1 kHz over the time basis given by the time vector from a.
c. Plot the sine from b. over the time vector from a.
d. Like c., however, only every second and fourth sampling values are to be plotted.
e. Save the generated data in a data file (hint: use the command "save")
f. Replace the first example loop in the section Loops and conditions by vector operations.

II. Using matrices

 1 2 3
 
a. Define the matrix A =  4 5 6
 
 7 8 9
and the vectors b = 1( 2 3) , c = (4 5 6) and d = (1 2) as line vectors.

b. Calculate AbT ; bA ; b * c as scalar product and (!) component-wise.

c. Solve the linear system of equations Ax=bT (Hint: "help slash") and verify the solution by insertion.

 5 6
d. Calculate A2 dT , where A2 =   (write A2 as submatrix of A!).
 8 9

III. 2D and 3D representation of matrices

a. Generate a two-dimensional Gaussian bell with a variance of 8 on a (25x25) matrix.

b. Plot the matrix in a 2D color representation with contours.
c. like b, however, in a 3D color representation.
d. Optional task: How can a. be solved without a FOR loop?

The examples in the menu ‚Help->Examples and Demos‘ are very instructive and show many tricks and
knacks of Matlab programming (the source texts for all examples are included!). It is recommended to
look through these examples.

14
3 Representation of numbers and numerical errors
Types of errors:

1. Errors in the input data (e.g. finite precision of measuring devices)

2. Roundoff errors
• because of finite precision of the number representation on computers
• because of finite precision of arithmetic/mathematical calculations on the computer
3. Range errors
4. Truncations errors
5. Error propagation (instable solutions in iterated algorithms)

Type 1 will not be treated here, because the measurement process is not a matter of Numerics. Before treating
types 2-5 in detail, common error measures will be introduced.

3.1 Error measures

Definition: x* be an approximation to x (x,x*∈ A, A set of numbers, e.g. real numbers). Then

∆x = x − x *
x − x * ∆x
r ( x) = =
x x

are the absolute error and the absolute relative error of x*.

Mostly, ∆x is not known, but an error bound ε can be estimated, which indicates the maximum possible
absolute error:

x * −ε ≤ x ≤ x * + ε ⇔ x = x * ±ε

In most cases the ”real“ or “true” number is not known, but only the approximation x*. Then the relative error is
approximated as well:

∆x ε
r ( x) ≈ ≤
x* x*

Example: If the measured value is 0.0123 ± 0.0005, then ∆x ≤ ε = 0.0005 . The measured value has two
significant digits. The relative precision is r ( g ) ≤ 0.0005 0.0123 ≈ 0.04 .

3.2 Representation of numbers and roundoff errors

Digital computers cannot represent infinitely many different numbers, because only a finite number of
information units (mostly stated in bits2) are available for representing a number. The set A of numbers that can
be represented is called the set of machine numbers and generally is a subset of the real numbers. A depends on
the exact kind of number representation on the computer and implicitly involves a rounding operation rd(x)
stating how real numbers are approximated by machine numbers. It can be defined generally:

2
1 Bit corresponds to a ‘yes-no’ decision (0 or 1).

15
 rd : ℜ → A,
rd ( x ) ∈ A: x − rd ( x ) ≤ x − g ∀g ∈ A

The absolute and relative approximation errors caused by the finite number representation are obtained by
setting x* = rd(x) in the formulas given above. In general, the result of an arithmetic operation is not an element
of A, even if all operands are elements of A. The result must again be subjected to the rounding operation rd(x)
in order to represent it. Thus, arithmetic operations can lead to additional roundoff errors, even if the result of
the operation was calculated with infinite precision:

x, y ∈ A ⇒ rd ( x + y ) ∈ A , but ( x + y ) not necessarily ∈ A

Digital computers work with different kinds of number representations and thus with different machine
numbers. Special common representations of numbers as used in most digital computers are shown below.

3.2.1 Integers

Integers are mostly represented by the two’s complement, representing each machine number z (z being in
decimal format) by a set of bits as follows:

z {α0 ,α1 ,...,α N −1} , α n ∈ {0,1}

with
N −2
z = −α N −1 ∗ 2 N −1 + ∑ α n 2n with α n ∈ {0,1}
n =0

An information unit of 1 bit („0“ or „1“) is required for each coefficient αn , such that N bits are necessary to
represent the number (common quantities for N being 8,16,32, and 64). Hence, the bit representation of several
numbers for e.g. N=8 is:

Number Bit representation

+1 0 0 0 0 0 0 0 1
+2 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0
-1 1 1 1 1 1 1 1 1
+117 0 1 1 1 0 1 0 1
-118 1 0 0 0 1 0 1 0
Significance -27 26 25 24 23 22 21 20

Key parameters of the two’s complement are:

N −1
Range of numbers: −2 ≤ z ≤ 2 N −1 − 1

Error bound: ∆z ≤ ε = 0.5

∆x 0.5
Relative error: r ( z ) ≈ ≤
z* z*

The representation of integers has a constant absolute error and hence a relative error depending on the value of
the number to be represented.

Remark: The name ”two’s complement“ is explained by the fact that the transition from a positive to the
negative number of equal absolute value results from performing a bit by bit complement in the bit
representation (each “1” becomes “0” and vice versa) and then adding 1.

16
Remark: The advantage of the two’s complement over other conceivable representations of integers is the fact
that the addition of two numbers is done by simply adding bit by bit in the bit representation. It is not necessary
to differentiate between the cases of positive and negative numbers. Example (note the carry-over):

+2 0 0 0 0 0 0 1 0
-1 1 1 1 1 1 1 1 1
1 0 0 0 0 0 0 0 1

Remark: Due to the significance of each single bit in the bit representation multiplying by 2 means simply
shifting all bits to the left by one position (“left shift”) and accordingly, dividing by 2 means a “right shift“.

3.2.2 Floating point numbers

A machine number g (g being in decimal format) is represented in the floating point binary format as follows:

g = ( −1) ∗ m ∗ B e , 1 B ≤ m < 1
s

with:

s sign bit (”0“ or ”1“

m mantissa
e exponent
B basis (fixed value)

Remark: The side condition 1/B≤ m< 1 makes the representation unique. This is called the normalized
representation. For example, the number 0.5 with the basis B=10 can be represented as 0.5*100 or 0.05*101. The
first representation is fixed by the side condition.

The 32bit IEEE Floating-Point Format is used on most digital computers. B=2 is used as a basis and the
exponent e is represented as an 8 bit two’s complement. A 23bit representation with the following bit
significances is used for the mantissa m:

Bit m0 m1 m2 ... m21 m22

Significance 2-1 2-2 2-3 ... 2-22 2-23

Hint: The normalization condition causes the bit m0 of the mantissa m to be always 1. Therefore, it is generally
not saved at all.

Example: The number 3 can be represented as follows:

(
3 = 2 − 1 + 2 − 2 ∗2 2 )
Thus, the bit representation of the number 3 in the IEEE Floating-Point Format is:

s e0 e1 E2 e3 e4 e5 e6 e7 m0 m1 m2 ... m21 m22

0 0 0 0 0 0 0 1 0 1 1 0 ... 0 0

Key parameters of this representatuion are:

Range of numbers (approx.): −10 ≤ g ≤ 1038

38

Be g=m*Be ∈ A a machine number:

17
 1 −M
Error bound: ∆g ≤ ε = ∗ 2 ∗ B e , M : Number of bits of the mantissa
2

∆g ε ε
Relative error: r ( g ) ≈ ≤ ≤ = 2− M
g * m∗ B e
1 2∗ B e

The floating point representation has a constant relative error. The error bound is called machine precision.

Hint: The machine precision is often defined as the smallest number ε > 0, which still yields a result >1.0 upon
addition (1+ε).

Hint: Besides the 32bit IEEE Floating-Point Format (”float“), a corresponding 64bit IEEE Floating-Point
Format (”double“) is available in most programming languages, which possesses more bits for exponent and
mantissa and thus a higher precision and a larger range of numbers. Matlab internally works with the 64bit-
format.

3.3 Range errors

The different representations of numbers cover a limited range of numbers (see above). This range can easily be
exceeded by arithmetic operations. An important example is the calculation of the factorial n!, which exceeds
the range of numbers of the 64bit floating point representation already for n>170. The range error caused by
exceeding the range of numbers is bigger than the roundoff error by orders of magnitudes and must be avoided
by choosing the appropriate representation of numbers and by checking the orders of magnitude of the numbers
examined (Example: Calculation of the faculty via application of the Stirling formula and logarithmic
representation).

Furthermore, it is possible that certain physical constants exceed or fall below the range of representable
numbers (cf. Exercises).

3.4 Truncation errors

Errors inherent to the algorithm, even if no rounding errors occur. Occur in iterated algorithms,when the
iteration is stopped after a certain number of iterations (Examples: Calculation of the exponential function by
finite sums and approximation of the integral by rectangles of finite width).

3.5 Error propagation in arithmetic operations

Addition:

x = a + b ; a, b ≥ 0 ,
xmax = a + ∆a + b + ∆b
xmin = a − ∆a + b − ∆b
∆x = ∆a + ∆b
∆x ∆a + ∆b
r ( x) = =
x a+b

This means that the absolute errors of the input values to the addition are added.

18
Subtraction:

x = a − b ; a, b ≥ 0 ,
xmax = a + ∆a − b + ∆b
xmin = a − ∆a − b − ∆b
∆x = ∆a + ∆b
∆x ∆a + ∆b
r ( x) = =
x a −b

When two approximately equal numbers are subtracted, the relative error becomes very big! This is called loss
of significance, catastrophic roundoff error or subtractive cancellation.

Example: Loss of significance in calculating the square formulas (→ Exercises).

Example: Look at the difference of the numbers p=3.1415926536 and q=3.1415957341, which are about equal
and are both significant on 11 decimals. The difference p-q=-0.0000030805 has only five significant digits left,
although the operation itself did not cause new errors.

Multiplication:

x = a ∗ b ; a, b ≥ 0 ,
xmax = (a + ∆a ) ∗ (b + ∆b ) ≈ a ∗ b + ∆a ∗ b + a ∗ ∆b , da ∆a ∗ ∆b much smaller than ∆a ∗ b
xmin = (a − ∆a ) ∗ (b − ∆b ) ≈ a ∗ b − ∆a ∗ b − a ∗ ∆b
∆x ≈ ∆a ∗ b + a ∗ ∆b
∆x ∆a ∗ b + a ∗ ∆b ∆a ∆b
r ( x) = ≈ = +
x a *b a b

Division:

x = a b ; a, b ≥ 0 ,

xmax = (a + ∆a ) (b − ∆b ) =
(a + ∆a )* (b + ∆b ) ≈ a b + ∆a ∗ b + a ∗ ∆b
(b − ∆b )* (b + ∆b ) b2

xmin = (a − ∆a ) (b + ∆b ) =
(a − ∆a )* (b − ∆b ) ≈ a b − ∆a ∗ b + a ∗ ∆b
(b + ∆b )* (b − ∆b ) b2
∆a ∗ b + a ∗ ∆b
∆x ≈
b2
∆x ∆x ∆a ∗ b + a ∗ ∆b ∆a ∆b
r ( x) = = ≈ = +
x ab b2 ∗ a b a b

The absolute errors of the operands add up in addition and subtraction, while the relative errors add up in
multiplication and division.

3.6 Error propagation in iterated algorithms

Many algorithms contain iterated reproductions, in which a set of calculations is repeated with the data
calculated in the preceding step as long as a certain break condition is reached. The errors produced in one step
propagate from step to step. Depending on the problem, the total error may increase linearly, increase
exponentially or decrease exponentially (error damping). The behavior of an algorithm can be deduced from

19
the operations used in each iteration and the rules of error propagation in the fundamental arithmetic operation
explained above.

Definition: E(n) be the total error following n iteration steps and ε the machine precision. If |E(n)| ≈ nε, the error
propagation is called linear. If |E(n)| ≈ Knε, the error propagation is exponential, which leads to error damping
for K<1.

Example: Iterated reproductions lead to a strong dependence on the initial values in case of exponential error
propagation and thus represent a path into deterministic chaos. An example for this is the logistic parabola.

20
3.7 Exercises

Remark: Save the solutions of all tasks in Script files ("M-Files").

I Roundoff errors

1. Determine the machine precision of Matlab and estimate the number of bits used to represent the mantissa.
(Hint: Start at eps=1 and divide eps by 2 until (1+eps) is not larger than 1 any longer).
2. Calculate the formula x=(10+h)-10 for h=B-n for n=0...21 and B=10. Calculate the absolute relative error
between x and h and plot it as a function of n. Repeat the same for B = 2. Why do we find distinct
differences here?

II Range errors

1. How can the range error be avoided in calculations with physical constants (atom physics, quantum
mechanics)?
2. Take eV (electron volt) as a unit for energy and the mass of the electron as a unit for mass. Convert 1kg, 1m,
and 1 s into these units.

III Truncation errors

Taylor’s series of the exponential function reads:

N
xi
e x = lim S ( x, N ) with S ( x, N ) = ∑
N →∞
i =0 i!
Calculate the subtotals S(x,N) up to N=60 and plot the absolute relative error relating to the true value (Matlab
function exp()) as a function of N. Test the program for x=10,2,-2 and -10. Why does the error increase for
negative numbers?

IV Error propagation

1. The square equation ax + bx + c = 0 has the solutions:

2

− b + b 2 − 4ac − b − b 2 − 4ac
x1 = and x2 =
2a 2a
a. Under which conditions do we find a loss of decimal places in the numerator?
b. Prove that the zeros can also be calculated according to the following equivalent formula:
− 2c − 2c
x1 = and x 2 =
b + b − 4ac
2
b − b 2 − 4ac
c. How can the loss of decimal places be avoided with b.?
d. Calculate the zeros for a=1, b=-10000.0001, and c=1 as well as for a=1, b=-100000.00001, and c=1 with
both formulas. Compare the results with the correct solution and verify them by inserting into the equation
(Hint: Use the command sprintf, otherwise Matlab will not represent enough decimal places (e.g.
sprintf('%5.10g',x1) )).

21
4 Numerical differentiation and integration
Numerical methods of differentiation and integration are applied when analytical solutions are not known. There
are many examples of functions that are not integrable analytically, e.g. the Gaussian error function

x
2
erf ( x ) = ∫ e − y dy ,
2

π
0

hence, the integral over the Gaussian distribution. Contrary to integration, differentiation in the analytical way is
virtually always feasible.Consequently, numerical differentiation itself is rarely applied directly. However, it is
the basis for developing algorithms to solve ordinary and partial differential equations and is therefore described
in detail below.

The most important numerical methods of differentiation and integration are explained in the following. The
procedures of developing algorithms and of estimating numerical errors are generally applicable and will be
found again in a complex form in later chapters.

4.1 Differentiation

The principle of numerical differentiation is to approximate the function to be differentiated by a polynomial,

which can then be derivated according to the known rules for polynomials. The value of the derivative of the
polynomial is used as an estimated value for the derivative of the function. In general the approximation and the
estimated value for the derivative are the more precise, the higher the order of the polynomial. But this is not a
law and there are exceptions, especially when the function to be derivated cannot be well approximated by
polynomials of high orders.

4.1.1 Right-hand formula (”naive“ ansatz)

We know from analysis that the derivative is calculated according to the following equation:

f ( x + h) − f ( x )
f ' ( x ) = lim
h →0
h

The limiting process cannot be done numerically, however, selecting an appropriate small h yields an
approximation of the derivative according to the following equation:

f ( x + h) − f ( x )
f '(x) = + εa ( h > 0),
h

εa being the (unknown) truncation error resulting from the omission of the limiting process (finite h). The
truncation error occurs even if the quotient has been calculated with arbitrary precision. Numerical calculation of
the quotient leads to the following additional sources of errors:

1. Round-off errors in calculating x+h.

2. Round-off errors in calculating the numerator (Attention: Loss of significance is possible.)

22
In order to avoid the first kind of errors h should be chosen such that x+h is a number that can be exactly
represented3. The second rounding error cannot be avoided. If the function f is calculated to the machine
precision ε m , the absolute error in the calculation of the derivative is as follows:

2εm
∆( f ' ( x )) = + εa ( h > 0)
h
For calculating the truncation error we expand the function into a Taylor’s series of the first order:

1 2 ''
f ( x + h ) = f ( x ) + hf ' ( x ) + h f (ζ ) , x≤z≤x+h
2

ζ is a value within the observed interval according to Taylor’s theorem. The equation above can be
'
immediately derived by resolving for f ( x ) :

f ( x + h) − f ( x ) 1 ''
f '(x) = − hf (ζ )
h 2
f ( x + h) − f ( x )
= + O( h)
h
Hence, the truncation error is linear in h. This is called the order O(h). The absolute error can then be stated as
follows:

∆( f ' ( x )) =
2ε m 1
h
+ Mh ,
2
M= max
x ≤ζ ≤ x + h
( f (ζ ) )
''

The first addend increases with decreasing h, while the second decreases with decreasing h. Thus, there is an
optimum h with minimal error, which depends on the maximum M of the second derivative in the observed
interval. Since M is usually unknown, it appears to be rather pointless to choose h arbitrarily small because of
the numerical errors. It is appropriate to start with a relatively large h and then reduce it until the estimation of
the derivative does not change any longer referring to a predetermined precision.

4.1.2 Centered formula

It would be possible to improve the precision, if the truncation error decreased more than linearly with h, e.g. by
the order O(h2) or higher. The object of developing such formulas is to use Taylor’s series of higher orders and
to incorporate the interval x-h as well:

1 2 '' 1
f ( x + h ) = f ( x ) + hf ' ( x ) + h f ( x )+ h 3 f ''' (ζ1 ) , x ≤ ζ1 ≤ x + h
2 6
1 1
f ( x − h ) = f ( x ) − hf ' ( x ) + h 2 f '' ( x ) − h 3 f ''' (ζ2 ) , x-h ≤ ζ2 ≤ x + h
2 6

'
The two formulas are subtracted from each other and resolved for f ( x ) . In doing so the second derivative is
dropped (just like all further even orders would be dropped):

~
x =x+h ~
3
First calculate the following temporary quantities: ~ . x, ~
x and h are exactly representable
~
h =x−x
numbers, as they were explicitly assigned to a variable in the workspace. Then calculate the formula with x , ~
x
~
and h .

23
 f ( x + h ) − f ( x − h ) 1 2 f (ζ1 ) + f (ζ2 )
''' '''
'
f (x) = + h
2h 6 2
f ( x + h) − f ( x − h)
=
2h
+O h 2 ( )
In this so-called centered formula the truncation error decreases with h2, while the roundoff error is
proportional to1/h in calculating the quotient, as already shown for the right-hand formula. In general this leads
to a smaller total error.
An estimation of the second derivative can be found in a similar way. For this purpose the Taylor expansion is
''
again performed in the intervals x+h and x-h and then resolved for f ( x ) by adding the two series:

f ( x + h ) = f ( x ) + hf ' ( x ) +
1 2 ''
2
1
h f ( x )+ h 3 f ''' ( x )+O h 4
6
( )
1 1
( )
f ( x − h ) = f ( x ) − hf ' ( x ) + h 2 f '' ( x ) − h 3 f ''' ( x )+O h 4
2 6
( )
f ( x + h ) + f ( x − h ) = 2 f ( x ) + h 2 f '' ( x )+O h 4

f ( x + h) + f ( x − h) − 2 f ( x)
⇒ f '' ( x ) =
h2
( )
+O h 2

Hence, the second derivative is calculated to the order O(h2) using the centered formula. However, the function f
must be evaluated at three instead of two points.

4.1.3 Richardson extrapolation

Of course, formulas of higher order can be found in a way similar to that shown in the preceding section, but the
formulas of higher orders can be calculated from the formulas of lower orders more easily. This procedure is
called extrapolation and is described in the following. For this we observe the entire Taylor expansion of
function f at the point x. An appropriate conversion yields the centered formula again, however, the truncation
error is calculated in all powers of h:

∞
f (k ) (x ) k
f (x + h ) = ∑ h
k =0 k!
∞
f (k ) (x ) k
f (x − h ) = ∑ (−1)
k
h ⇒
k =0 k!
f (2k +1) (x ) (2k +1)
∞
f (x + h ) − f (x − h ) = 2∑ h ⇒
k = 0 (2k + 1) !

f (x + h ) − f (x − h ) ∞
f (2k +1) (x )
= f ' (x ) + ∑ αkh 2k , αk =
2h k =1 (2k + 1) !

The series represents the (unknown) truncation error. The quotient is defined as D0 ( h) 4 in the left part of the
last equation and is evaluated for the step sizes h and 2h:

4
D0 ( h) corresponds to the centered formula.

24
 ∞
D0 ( h) = f ' ( x ) + ∑ α k h2 k
k =1
∞
D0 (2h) = f ' ( x ) + ∑ 4 k α k h 2 k
k =1

The next higher power of h in the truncation error (k=1) can now be eliminated by appropriately
combining D0 ( h) and D0 ( 2h) :

∞ ∞
4D0 (h ) − D0 (2h ) = (4 − 1) f '
(x ) + 4 ∑ αkh − ∑ 4k αkh 2k
2k

k =2 k =2
∞
= 3 f ' (x ) + ∑ αkh 2k
k =2

The next iteration to the approximation of f

'
( x ) is thus obtained as follows:
4D0 (h ) − D0 (2h )
D1 (h ) ≡
3
D (h ) − D0 (2h )
= D0 (h ) + 0
3
∞
1
= f ' (x ) + ∑ αkh 2k
3 k =2

Since the series for the truncation error starts only at k=2, D1 ( h) has a truncation error of the order O(h4) and is
calculated from the two quantities D0 ( h) and D0 ( 2h) , which are only of the order O(h2) themselves. This can
be continued in general in order to eliminate higher orders:

4 k Dk −1 ( h) − Dk −1 ( 2h)
Dk ( h) =
4k − 1
D ( h) − D ( 2 h)
= Dk −1 ( h) + k −1 k k −1
4 −1
This recursion formula is called Richardson extrapolation and can be used whenever the truncation error
shows only even powers in h. The numerical differentiation is then obtained by calculating the following table:

D0 (h) . . .
D0 ( h 2) D1 (h 2) . .
D0 ( h 4) D1 (h 4) D2 ( h 4) .
D0 (h 8) D1 (h 8) D2 (h 8) D3 (h 8)

The table can only be calculated row by row, except for the first column which is directly yielded according to
the centered formula. The following columns then result from the extrapolation formula. The highest-order
result is then found as the last value in each row (diagonal element). Calculation of the rows is stopped as soon
as the row result does not change any longer (referring to a predetermined precision).

25
4.2 Integration
Consider the definite integral
b

∫ f ( x )dx
a
The essential technique for numerically calculating definite integrals is to approximate the function f in the
observed interval [a,b] by a polynomial P, which is then integrated according to the simple rules of polynomial
integration. The integral of the polynomial is then used as an estimate of the integral of the function:

N
f ( x ) [ a , b] ≈ P( x ) , P( x ) = ∑ αn x n
n=0
b b

∫ f ( x )dx ≈ ∫ P( x )dx
a a

The higher the order N of the polynomial, the better the assessment of the integral value in general. However,
this is only true for functions which can be approximated by polynomials of high orders (which is not always
possible, e.g. in case of discontinuities). Thus, the problem is reduced to performing an appropriate polynomial
approximation. Generally speaking, a polynomial of the N-th order can be determined from N+1 functional
values. For approximation of the function f by a polynomial of the N-th order it is sufficient to evaluate f at N+1
points in the observed interval [a,b].

4.2.1 Trapezoidal rule

Let us assume that we divide the interval [a,b] into N intervals. For this N+1 points within the interval need to be
defined:

x0 = a, x N = b, x 0 < x1 <… < x N −1 < x N

The function f is evaluated at these points: f i ≡ f x i , ( ) 0<i< N

The simplest method of integration is to linearly connect the functional values at these points piecewise. Thus,
we obtain a trapezoid with the following area in each interval:

1
Ti = ( x − xi −1 )( f i + f i −1 ),
2 i
1< i < N

N
The total area is then: T = ∑ Ti
i =1

This formula is simplified, if the points are evenly distributed in the interval. Be h = ( b − a ) N the width of
( )
each interval, hence x i = a + ih . The area of each interval is then Ti = 1 2 h f i + f i −1 . Except the
marginal points x0 and xN1, each f i occurs twice in the formula for the total area:

N −1  N 
T ( h ) = 1 2 h( f 0 + f N ) + h ∑ fi ≡


∑ wi f i 

i =1 i=0

According to this trapezoidal rule the integral can be approximated the better, the smaller h is (the occurring
truncation error due to the finite interval is stated later on) (Remark: The estimate of the integral results as a
weighted sum of the sampling points of the function).

26
4.2.2 Recursive trapezoidal rule
Now we want to iteratively calculate the aproximation of the integral according to the trapezoidal rule by
systematically reducing the interval h. It is useful to bisect h each time proceeding from the complete interval
(h=(b-a)):

1
hn = ( b − a ), n = 0,1, …
2n

The functional values f i calculated in the preceding iteration step can be used again in this case. By simple
consideration we obtain:

1
T ( h0 ) = h ( f (a ) + f (b))
2 0
2n
1
T ( hn +1 ) = T ( hn ) + hn +1
2 ∑ f (a + (2i − 1)hn +1 )
i =1

The series T hn ( ) converges towards the value of the integral.

4.2.3 Error analysis and Romberg integration
A procedure similar to that one demonstrated in section 4.1.3 shows that the truncation error of the trapezoidal
rule contains only even powers of h:

b ∞
T ( hn ) = ∫ f ( x) + ∑
k =1
α k hn 2 k
a

Hence, the Richardson extrapolation for raising the order of errors is applicable again:

T0 (hn ) ≡ T (hn ) (recursive trapezoidal rule)

4 k Tk −1 (hn ) − Tk −1 (hn −1 )
Tk (hn ) =
4k −1
T (h ) − T (h )
= Tk −1 (hn ) + k −1 n k k −1 n −1
4 −1
Numerical integration is then obtained by calculating the following table:

T0 ( h0 ) . . .
T0 ( h1 ) T1 ( h1 ) . .
T0 ( h2 ) T1 ( h2 ) T2 ( h2 ) .
T0 ( h3 ) T1 ( h3 ) T2 ( h3 ) T3 ( h3 )

The table can only be calculated row by row, except for the first column which directly results from the
recursive trapezoidal rule. The columns result from the extrapolation formula. The highest-order result is then
found as the last value in each row. Calculation of the rows is stopped as soon as the result of the row does not
change any more (referring to a predetermined precision).

Calculation of the integral according to the recursive trapezoidal rule including additional Richardson
extrapolation is called Romberg integration.

27
4.2.4 Gaussian integration
The trapezoidal rule includes the evaluation of the function f at points xi, which are equidistantly distributed in
the interval [a,b]. The estimated value of the integral then reads:

b N −1

∫ f ( x)dx ≈T (h) = 1 2 h( f 0 + f N ) + h ∑
i =1
fi , f i = f ( xi )
a

Hence, it represents a weighted sum of values of the function fi, the weights wi being 0.5 in this case for the end
points and 1 for the other points. The question arises of whether a nonuniform distribution of the values xi yields
a higher precision and if so, which weights are to be used for summing up. Here, just as in many other
mathematical problems, Gauss provides us with an answer. He succeeded in showing that a certain selection of
weights and sampling points xi yields an approximation of f by a polynomial of the order 2N, if the function is
evaluated at N points only. Thereby the order doubles as compared to the trapezoidal rule. The values of wi and
xi are found in tabular form for different N in mathematical tables of formulas. As mentioned above, a higher
order of polynomials does not necessarily mean a higher precision. The former is only true for functions f that
can be approximated by polynomials of high orders.

Extending the Gaussian integration allows for the integration of certain weighting functions into the integrands:

b N

∫ f ( x ) g( x )dx ≈∑
i=0
wi f i , f i = f ( xi )
a

The values of wi and xi are found in tabular form for different weighting functions g(x) and different N in
mathematical tables of formulas. The following weighting functions are commonly in use:

Integration formula Weighting function

Gauss-Legendre g( x ) = 1

( )
Gauss-Chebychev −1 2
g( x ) = 1 − x 2
Gauss-Hermite
g( x ) = e− x
2

Gauss-Laguerre g( x ) = x α e− x
Gauss-Jacobi
g( x ) = (1 − x ) (1 + x )
α β

The advantage of these rules is that they are precise for integrands which can be written as a product of a
polynomial and one of the weighting functions listed above.

28
4.3 Exercises5

1. Write a Matlab function deriv for numerically calculating the derivative of any arbitrary function f at a
predetermined point x according to the right-hand and centered formulas. Bisect the interval h proceeding
from h=1, until the value does not change any longer referring to a predetermined precision (e.g. 10-5).
(Hint: Write f as an independent Matlab function and call it to calculate the values of the function f within
deriv).
2. Calculate the derivation of the following function numerically using deriv:

a. f(x)=sin(x) , x = 5/4*π
b. f(x)=e-x/x2 , x = 0.5

Plot the absolute error between the numerically calculated and the exact derivatives as a function of the
interval h. Check whether the absolute error decreases with h (right-hand formula) and h2 (centered formula),
respectively (Hint: Use a double logarithmic representation).
3. (*) Improve the function deriv of Exercise 1. using Richardson extrapolation and calculate the derivatives
2.a and b anew. At which interval h is the precision equal to that yielded by the centered formula?

Use the Matlab function rombf for Romberg integration for the next exercises:

4. Determine π by calculating the integral over a quarter unit circle with a precision of 20 digits.

5. (*) In Fresnel’s diffraction the Fresnel integrals play an important part:

C( w) = ∫ cos(1 2π ⋅ x 2 )dx S( w) = ∫ sin(1 2π ⋅ x 2 )dx
w w

0 0
Write a function for calculating C(w) and S(w). Plot S against C for w = 0.0, 0.1, 0.2,..., 5.0. What does the
plot represent?

5
Exercises with (*) are optional.

29
5 Ordinary differential equations (ODE)
Differential equations (DE) are of special relevance in physics since Newton dropped the apple and described
this event (and thus gravitation) by means of the differential calculus, which was new at that time. Ordinary
differential equations (ODE) contain only derivatives for one quantity, e.g. the time t. Differential equations
containing derivatives for several different quantities are called partial differential equations (PDE). Spring
pendulum and planetary motion are examples of ODE. Schrödinger’s equation and the wave equation are
examples of PDE. Many DE cannot be solved analytically so that numerical methods are required again.

5.1 Mathematical formulation of the initial value problem

For deriving the general form of ODE we proceed from Newton’s law F = ma , which states a relation
between force and acceleration6. For a given force we can derive the equations of motion:

d2r F
a= 2 =
dt m

This vectorial equation of the 2nd order (2nd derivation occurs) can be rewritten into a system of two equations
of the 1st order by inserting a new intermediate variable. This is always feasible, but not plain in general, as there
are several possibilities to define the variable. In this case it is advisable to use velocity as an intermediate
variable:

dr
=v
dt
d v F ( r ( t ), v ( t ), t )
=
dt m

r ( t = t0 ) = r 0
v ( t = t0 ) = v 0

Now we have two equations of the first order. By integration of the two derivatives we obtain two integration
constants, which are determined by selecting the initial conditions. The force can depend on position and
velocity (and hence implicitly on time) and also explicitly on time. Let us assume that the problem is two-
dimensional in the Cartesian coordinates x and y. Then we can convert the problem into a general form by
introducing further variables:

 z1   rx   H1   vx 
       
 z   ry   H 2   vy 
z = 2≡  and H = ≡
z v H F (z ) m 
 3  x  3  x 
 z  v   H   F (z ) m 
 4  y  4  y 

The explicit dependence on time of all quantities has been omitted. According to this definition the general form
of an ODE as an initial value problem is then:

dz
dt
= H ( z ( t ), t ) ( )
with z t 0 = z 0

6
All vectorial quantities are underlined.

30
In the above equation, z is the required solution vector (in the simplest case it is a scalar function of t) and H, a
linear or nonlinear vectorial function in the arguments z and t. z0 is the initial value, t mostly represents the time.
The function H can be interpreted geometrically as a vector field, which predetermines a direction for z at each
point of the space spanned by z and t.
A given problem should be converted into the general form shown above. Different formulations may arise
depending on the intermediate variables chosen. Appropriate intermediate variables are not always easy to find,
as numerical problems may arise according to the individual choice which are not always assessable in the
beginning. In case of doubt several formulations have to be found and calculated.

5.2 Simple methods

Simple approaches directly resulting from the derived formulas for numerical differentiation are shown in the
following. The different formulations are very easy to calculate and they appear to differ relatively little.
However, experiments show that they are optimally applicable to different problems. In general, it is not clear
beforehand which method is suitable in which case, as the numerical errors are difficult to assess. Only intuition
and experience can help here. The idea behind all of these methods is to discretize the derivatives by appropriate
approximations and thus to calculate the values of the variables within a time step from the values of the
variables in the preceding time step.

5.2.1 Euler’s method

Let us proceed from the simple equations of motion again:

dr
=v
dt
d v F ( r ( t ), v ( t ) , t )
= = a( r ( t ), v ( t ), t )
dt m
Approximating the derivative according to the right-hand formula using a time step τ, we obtain:

r( t + τ ) − r( t )
+ O(τ ) = v ( t )
τ
v( t + τ ) − v( t )
+ O(τ ) = a(r ( t ), v ( t ), t )
τ
or

r( t + τ ) = r( t ) + τ v( t ) + O τ 2 ( )
( )
v ( t + τ ) = v ( t ) + τ a ( r ( t ), v ( t ), t ) + O τ 2

Note that the order of the truncation error increases after multiplication by τ. Introducing the following notation
simplifies the equations:

f n ≡ f ( nτ ) ; n = 0,1,2,...

Hence it follows:

( )
r n +1 = r n + τ v n + O τ 2

v n +1 = v n + τ a n + O(τ 2 )

For the one-dimensional case this Euler step can be interpreted geometrically (Figure 1). Proceeding from the
current state rn the gradient of the derivative field v is determined at point rn. With this gradient the new value
rn+1 is determined by linear continuation to the next sampling point.

31
 Euler-Verfahren
0.6

0.55 Geschwindigkeitsfeld
0.5

0.45

0.4
Ort r

0.35 v
n
wahre Lösung
0.3 r
n
0.25
r
n+1
0.2

0.15

0.1
2 3 4 5 6 7 8 9 10 11 12
nτ Zeit t (n+1)τ
Figure 1: Geometrical interpretation of Euler’s method (one-dimensional) on the basis
of a simple system of solutions (exponential functions).

In order to calculate the trajectory as a series of discrete sampling points nτ, we obtain the following algorithm:

1. Predefine the initial conditions r0 and v0.

2. Select a time stepτ.
3. Calculate the acceleration from the current values of r and v.
4. Apply Euler’s method for calculating the new values of r and v.
5. Return to step 3.

5.2.2 Euler-Cromer method

The following equations are a simple modification of Euler’s method (not derived here):

( )
v n +1 = v n + τ a n + O τ 2

r n +1 = r n + τ v n +1 + O(τ 2 )

The small alteration lies in the fact that the velocity newly calculated in the observed time step is used for
calculating the new position vector. Note that the equations can only be calculated in the stated succession for
this reason. This method has the same truncation error as Euler’s method, but it works much better in some
cases.

5.2.3 Midpoint method

The midpoint method is a mixture of Euler’s method and the Euler-Cromer method:

v n +1 = v n + τ a n + O τ 2 ( )
v n +1 + v n
r n +1 = r n + τ
2
+ O τ3 ( )

32
Here, the mean of the two velocity values are taken which raises the order of the position equation by 1.
Inserting the first into the second equation yields:

1
r n +1 = r n + v n τ + a τ2
2 n
This equation is exact if the acceleration a is constant. This makes the method interesting in case of problems
with a steady and slowly changing acceleration.

5.2.4 Verlet method

The Verlet method has the order 4 in the position equation and thus is especially suitable for problems in which
this variable is of particular relevance. Moreover, the position equation can also be solved alone, if the force or
acceleration depends on the position exclusively and the velocity is of no interest. For deriving the Verlet
method we calculate with the first and second derivatives of the position and approximate them according to the
centered formula for the first and second derivative, respectively:

dr r n +1 − r n −1
dt
=v
2τ
( )
+ O τ 2 = vn
⇒
d2r r n +1 + r n −1 − 2r n
dt 2
=a
τ 2 ( )
+ O τ 2 = an

A conversion yields:

r n +1 − r n −1
vn =
2τ
+O τ2 ( )
( )
r n +1 = 2r n − r n −1 + τ 2 a n + O τ 4

The disadvantage of this method is that it is not ”self-starting“, i.e., knowing the initial values r0 and v0 is not
sufficient to iterate the trajectory. On the contrary, r0 and r1 must be known, which generally is not the case. In
order to start the iteration, r1 can be calculated according to Euler’s method and then be further iterated using the
Verlet method.
Remark: If the acceleration does not depend on the velocity, the position equation can be iterated alone without
calculating the velocity.

5.3 Global and local truncation error

Depending on the method applied, a truncation error occurs at each time step, e.g. of order O(τ2) for Euler’s
method. This error is called local error in the following, as it occurs per time step. Considering a time interval
of T= Nτ for calculating the trajectory (hence N time steps), we maximally obtain the following global error in
the observed variable at the end of the time interval:

global error ∝ N × local error

( ) ( )
= NO τ n = (T τ )O τ n = TO τ n −1 ( )
For example, Euler’s method has a local truncation error of the order O(τ2), but a global truncation error of only
O(τ). This equation yields only an estimate of the global error, because it is not clear, whether the local errors
accumulate or cancel (i.e. interfer constructively or destructively). Obviously this depends on whether the vector
field of the derivatives in the observed range of variables has points of inflection and if so, how many of them
and also on the kind of approximation (thus on the method applied). This missing knowledge of the error
accumulation/prpagation implies that it is generally impossible to indicate, whether a certain method is suitable
for a certain problem.

5.4 Runge-Kutta method of the 4th order

To get solution formulas of a higher order we have to leave the empirical path of finding simple solutions. For
this purpose we proceed from the general form of the initial value problem:

33
 dz
dt
= H ( z ( t ), t ) ( )
with z t 0 = z 0

The simple approaches estimate the multi-dimensional derivative vector H from the current state z at time t once
and then calculate the state at the next sampling time t+τ. For the one-dimensional case

dz
= H ( z ( t ), t )
dt
this can be interpreted geometrically as shown above. It is imaginable to calculate several possible values of H in
the following way:

H1 = H ( z , t )
(
H 2 = H z + 12 τ H1 , t + 12 τ )
H 3 = H ( z + 12 τ H 2 , t + 12 τ )
H 4 = H ( z + τ H3 , t + τ )
The values Hi can only be calculated in the given succession. For this, an Euler step is performed twice at half
the step size. In the second step, the derivative assessed in the first half Euler step is used. Additionally, a full
Euler step is performed using the gradient assessed in the second half step. The values Hi can be interpreted
geometrically (as shown in Figure 1 for the simple Euler step) (cf. Fig. 2).

Fig. 2: Calculation of the estimates of the derivatives for the

Runge-Kutta method of the 4th order (cf. text).

The question arises of how the state at the next sampling time (t+τ) is calculated from these four assessments of
the derivative. It seems to be appropriate to take the weighted mean of the respective derivatives and then to
perform an Euler step with this mean gradient. By means of the Taylor expansion (cf. textbooks) we can prove
that the following weights are optimal, i.e., they yield a maximum order:

34
 4
z (t + τ ) = z (t ) + τ ∑ wi Hi + O(τ 5 )
i =1

wi = { 16 , 13 , 13 , 16}
Hence, the sum of weights is 1 as required. The gradients (H2 and H3) are weighted double. Altogether this
procedure can be described vectorially for the general form of the initial value problem as follows:

( )
z (t + τ ) = z (t ) + 16 τ ( H 1 + 2 H 2 + 2 H 3 + H 4 ) + O τ 5

H 1 = H ( z, t )

(
H 2 = H z + 12 τ H 1 , t + 12 τ )
H 3 = H ( z + 12 τ H 2 , t + 12 τ )
H 4 = H ( z + τ H 3, t + τ )

The Runge-Kutta method described here is of the 4th order in each variable (i.e. in each component of z) and
therefore has a local truncation error of the order τ5. Of course, formulas of even higher orders could be found in
the same way, but they are more complex, i.e. require more evaluations of the function H per time step. Runge-
Kutta of the 4th order has turned out to be the best compromise between precision and complexity. All textbooks
agree that it represents the „workhorse“ for solving the ODE, especially when it is combined with an adaptive
step size selection with a time step τ varying according to the assessed local truncation error (cf. next section).
Furthermore, the so-called predictor-corrector methods enable the local errors and the step sizes to be
controlled more precisely and also enable a correction term to be calculated for a higher precision. These
methods will not be described in this context. Please refer to relevant textbooks.

5.4.1 Adaptive step size selection

The step size should be selected such that it is a fraction of the system-inherent time scales (e.g. 1/50th of the
oscillation period of the system)7. Nevertheless, the system can reach states with strong differences in the
absolute value of the derivative field in the course of time development. In case of large derivatives (i.e. large
changes of the state of a system), a small time step has to be used in order to keep the local error below a
predetermined bound. Within ranges of slowly changing states of a system, larger step sizes are sufficient. The
number of states to be calculated can be minimized by adaptive step size selection. The additional time required
for calculating the step sizes is less than the saved calculation expenditure in most cases.
The object of the adaptive calculation of step sizes is to assess the local truncation error for each time step anew
and to select the step size such that it keeps below a predetermined bound. For this purpose the new state of a
system is calculated twice, once via one step with the current step size and once via two steps with half the step
size:

z( t ) → → → z1 ( t + τ )

( )
z ( t ) → z t + 12 τ → z 2 ( t + τ )

The local truncation error is assessed as the difference of both values:

∆ z ≈ z 2 − z1

Then the relative error can be assessed as follows:

7
Systems with strongly differing time scales are problematic. The shortest time scale always determines the step
size so that short step sizes have to be used, even if the long time scales predominate. These so-called stiff ODE
are often resolved with implicit schemes (cf. textbooks).

35
 ∆z z 2 − z1
r( z ) = ≈
z 1
2 ( z 2 + z1 )
Since the local truncation error is proportional to the fifth power of the step size τ (Runge-Kutta 4th order), the
current values of step size and truncation error can be related to the values to be expected upon selection of a
new step size τnew:

∆ z r (z )
5
 τ 
  = =
 τ new  ∆ r

This equation can be written equally for the relative errors, because the scaling factors are cancelled.With a
predetermined relative error bound r (e.g. 10-3) the new step size can thus be calculated, which keeps within the
predetermined error bound according to the relative error assessed before. Therefore, the equation is resolved for
τnew:

15
 r 
τ new = Sτ  
 r (z ) 

An additional factor of safety S≅0.9 is incorporated such that the change will not become too large. Moreover,
the new step size should not shift to higher or lower values by more than a predetermined factor (factor about 4
and ¼, respectively). Such a limitation is necessary, because the assessment and the quotient in the equation
above may become singular.
The step size is repeatedly reduced per time step according to the given formula (including limitation), until the
truncation error assessed with the new step size becomes smaller than the given truncation error. Vice versa, the
step size is increased once only per time step. Thus, we reduce the danger that the time step becomes too large
due to inaccurate assessment of the truncation error, which may extremely impair the precision. Assessing the
time step too small will lead to a higher expenditure of work, but will not reduce the precision.

5.5 Application of various methods to solve initial value problems

The presented methods can be applied to different simple physical systems. Some of them are described in the
following. The question of how these different methods behave when solving physical problems is investigated
(partly by means of exercises).

5.5.1 Projectile motion

 x
The trajectory of a ball is studied in two dimensions   . We assume a gravitational force in negative y-
 y
direction and a frictional force proportional to the square of the velocity and directed against the current
direction of motion. The equations of motion then read:

36
 dr
=v
dt
d v F (r (t ), v(t ), t )
=
dt m
F a (v ) + F g
=
m
with
1
F a (v ) = − c w ρ A v v
2
 0
F g = −mg  
 1

with:

cw ≅ 0.35 (drag coefficient)

kg
ρ ≅ 12
. 3 (density of air)
m
[ ]
A m2 and m[ kg ] as the cross-sectional area and mass of the ball
 kg ⋅ m 
g  2  gravitational acceleration.
 s 

If the frictional force is negligible, the analytical solutions are known. If the ball is thrown with an original
velocity v0 at an anagle Θ to the horizontal line, the throwing range, maximal height, and flying time are:

2v 02
x max = sin Θ cos Θ
g
v 02
y max = sin 2 Θ
2g
2v
T = 0 sin Θ
g

This limiting case can be used to verify the numerical solution. Moreover, a physically reasonable step size can
be selected as a fraction of the expected flying time. The projectile motion can be integrated using Euler’s
method (cf. Exercises).

5.5.2 Physical pendulum

We look at a pendulum, which is modeled as an ideal mass point of the mass m and oscillates around an axis by
means of a massless arm of the length L. The equation of motion for this system reads (see textbooks):

dΘ
=ω
dt
dω g
= − sin Θ
dt L
g being the gravitational acceleration. The angular acceleration is a nonlinear function of the elongation and
there is no analytical solution. For small excursion angles Θ, however, the equation can be linearized:

37
 dω g
sin Θ ≈ Θ ⇒ =− Θ
dt L
For this mathematical pendulum we obtain the solution:

 t  L
Θ( t ) = exp i2π  , T0 = 2π
 T0  g

If the ratio L/g is normalized to 1, the oscillation period T0 equals 2π. The step size can then be adjusted as a
physically useful fraction of the oscillation period. Using the law of conservation of energy it can be checked,
whether the numerical solution is physically correct. For this purpose the total energy

1 2 2
E ges = mL ω − mgL cos Θ
2
is numerically calculated for each time step and plotted as a function of time (see Exercises).

Nonlinear effects are found at large elongation angles as a deflection of the oscillation form in relation to the
sine function. In this case the oscillation period also depends on the maximum excursion angle θm (for derivation
see next subsection):

L   1 
T0 = 4 ⋅ K  sin  θ m  
g   2 
L 1 2 
≈ 2π 1 + θ m + … 
g  16 

For smaller θm the formula turns into the above-mentioned formula for the linear case. K is the complete elliptic
integral of the 1st category:

π 2
1
K (x ) = ∫ dz
0 1 − x 2 sin 2 z

If a frictional force proportional to the angular velocity as well as a harmonic excitation are introduced in
addition to the nonlinear term of the restoring force, a path into deterministic chaos can be demonstrated by this
simple system:

dΘ
=ω
dt
dω g
= − sin Θ − αω + A0 sin(ω0 t )
dt L
If the excitation amplitude is sufficient for the pendulum to overturn, this leads to a chaotic motion depending on
damping and excitation. For further analysis, the trajectory of motion can be observed in the phase space, i.e.
the space spanned by Θ and ω. If only points of the trajectory are plotted, which belong to a certain phase of
excitation (e.g. zero phase), this leads to the so-called Poincaré section. If the motion is periodic there are only
a few entries in the Poincaré section, the number of points indicating the multiplicity of the oscillation period
relative to the excitation period. The path into chaos can then be understood with the help of the so-called period
duplication (bifurcation): Plotting the Poincaré points (only Θ components) as a histogram as a function of the
excitation amplitude, we obtain windows with periodic motions. With increasing excitation amplitude, the
multiplicity of the period doubles. A multitude of such period duplications then leads to chaotic motions. This is
an example of deterministic chaos: The underlying equation is well-defined and deterministic. Nevertheless, the
motion is not predictable, because it extremely depends on the initial conditions. Additionally, the truncation and
roundoff errors accumulate when pursuing the numerical calculation, leading to unpredictable global errors.

38
5.5.2.1 Derivation of the formula for the oscillation period
The undamped physical pendulum oscillates periodically. The maximum excursion be θm. At the stationary
point, i.e. at θ = θm , the angular velocity ω equals 0. The total energy at this point then is
E ges = −mgL cos Θm .
Θ= Θm
Due to the constancy of the total energy we obtain for all time points

1 2 2
mL ω − mgL cos Θ = − mgL cos Θm
2
2g
⇒ ω2 = (cos Θ − cos Θm )
L
dΘ 2g
⇒ = (cos Θ − cos Θm )
dt L
L dΘ
dt =
2 g cos Θ − cos Θm

In the time period from t = 0 to t = T0/4 (1/4 of the required oscillation period) the pendulum oscillates from θ =
0 to θ = θm. Thus, both sides of the equation can be integrated:

T0 4 Θm
L dΘ
∫0 dt = 2 g ∫
0 cos Θ − cos Θ m
Θm
L dΘ
⇒ T0 4 = 2
g ∫
0 sin (Θ m 2 ) − sin 2 (Θ 2 )
2

with the replacement

cos 2Θ = 1 − 2 sin 2 Θ

An evaluation of the expression by means of the above-mentioned complete elliptic integral of the 1st category
yields the above-mentioned formula for the oscillation period T0.

5.5.3 Planetary motion

We observe the motion of a body in a central potential, e.g. a comet in the gravitation field of the sun.
Neglecting all further forces (such as planets, solar wind etc.) the force affecting the body in the coordinate
system of the sun is:

GmM
F (r ) = − 3
r
r

m being the mass of the body, M the mass of the sun, and G the gravitation constant. It is useful to calculate in
astronomical units. Then, GM=4π2 AU3/yr2 with AU=1.496x1011m being the average distance between sun and
earth, and yr being one year. Furthermore, the mass m of the body can be set to 1. The equation of motion can
now be integrated numerically. From Kepler’s laws we know that the solutions are ellipses. Furthermore, total
energy and angular momentum

1 GmM
L = r × ( mv )
2
E= mv − ,
2 r

are conserved. These physical boundary conditions allow for a verification of the numerical solution.
Conservation of energy and angular momentum are also applicable to the verification of numerical solutions of
many-body problems, which cannot be solved analytically anymore.
Applying the different methods, we find that Euler’s method fails for (approximately circular) orbits, while the
Euler-Cromer method is particularly suitable. In case of eccentric motion the Euler-Cromer method fails, too. In

39
this case, the adaptive Runge-Kutta method of the 4th order is the favourable method, which also applies to
three- or many-body problems.

5.5.4 Lorenz model

The weather was previously assumed to be unpredictable only because of the great number of variables.
Consequently, the development of efficient computers nourished hopes that they would enable long-term
forecasts. This hope had to be abandoned. We know today that the weather is intrinsically unpredictable,
because the underlying equations are nonlinear and their solution reacts extremely sensitively to changes of the
initial conditions (“butterfly effect”). This can already be demonstrated by means of systems with few variables.
For example, E. Lorenz investigated a model with three variables which describes the convection of a liquid
between two parallel plates with a constant difference in temperature:

dx
= σ ( y − x)
dt
dy
= rx − y − xz
dt
dz
= xy − bz
dt
The derivation of the model is found in textbooks (e.g. H.G. Schuster: ”Deterministic Chaos“). Please note that x
represents approximately the circular flow rate, y the temperature gradient between rising and dropping liquid
volumes, and z the deviation of the course of temperature from the equilibrium (thermal conduction only). The
parameters σ and b depend on the properties of the liquid as well as on the geometry of the plates. Typical
values are σ = 10 and b = 8/3. The parameter r is proportional to the imposed temperature gradient.
In order to investigate the sensitivity of the system to the initial conditions, it is the obvious thing to integrate the
equations with two slightly different initial conditions (Runge-Kutta with adaptive step size) and to plot the time
courses of the variables for comparison (cf. Exercises). Moreover, it is an obvious choice again to represent the
trajectory in the phase space. It moves on an aperiodic orbit which is relatively well localized in the phase space,
i.e. on an attractor (cf. Exercises).

5.6 Boundary value problems

So far, only the initial value problem has been treated: From a completely known initial state the temporal
development is iterated. Of course, as many initial values as there are integration constants have to be specified.
However, they do not need to be determined at a fixed time t0. This leads to the so-called boundary value
problems, in which part of the initial state at t0 and part of the final state after a time T are given. For example,
the locations at times t=0 and t=T may be predetermined for a projectile motion, but the velocity may be
unknown. Then the above-mentioned formulas for calculating the temporal development cannot be directly
applied any longer. Instead, assumptions about the complete initial state would have to be made, it would have
to be iterated according to the above-mentioned formulas and varied until all boundary conditions were fulfilled.
Of course, this means that the trajectory has to be calculated many times, until the right one is found. Relaxation
methods are better suited to resolve boundary value problems. They are explained when treating partial
differential equations (see below).

40
5.7 Exercises

1. Write a function ball, that numerically calculates the trajectory of a volleyball in two dimensions. Apply
Euler’s method to this problem. The initial conditions of the problem x ( 0) and v ( 0) should be stated as
parameters. The trajectory should be represented graphically and the time of flight and throwing range
should be calculated. Hint: A moving ball is decelerated by air friction as follows:
A
a friction = −0.5c w ρ v v with
m
cw ≅ 0.35 (drag coefficient)
kg
ρ ≅ 12
. 3 (density of air)
m
[ ]
A m2 and m[ kg ] as the cross-sectional area and the mass of the ball.

Gravitation should not be neglected, either.

2. Calculate the oscillatory behaviour of a physical pendulum using Euler’s and Verlet’s methods. The function
pendulum should expect the initial conditions as well as step size and kind of solution (Euler / Verlet). The
elongation angle and the total energy in dependence on time are to be plotted.
Recall the equations of motion:
dω g
= − ⋅ sin( Θ)
dt L
dΘ
=ω
dt
The total energy of the system is:
1
E= m L2ω 2 − m g L cos( Θ)
2
Hint:
g E
Set = 1 to simplify matters and observe the standardized energy .
L m ⋅ L2
[ °
] [ °
]
Examine the cases Θ 0 = 10 ; τ = 0,1 and Θ 0 = 10 ; τ = 0,05 using Euler’s method and

[Θ 0 ] [ ]
= 10 ° ; τ = 0,1 and Θ 0 = 170° ; τ = 0,1 using Verlet’s method!
What can we say about the suitability of the methods for this problem ?
3. Observe the harmonically driven, damped physical pendulum:

dΘ
=ω
dt
dω g
= − sin Θ − αω + A0 sin(ω0 t )
dt L
g
a.) Integrate the system using Verlet’s method. Set ω0 = ≡ 1 and calculate about 80 periods
L
(t=0 - 80*2π). Write a plot program drawing the trajectory in the phase space (ω over Θ) as well as the
elongation Θ and the total energy as a function of time.
b.) Generate the following oscillatory behaviour by an appropriate selection of parameters:
(i) Oscillation of a damped pendulum without excitation (Hint: Start with
A0 = 0 and vary the damping α until the pendulum dies off after about 5-10
periods from a starting position of Θ =90 Grad).
(ii) Limit cycle (Hint: low damping and low excitation, oscillation without
overturn).

41
 (iii) Beating (Hint: no (!) damping and low excitation, oscillation
without overturn).
(iv) Chaotic motion (Hint: Increase excitation up to and beyond an
overturn).
Explain the different types of oscillations and document the sets of parameters found.
c.) Proceeding from b.), generate a Poincaré section by plotting only those points in the phase
space that belong to the zero phase of excitation:
(
A = A0 sin ω0 t ∗ ) , ω0 t ∗ = n2π , n ∈N .
Reject those points that are taken during the build-up time (about 20 periods, i.e. until
t=40*π). What can we tell about the four cases of b.)?
d.) Plot the Poincaré points (Θ-components) as a histogram as a function of the excitation amplitude
(“bifurcation diagram”). Generate a bifurcation diagram with excitation amplitudes
between approx. 0 and 5 rad/s with a step size of 0.5 rad/s. Use the value of damping
from exercise 1.b.i). Identify the interesting range of transition from periodic to chaotic
motion and calculate this area with a solution of von 0.1rad/s anew. Describe the result.
4. Observe the body in the central potential. It is affected by a gravitational force:

GmM
F (r ) = − r
r3

a.) Program the adaptive Runge-Kutta method of the 4th order. For this write a function
[x, t, tau] = rka(x,t,tau,err,derivsRK,param) for calculating a time step, with:
x = current value of the dependent variable
t = independent variable (usually time)
tau = current step size
err = desired local truncation error (appr. 10^-3)
derivsRK = right side of DE (name of function returning dx/dt
(calling format: derivsRK(t,x,param)).
param = extra parameter of derivsRK
x = new value of dependent variable
t = new value of independent variable
tau = recommended value for tau for next call of rka
b.) Integrate the system with Euler-Cromer and rka for a circular and an eccentric path. Plot
each trajectory in the position space for 10 periods as well as the course of the total energy.
How do the two methods differ?
c.) Draw a log-log plot of the step size determined for each step by rka as a function of the related radial
distance from the origin for the eccentric path. Which law can be derived? (Hint: Kepler’s laws).
5. Observe the Lorenz model

dx
= σ ( y − x)
dt
dy
= rx − y − xz
dt
dz
= xy − bz
dt
Integrate the system with rka. Plot x and z as a function of t as well as the trajectory in the phase space (zx-
and yx-projection). Use the parameters σ =10, b=8/3 and r=28 as well as the initial conditions:
a.) x=1, y=1, z=20
b.) x=1, y=1, z=20.01

How do the paths belonging to the different initial conditions differ?

42
6 Solving systems of equations
The preceding section dealt with the solution of ordinary DE

dz
dt
= H ( z ( t ), t ) ( )
with the initial state z t 0 = z 0

Starting from the initial state, the temporal development of the system can be determined as a series of states
z(nτ ) , e.g. by means of the Runge-Kutta method. So far, only autonomous systems have been treated in
( ) ( )
which the force and thus H does not explicitly depend on time, i.e. H z( t ), t = H z( t ) . For this class of

systems there often exists a class of initial conditions z 0 for which z( nτ ) ≡ z 0 . These states in the N-
∗ ∗

dimensional phase space (state space) are called stationary. The following equivalence is easy to realize:

∗
z 0 stationary state ⇔ ( )
H z ≡0 ,
∗

because the latter condition means

dz
dt
( )
≡ 0 . The vector equation H z ∗ = 0 contains N equations with N
unknown quantities. Therefore, the roots of this system of equations have to be found in order to find the
stationary states. This problem is divided into two classes. If H is a linear function in z or can be linearized in the
surrounding of expected stationary states, methods of linear algebra can be used to solve the problem. If H is
nonlinear (e.g. the physical pendulum comprises the nonlinear cosine-function), Newton’s gradient method can
be applied. Both classes of equations are explained in the following. Additionally, there are special methods in
case H is a polynomial in z, which will not be treated in this context (please refer to textbooks and to the
function roots in Matlab, which finds zeros of polynomials).

6.1 Linear systems of equations

A linear system of equations generally consists of M equations with N unknown quantities:

a11 x1 + a12 x2 + ... + a1N x N − b1 = 0

a M 1 x1 + a M 2 x2 + ... + a MN x N − bM = 0

xi being the unknown quantities to be determined and aij and bi the (fixed) coefficients. This system can be
written in a matrix form8:

A x − b = 0 or A x = b ,

with the definitions:

 x1   b1 
 a11 a1N     
   x2  b 
A= , x =  und b= 2 
a a MN     
 M1 x  b 
 M   M 

A being the (MxN) matrix of the coefficients, x, the (unknown) solution vector, and b, the vector of the
coefficients of the zeroth order. Example:

8
Matrices are marked by a double underscore.
43
  x1 
 
2 x1 + x2 − 4 = 0  x2 
or
4 x1 − 2 x2 − 2 = 0 2 1   4
  =  
 4 − 2  2

Regarding the dimensions of the problem (MxN) three cases are distinguished:

1. M>N: The system is overdetermined, i.e. there are more equations than unknown quantities. In this case –
if the equations are linearly independent – only an approximate solution is possible (see next chapter).
2. M<N: The system is undetermined, i.e. there are less equations than unknown quantities. The solutions x
lie within a solution space of the dimension N -M.
3. M=N: The system has a well-defined solution, if all equations are linearly dependent (determinant of A
being not equal to zero). If equations are linearly dependent, case 2 applies.

In the following all three cases are treated. For case 3 the elementary Gaussian and Gauss-Jordan methods are
applied, while the method of singular value decomposition is presented for cases 1 and 2. The latter is widely
applied in linear algebra.

6.1.1 Elementary operations

Basic algorithms for solving systems of equations (Gaussian and Gauss-Jordan methods, respectively), use
operations which do not alter the result (i.e. the solution vector x). These operations are called elementary
operations. The following elementary operations are feasible:

Operations Meaning
1 Exchange of rows in A and b Exchange the sequence of two equations
2 Form linear combinations of rows Form linear combinations of equations
in A and b
3 Exchange columns in A Exchange variables
(see remark)

Remark: Operation 3 means that the sequence of components in the solution vector x is permuted. The
permutations have to be made undone in order to obtain the original solution (exchanging, for example, the
components position and velocity would be a fatal error). If several operations of type 3 are performed, they
must be made undone in the reverse order. This requires a corresponding “bookkeeping“, i.e. memorizing the
operations of type 3 performed.

Example:

2 x1 + x2 = 4 4 x1 − 2 x2 = 2
⇔ 1st operation
4 x1 − 2 x2 = 2 2 x1 + x2 = 4
2 x1 + x2 = 4
⇔ 2nd operation (2nd row equals 1st row plus 2nd row)
6 x1 − x2 = 6
x1 + 2 x2 = 4
⇔ 3rd operation (exchange x1 , x2 in the solution vector)
− 2 x1 + 4 x2 = 2

All systems of equations on the right-hand side have the same solution (considering the exchange of the solution
components in case 3).

6.1.2 Gaussian elimination with backsubstitution

This method is suitable for case 1 (M=N). The underlying idea is to first convert the system by means of
elementary operations such that the matrix A takes on an upper triangular shape:

44
 ′
 a11 a1′N   x1   b1′ 
     
 0   x 2   b2′ 
 ⋅  =  
     
 0 0 a ′NN   x N   bN′ 


Remark: The components of A and b are written here with apostrophe, since the numerical values have changed
due to the elementary operations used for finding the triangular shape.

When the system is in the stated triangular form, the components of the solution can be iteratively determined by
means of backsubstitution:

x N = bN′ a′NN , to be calculated directly from the last row

1
x N −1 = (bN′ −1 − a′N −1,N x N ) , to be calculated from the second last row
a′N −1N −1
1  N 
xi =  bi′ − ∑ aij′ x j  , i = N − 1,N − 2,… ,1 , general form
aii′  j =i +1 

All components of the solution can be calculated from the first equation and by iteration of the last equation.
Proceeding from xN , xN-1 is calculated etc.

6.1.2.1 Generating the triangular shape and pivoting

The triangular shape required for the backsubstitution can be generated as follows:

1. Start out from the diagonal element a11.

2. Subtract the a1j / a11-fold of the first row from each row j below the diagonal element. This results in zeros
in the 1st column below the diagonal element.
3. Proceed to the next column and continue at with 1, however, starting out from the diagonal element of the
actual column.

The problem inherent to this method is caused by the numerical errors occurring upon division by the respective
diagonal element. The element by which we divide is called the pivot and the (empirically founded) method for
minimizing the error is called pivoting. Pivoting means that the rows below and the columns to the right of the
currently observed diagonal element (including the row and column with the current pivot) are exchanged such
that the biggest element becomes the pivot (elementary operations 1 and 3). Partial pivoting means that only
rows are exchanged. Here we save the necessary “bookkeeping“ when exchanging columns. Partial pivoting is
empirically found to be sufficient in general.

Remark: If the pivot (following pivoting) reaches the range of machine precision, the procedure is numerically
instable. If the pivot equals zero, this may indicate that the system of equations is really singular (determinant =
0), or that the value 0 has accidentally resulted from numerical inaccuracies. Hence, the problem cannot be
clearly analysed on the basis of the quantity of the pivot. In case of a small pivot we can only state that the
solution is instable or that the system may be even singular. The solution should be verified by evaluating the
equation using the estimated solution in that case.

6.1.2.2 Iterative improvement of the solution

The described Gaussian method comprises many calculation steps so that calculation errors may accumulate.
Therefore, the numerically found solution will generally deviate from the true solution by an amount larger than
the machine precision. In the following a method is described that allows for an iterative improvement of the
solution up to the order of magnitude of the machine precision. The approximate solution for the system of
equations

45
Ax = b

found using the Gaussian method be x‘ and x be the (unknown) true solution. By insertion the solution can be
verified:

′ ′
Ax = b

We define the numerical errors as:

δ x = x′ − x
δ b = b′ − b

δb is known (and is, as mentioned above, generally larger than the machine precision), while δx is unknown,
because the true solution x is unknown. δx can, however, be calculated proceeding from the approximate
solution and utilizing linearity:

′ ′
Ax = b
⇔ A( x + δ x ) = (b + δ b )
⇔ A x + Aδ x = b + δ b
⇔ Aδ x = δ b

The linearity of the system is utilized in the third equation. Because x is assumed to be the true solution, the first
addends on both sides of equation 3 can be dropped. Resolving the system of equations written in the fourth row
(applying the Gaussian method again), we obtain an approximate solution δx‘ for the true error δx. The
improved estimate of the solution is then:

″ ′ ′
x = x −δ x

This method can be iterated until the error δb reaches the order of magnitude of the machine precision. For this
purpose the improved solution x‘‘ is again inserted into the system of equations, a new δb is calculated, a new δx
is determined etc. The drawback of the described method is that the system of equations has to be resolved
several times. This expenditure of time is only justified in case a particularly precise solution is required.

6.1.3 Gauss-Jordan method

Gaussian elimination with backsubstitution is the most favourable method for calculating the solution in case 1
(M=N) regarding the expenditure of time. However, if the inverse matrix A-1 is to be calculated in addition, an
extended procedure known as the Gauss-Jordan method is more appropriate. For this purpose the following
matrix equation is defined:

[ ] [ ]
A x ∪Y = b ∪1
with:

 b1 1 0 0
   x1 y11 y1N 
0  
[b∪1] = 
0
and [x ∪Y ] =  
  x y yNN 
b 0 
0 1  N N1
 N
46
We can explain this formulation as follows: N+1 systems of equations are established simultaneously, each
column in [b∪1] and [x∪Y] representing an individual, independent system of equations. These systems can be
resolved simultaneously by converting matrix A into a unit matrix using elementary operations (Remark: Each
elementary operation then influences each column of [b∪1]) separately, as described above for the case of one
column. The right side of the equation then represents the unknown x which solves

Ax = b

and the unknown Y which solves

AY = 1 .

Thus, Y is the inverse of A.

6.1.3.1 Generating the unit matrix

The unit matrix can be generated with the following scheme which strongly resembles the Gaussian method:

1. Go through matrix A column-wise, starting out from column 1. i be the current column.
2. Take the diagonal element aii as pivot and perform pivoting.
3. Divide row i by the pivot. Hence, the diagonal element assumes the value 1.
4. Set the values of all other rows in column i to zero by subtracting the aji–fold of the i-th row from the j-th
row (for all rows j except the i-th row itself).
5. Proceed to the next column i+1 and continue with 2.

The notes from 6.1.1.1 on the numerical stability hold also for this extended method.

6.1.4 Singular value decomposition

As indicated above, the Gaussian and Gauss-Jordan methods are not suitable for systems of equations which are
almost singular because of the repeated division by the respective pivot. Moreover, it is not suitable for
overdetermined and underdetermined systems. In the following a method is presented that is numerically very
stable and is suitable for all classes of systems of equations. Furthermore, it enables a clear diagnosis of the
system’s numerical stability. In the literature this method is mostly called SVD (singular value decomposition).
SVD is based on a theorem of linear algebra the proof of which is not given here (refer to textbooks of linear
algebra). It says that each MxN matrix A with M≥N can be decomposed as follows:

w1 
 
A = U ⋅   ⋅V = U ⋅ diag (w j ) ⋅V
T T

 
 w N 
with:
U : MxN -matrix
V : NxN -matrix
wj : singular values

U and V being column-wise orthonormal, i.e.:

47
 U T ⋅U = V T ⋅V = 1
M 1≤k ≤N
⇔ ∑U ikU in = δ kn
i =1 1≤n ≤N
N 1≤k ≤N
∑VjkVjn = δ kn
j =1 1≤n ≤N

We will not deal with the details of how the decomposition is calculated. It is important to know that stable and
efficient algorithms are available for this purpose. In the following it is shown, how SVD can be used to
calculate the solution of linear systems of equations in the different cases.

Remark: The command SVD is available in Matlab. Attention: There are several variants of singular value
decomposition. The variant presented here is calculated by the SVD command, if 0 is given as additional
parameter: [U,S,V]=svd(A,0)).

Remark: SVD is also feasible in the case of M<N (underdetermined system), however, N-M of the singular
values are then equal to zero and the column-wise orthonormality of U and V does not apply to those columns
containing the zeros in the matrix of the singular values.

6.1.4.1 SVD of a quadratic matrix (case M=N)

In the case of M=N all matrices involved are quadratic. Then, row-wise orthonormality also applies to U and V,
i.e.:

T T
U ⋅ U = V ⋅ V = 1 is valid.

U and V thus are orthogonal and

T −1
U =U
T −1
V =V

Then, the inverse of matrix A is yielded as

1 w1 
  T
 ⋅ U = V ⋅ diag (1 w j ) ⋅ U
−1
A = V ⋅ T

 1 wN 


and the solution of the system of equations is:

A x = b ⇒ x = A b = V ⋅ diag (1 w j ) ⋅ U ⋅ b .
−1 T

This method will only go wrong, if one or several of the singular values become zero or lie within the order of
magnitude of machine precision. The more singular values get small the more instable the solution. In the first
place SVD enables a clear diagnosis of the situation owing to the analysis of singular values. Two cases are to
be distinguished (εm: machine prescision):

 
min  w 
j ∈1,N  j 
(i) The matrix is ill-conditioned, i.e. < εm
 
max  w 
j ∈1,N 
 j
(ii) The matrix is singular, i.e. ∃j : w j = 0 , j ∈ 1, N
48
6.1.4.1.1 Case (i): „ill-conditioned“
Formally, there is a unique solution in case (i). However, the Gaussian and Gauss-Jordan methods are
numerically instable in this case, whereas SVD yields a generally more robust and more precise solution in the
following way (without proof):

′
x = V ⋅ diag 1 w j  ⋅ U ⋅ b , with
T

 

′  0 wj < εm
1 wj = 
1 w j sonst

i.e., we simply eliminate those singular values the value of which is not exactly known. This is illustrated by the
fact that we throw out those equations and their combinations which more or less contain round-off errors only
by setting the 1/wj to zero. The solution x found in this way is generally better in the least-squares sense than the
solutions yielded by other methods (Gauss elimination, Gauss-Jordan), i.e.

2
r 2 = Ax − b

becomes smaller. The proof is not given here.

6.1.4.1.2 Case (ii): singular

The case (ii) is mathematically more complex. The inverse matrix A-1 does not exist in this case, however, the
system of equations

Ax = b

can have a solution space as solution. The equation

′
x = V ⋅ diag 1 w j  ⋅U ⋅ b , with
T

 
′  0 w j = 0
1 wj = 
1 w j otherwise

yields an “as optimal as possible“ special solution (without proof). “As optimal as possible” means that

2
r 2 = Ax − b

becomes minimal and that x has the minimum length of all vectors within the solution space. The solution space
L is yielded as the sum of the so-called nullspace and the special solution:

{′ ′
L= x+x :x ∈N }
The nullspace N is defined as follows:

N = {x : A x = 0}

Thus, the nullspace contains all vectors mapped on zero.

Remark: Those columns of V containing zeros in the related matrix of singular values (diag(wj)) span the
nullspace (without proof).

Remark: Those columns of U containing zeros in the related matrix of singular values (diag(wj)) span the range
(without proof). The range is the set of all vectors b, for which the system of equations Ax=b has at least one

49
exact solution. If b lies within the range of A, the above-mentioned special solution is exact. If b lies outside the
range, the least-squares condition applies.

Remark: The case (ii) mathematically corresponds to the underdetermined system (M<N), which therefore can
be resolved in the same way.

6.1.4.2 SVD for overdetermined systems (case M>N)

In case there are more equations than unknown quantities there is generally no exact solution to the system of
equations. SVD, however, finds the “optimal” approximate solution again (without proof):

x = V ⋅ diag (1 w j )⋅ U ⋅ b
T

2
with : r 2 = A x − b being minimal

Remark: Again, singular values smaller than the machine precision may exist here. Replace them by zeros
again in the matrix of inverse singular values.

6.2 Nonlinear systems of equations

When a system of equations cannot be linearized, the algebraic methods of linear algebra described above
cannot be applied. There is no choice but to apply iterative methods, of which Newton’s gradient method is
described below as a simple example.

6.2.1 Newton’s method

The gradient method searches the zeros of a system of equations by taking steps in the direction of the gradient
starting from the current estimate of the zero. Let us assume that z* is the solution of the system of equations

H ( z ) = 0 with H = (h1 (z ) … hN ( z )) and z = ( z1 … z N )

(Remark: H and z are row vectors here). Let us further assume that z1 is the estimate of the solution, which
deviates from the true solution by δz, i.e.

δ z = z1 − z ∗ .

With this definition and assuming z* to be the true zero, a Taylor expansion of the 1st order yields:

( )
0 ≡ H z = H (z 1 − δ z ) = H (z 1 ) − δ z ⋅ D(z 1 ) + O δ z
∗
( )
2

with the Jacobi matrix D:

∂h j ( z )
D( z ) = {Dij } , Dij =
∂z i

Neglecting the terms of higher order we thus obtain:

( )
0 = H ( z 1 ) − δ z ⋅ D( z 1 ) + O δ z
2

⇒ H ( z 1 ) ≈ δ z ⋅ D( z 1 )
⇒ δ z ≈ H ( z 1 ) ⋅ D −1 ( z 1 )
z 2 = z 1 − δ z = z 1 − H (z 1 ) ⋅ D (z 1 )
−1
⇒

z2 is a new estimate of z*, which can be iteratively inserted into the last equation again.

50
6.3 Exercises

1 2 3   1
   
1. The matrix A =  7 3 4  and the vector b =  2 be given. Solve the linear system of
   
 8 9 15  4
−1
equations Ax = b and calculate the inverse matrix A . Use the supplied Matlab script gaussj.m, which
performs the Gauss-Jordan elimination procedure. Try to understand the program and compare the result
with the installed Matlab functions inv (calculation of the inverse) and \ (calculation of the solution of the
linear system of equations; type ”help slash“ for hints on use).

2. The matrices
 4 −6 −2 5 
 1 − 5 3 10 − 8 − 9  
     −8 7 5 7 
A =  2 − 10 0 , B =  6 10 2  and C = 
    −1 4 0 − 10
2 8 8 1 7 3  
− 6 −1 − 4 0 
be given. Find the inverse matrices using the script gaussj.m. Why does this go wrong in all cases? In which
cases should a solution exist (calculate the respective determinants) and how could we still obtain a result
with this script in these cases? (Hint: ”pivoting“)

3. The Lorenz model

dx
= σ ( y − x)
dt
dy
= rx − y − xz
dt
dz
= xy − bz
dt
be given.
a) Show that it has the following stationary states:
x ∗ = y ∗ = z ∗ = 0 and x ∗ = y ∗ = ± b(r − 1), z ∗ = r − 1, ifσ ≠ 0
b) Search the stationary states for r=28, σ =10 and b=8/3 using Newton’s method. Perform the calculation
for several different initial estimates of the stationary states. How does the solution depend on it?
Hint: Use the supplied routine newtn, which implements Newton’s method. Then write a
routine [H,D]=fnewt(z,p), which calculates the function H as well as the Jacobi matrix D
from the current state z=(x,y,z) and the set of parameters p=(r,σ,b).

51
7 Modeling of data
Modeling of data is of great relevance when validating physical theories and models. In most cases some
parameters of the physical models are to be fitted to measured data and it is to be investigated to what degree the
fitted model describes the data. This procedure is illustrated by means of simple examples in the following. In
particular the techniques of fitting parameters and of assessing the goodness of fit will be described.

7.1 General mathematical formulation of the fitting problem

Let us assume N samples of data pairs (xi,yi) and a model to be fitted which predicts a functional relationship
between the dependent quantity y and the independent quantity x:

 a1 
 
y = f ( x, a ) with a =   ,
a 
 M

a being a vector consisting of M model parameters9. Of course, there should be clearly more data pairs than
model parameters, i.e. N >> M.

The activity of a radioactive isotope be given as an example:

a 
y = f ( x, a ) = a1 ⋅ e − a2 x with a =  1 
 a2 

y corresponding to the activity, x is the time and the parameters a1 and a2 describe the scaling of activity and
decay rate. The model describes an exponential decrease in activity.

To simplify the expressions used later on we also write the data in vectors:

 x1   y1 
   
x=  , y= 
x  y 
 N  N

The general procedure for fitting the model is as follows:

1. Define a ”merit function“ K (or: “cost function”) describing the deviation of the model from the data. The
larger the deviation of the model from the data, the larger the merit function. This monotonous relationship
applies to the least-squares merit function for example, which is defined as the sum over all data pairs of the
squared deviations of the values of dependent variables assessed by the model from the actual data, i.e.:
2

K (x, y , a ) = ∑ yi − f ( xi , a ) .
N

i =1

2. Minimize the merit function for a given data set by varying the parameter a.

The result of the fitting should then be:

1. Statement of the best set of parameters aopt (best fit) corresponding to the minimum of K.
2. Estimation of the error of parameters, i.e. statement of a range in which the ‘true’ parameters are found with
a given probability.

9
The problem is formulated here without limitation for one independent and one dependent variable. In case of
several independent variables, all of them enter the model function as arguments. Likewise, the model may
predict several dependent variables so that the model function would be a vector-valued function.

52
3. Statistical statement about the quality of the model (”goodness of fit“), i.e. clarification of the question of
whether the model is suited to describe the data sufficiently.

The last two points are often omitted, maybe, because they are time-consuming and difficult to understand.
However, they belong to a systematic model analysis. In the following methods of model fitting are described.
In each case the questions of which is the best merit function, how the merit function can be minimized, and
how the fit is to be assessed will be clarified.

7.2 Establishing the merit function: Maximum-Likelihood Methods

The Maximum-Likelihood concept (ML method) allows to derive well-defined merit functions K, which are
adapted to the respective fitting problem. The idea is to calculate the probability W that the observed data are
generated by the model starting out from a statistical model of the measuring process. The negative logarithm of
this probability is then used as the merit function K:

K (x , y, a ) = − log (W )
W : probability that the observed data are generated by the model
The set of parameters that minimizes K is stated as the best fit aopt:

a opt = argmin(K (x,y,a ))  = argmax (W (x,y,a ))

 
a  a 

The minimization of K corresponds to the maximization of the probability W due to the monotony of the –log(x)
function. Why is this procedure not called “Method of maximal probability” instead of “Maximum-Likelihood
method“? The answer to this question is that the statement of the probability W refers to data and not to the
model parameters: If the model is correct, there is only one “true” set of parameters and to state the probability
of the set of parameters makes no sense. Therefore, the expression likelihood of a set of parameters is defined as
follows:

The likelihood of a set of parameters a corresponds to the probability W that the

observed data are generated by the model with this set of parameters

In the following, examples are given of how W and K, respectively, can be calculated from assumptions on the
statistics of the measuring process. In particular, the statistical prerequisites of the measuring process are
explained which render the least-squares merit function optimal in the sense of the ML concept. Starting from
that fact, extensions of the least-squares method are stated based on slightly changed statistical assumptions on
the measuring process (Chi-square method and robust fit methods).

Remark: The ML method does not detect systematic errors of the measurement process, but only models
random errors. Systematic errors lead to a misfit of parameters.

7.2.1 Least-Squares method

The Least-Squares method is derived on the basis of the following assumptions:

1. The measurement errors are distributed as a Gaussian distribution with a fixed variance σ, i.e. they vary
around the “true“ value predicted by the model y = f ( x, a ) .
2. All measurement errors are statistically independent, i.e. the random error of one observation is independent
of the random error of another observation.

In this case the probability that an observation with a certain deviation from the predicted value is generated is:

 1  y i − f ( xi , a )  2 
Wi = exp −    ⋅ ∆y .
 2 σ  


53
On the basis of the statistical independence of the observations the total probability that all observations are
generated by the given model then is the product of the probabilities of each observation:

  1  y − f ( x , a ) 2  
W = ∏ Wi = ∏  exp  −  i i
  ⋅ ∆ y  .
i   2 σ   
i
  

Remark: The term ∆y is a constant term which turns the probability density (1st term of equation) given by the
Gaussian distribution into a probability. It is constant and can be taken as equal for all samples.

The merit function then results as

K (x, y , a ) = − log(W ) = ∑
N
( yi − f (xi , a ))2 − N log(∆y )
i =1 2σ 2

Since the scaling of the merit function is arbitrary (no change of the position of the minimum), the constant
factor 2σ2 and the constant 2nd addend -N⋅log(∆y) can be omitted. As expected, we obtain the merit function of
the Least-Squares method and the optimum set of parameters derived from it:

K (x, y , a ) = ∑ ( yi − f ( xi , a ))
N
2

i =1

a opt = argmin(K (x, y , a ))

a

Remark: This derivation demonstrates the statistical prerequisites underlying the Least-Squares method which
is often presented descriptively. They are often, but not always, fulfilled in measurements. Therefore, this
method should never be applied unchecked!

7.2.2 Chi-Square method

The Least-Squares method is generalized by admitting different variances for the different observations. Of
course, they must be known beforehand10. Thus, there are triplets (xi,yi,σi) for each sample. In this case, the
variance cannot be extracted as a constant factor in the formula for the Least-Squares method given above.
Hence, we obtain as the function to be minimized and optimum set of parameters:

 y − f ( xi , a ) 
2

( )
N
χ x, y, σ , a = ∑  i
2

i =1  σi  ,
( (
a opt = argmin χ x, y, σ , a
a
2
))
σ being a vector of the variances corresponding to the observations. The merit function is called the Chi-
Square. Since the Chi-Square method includes the Least-Squares method as a special case with constant
variance, we will deal with the Chi-Square exclusively in the following.

7.2.3 Robust Fit methods

The assumption that the random errors are normally distributed represents a rather “rigid” marginal condition.
The likelihood that a sample deviates from the model value by more than 5σ is infinitely low in this case.
Therefore, the Chi-Square method will fail, if there is one single “outlier” among the observations. The merit
function increases significantly, because the likelihood for this event is considered to be very low based on the
assumption of the Gaussian distribution. The minimizing procedure will attempt to decrease the costs, i.e. to
include the outlier in the model. Thus, its weight becomes too large in the merit function and the fit method does
not react to outliers robustly. If it is not sure, whether the Gaussian distribution of random errors can be taken as

10
For example, the random error may be proportional to the quantity y upon reading the measuring instrument.

54
a basis, it should be considered which other distribution may be more suitable. Having found one, we can
determine the likelihood and thus the merit function according to the above-described scheme.

E.g., if outliers have to be expected, we can use a distribution that does not approach zero as fast as the Gaussian
distribution. Then the likelihood of large deviations from the mean is not too low any longer. The Lorentz
distribution is the obvious choice for this purpose:

1
Wi = ⋅ ∆y
1  y f (x i , a ) 
2

1 +  i− 
2 σ 

It resembles the Gaussian distribution within the range ±2σ, but its slopes are not so steep. As a merit function
for the ML method based on the Lorentz distribution we obtain:

( )  
K x, y, a = − log(W ) = − log ∏ Wi 
 i 
,
N  1  y i − f ( xi , a )  2 
∝ ∑ log1 +   
 2  σ  
i =1 

The constant addend -N⋅log(∆y) has been omitted again. This merit function is significantly less sensitive to
outliers than the Chi-Square method and the Least-Squares method, respectively (cf. Exercises), however, it is
based on an assumption about the distribution of measuring errors (Lorentz distribution), which is only
empirically founded. At any rate, it is better to use theoretical arguments for choosing a distribution. For
example, the Poisson distribution can be assumed to be the theoretically correct distribution for the random error
of radioactive decay measurements (see Exercises).

Remark: If possible, the assumption of a Gaussian distribution should not be rejected carelessly, because it has
the advantage that steps 2. and 3. of the fitting process (statement of a parameter interval where the “true”
parameters are probably found and statement of fit quality, i.e. suitability of the model) can be solved
statistically more “neatly” (see section further below in this chapter).

7.3 Minimization of functions

The merit function being determined, the question arises of how to minimize it. Two cases are presented here,
which allow the majority of practical cases to be solved:

1. The Chi-Square is used as merit function and

a) the model f ( x, a ) is linear in the parameters a (not necessarily linear in the data!). In this case, the
minimization problem can be directly solved without time-consuming iterative algorithms (in the
literature this case is usually called ”general linear least-squares“)

b) the model f ( x, a ) is nonlinear in the parameters a. Then only iterative methods can be used to
minimize the function.
2. The Chi-Square is not used as a merit function, i.e. distributions of the random error other than the Gaussian
distribution are employed. In this case there is nothing left but iterative methods independent of whether the
model f ( x, a ) is linear in the parameters a or is not.

All cases are covered in the following.

7.3.1 ”General Least-Squares“ method

The model is assumed to be linear in the parameters a for this method:

55
 M
f ( x, a ) = ∑ a k f k ( x ) ,
k =1

f k ( x ) being an arbitrary function of the independent variable x. It can also be nonlinear, e.g. a power function

f k ( x ) = x k −1 .

Then the model function is a polynomial in x

f ( x, a ) = a1 + a 2 ⋅ x + a3 ⋅ x 2 + … + a M ⋅ x M −1 .

With the general formulation of a linear function in the parameters, the Chi-Square is:

2
 M

 y i − ∑ a k f k ( xi ) 
χ 2 (x, y , σ , a ) = ∑ 
N
k =1  .
i =1
 σi 
 
 

With the definitions

f j (x i )
"design matrix" A = {Aij } with Aij = (N × M ) − matrix
σi
 y1 σ 1 
  
vector of normalized observations b =  
 
yN σ N 
 a1 
 

parameter vector a =  
 
aM 

the Chi-Square can be rewritten as follows:

χ 2 (x, y, σ , a ) = A ⋅ a − b
2
.

The Chi-Square can thus be converted in this case such that it corresponds to the quadratic distance between the
vectors A a and b. The optimum parameter vector

( ( ))
a opt = argmin χ 2 x, y, σ , a = argmin A ⋅ a − b
a a
( 2
)
exactly corresponds to the Least-Squares solution of the system of equations A a = b, because this solution just
minimizes the quadratic distance between A a and b. The solution is known using SVD:

a opt = V ⋅ diag (1 w j ) ⋅ U ⋅ b ,
T

A = U ⋅ diag (w j ) ⋅ V
T
being the SVD of the matrix A.

56
Using SVD it is possible to solve the minimization problem in one step without iteration. Therefore, this method
is particularly efficient and the absolute minimum is found. Additionally, any possible numerical problems can
be solved by treating the singular values (cf. chapter on systems of equations: Setting the reciprocal value to
zero, if the singular value is too small).

Remark: The so-called normal equations are often found in textbooks as the solution of the “General-Least-
Squares”. They result from setting the derivatives of the Chi-Square to zero with respect to the parameters.

! ∂χ 2
0= ∀k = 1, … , M
∂a k
! M
1  N 
⇒ 0=∑ y
2  i
− ∑ a j f j ( xi ) f k ( xi ) ∀k = 1, … , M
i =1 σ i  j =1 
(A A)⋅ a = A
!
T T
⇔ ⋅b

The equation in the third row describes the M normal equations (second row) in a matrix form, using the above-
mentioned definitions of the design matrix, the parameter vector, and the vector of the normalized data. The
optimum parameter vector results from matrix inversion:

(
a opt = A A
T
) −1 T
⋅ A ⋅b .

By this method an extremum of the Chi-square function is found. It has to be verified that it is really a minimum,
which means an additional calculation step. Altogether, a solution via SVD appears to be more suitable than a
solution via normal equations, because the former directly yields a minimum and is numerically very robust.

7.3.2 Methods for nonlinear fitting: The ”Simplex“ method

If the model is nonlinear in the parameters and/or if the merit function is not the Chi-square (cases 1.b. and 2.),
there is nothing left but iterative methods to find the minimum in “mountains” of the merit function. There exist
several methods which search the minimum proceeding from an initial set of parameters. The simplest method is
the Newtonian gradient method again. Its disadvantage is that the partial derivatives of the model function f have
to be calculated (which is not always feasible) and that the method, depending on the initial values, certainly
ends up in the nearest minimum, which may be a local minimum and not the absolute minimum. The Simplex
method is presented below as one of the methods that neither require a partial derivative nor inevitably end up
in the nearest local minimum.

Remark: In Matlab, the function FMINSEARCH implements the Simplex method.

Remark: The Simplex method and related methods can probably ”jump across“ the nearest local minimum.
None of these iterative methods guarantees that the absolute minimum of the merit function is found!

The Simplex method works as follows:

1. Generate a Simplex in the N-dimensional parameter space. A Simplex generally is a geometrical figure with
(N+1) interconnected points (or vertices) in the N-dimensional space, e.g. a triangle in two dimensions.
Distribute the vertices of the Simplex around an initial value of the parameters in the parameter space.
2. Calculate the merit function K on the vertices of the Simplex.
3. Depending on the result, apply one of the elementary operations described below, which cause a move of
the Simplex in the parameter space.
4. Continue with 2., unless the termination criteria are fulfilled. The termination criteria are generally fulfilled,
if the values of the merit function on the vertices of Simplex do not differ by more than a predetermined
difference or if a given maximum number of iterations has been reached.

The following elementary operations can be performed:

1. Moving the vertex where the function is largest through the opposite face of the Simplex (reflection). Thus,
the Simplex moves away from the maximum.

57
2. Like 1., however, with additional contraction of the Simplex in the dimension perpendicular to the
reflection face.
3. Like 2., but expansion instead of contraction.
4. Contraction of the Simplex in one or more dimensions.
5. Expansion of the Simplex in one or more dimensions.

With each iteration the Simplex changes its position and size in the parameter space due to the elementary
operations and thus moves “amoebia-like” in the direction of the minimum. The Matlab function
FMINSEARCH uses empirical rules to choose the elementary operations depending on the distribution of the
values of the merit function on the corners of the Simplex. The following figure illustrates the rules for the
example of a tetrahedron (Simplex in three dimensions). Starting out from a current position of the Simplex, it is
calculated what the Simplex looks like in the next step:

58
7.4 Error estimation and confidence limits
The optimal set of parameters aopt which corresponds to the minimum of the merit function is determined by
establishing and minimizing the merit function. Now the question of the error, or uncertainties, in a set of
estimated parameters arises, i.e., how far does the estimated set of parameters aopt deviate from the “true” set of
parameters atrue. For this purpose confidence regions are usually calculated stating a region within the parameter
space, in which the true set of parameters is found with a predetermined probability. The limits of the confidence
region are the confidence limits. For calculating the confidence regions, again different cases are to be
distinguished (as for minimizing the merit function):

1. The Chi-square is used as a merit function and

c) the model f ( x, a ) is linear in the parameters a (”general linear least-squares“). In this case the
confidence regions can be stated directly from the SVD solution.

d) the model f ( x, a ) is nonlinear in the parameters a. Then the confidence regions must be calculated
via repeated calculation of the Chi-square in the neighbourhood of the optimal set of parameters (iso-
contours of the merit function).
2. The Chi-square is not used as merit function, i.e., distributions of the measurement error other than the
Gaussian distribution are used. In this case only statistical methods can be applied to determine the
confidence regions, which require the optimal set of parameters to be calculated from “synthetic”data sets
(bootstrap method).

The different methods are explained in the following.

7.4.1 Chi-square fitting: Confidence regions of parameters

7.4.1.1 Model function f ( x, a ) is nonlinear in the parameters

In this case the χ2 function is evaluated on a grid within the neighbourhood of the optimum set of parameters and
the iso-χ2 curves are determined which correspond to a change of the χ2 as compared to the minimum by
different ∆χ2. The following table shows ∆χ2 values for several probabilities p of the true set of parameters to be
found within the region delimited by the related iso-curve and for different numbers of parameters (without
proof):

M
p/% 1 2 3 4
68.3 1.00 2.30 3.53 4.72
95.4 4.00 6.17 8.02 9.70
99 6.63 9.21 11.3 13.3

Example: The probability that the true set of parameters is found within the region delimited by the iso-curve
with ∆χ2 =6.17, i.e. the curve on which the value of χ2 is by 6.17 larger than in the minimum, is 95.4% with two
parameters (M=2). This region is called the confidence region on the confidence level p = 95.4%.

Remark: There are many statements about confidence intervals in the literature. They represent the projection
of confidence regions on the axes of the parameter space. The ∆χ2 values given in the table for M=1 are always
valid for the confidence intervals (without proof). Example: A confidence interval lies on the 99% level, if it has
been generated by projection of the confidence region with ∆χ2 =6.63, independent of the dimension M of the
parameter space.

Remark: The confidence regions are always ellipsoids in this case (see below).

7.4.1.2 ”General linear least squares“ case

In this case an analytical formula for the limits of the confidence regions can be given, i.e. it is not necessary to
extract the iso-curves by calculating the Chi-square on a grid. The table given above is still valid for the
confidence levels, however, there is an easier way to calculate the limits. If a variation δa of the parameter

59
vector proceeding from the optimum vector aopt is given, the change of the Chi-square is calculated according to
the following formula (without proof):

( ( ))
M
∆χ 2 = ∑ wi V (i ) ⋅ δ a
2 2

i =1

with
V (i ) : i - th column of the matrix V
δ a : variation of the parameter vector proceeding from the optimum set of parameters a opt

Thus, the columns of the matrix V from the SVD of the design matrix A form the axes of an ellipsoid
representing the margins of the confidence regions (iso-χ2 curves). The axis intercepts for the iso-curve at ∆χ2 =
1 just represent the variance of the parameters (without proof):

2
V
M

σ (a k ) = ∑  ki
2
 , k = 1...M
i =1  wi 

These values represent the confidence interval of the parameter ak on the 68.3% level.

7.4.2 The “Bootstrap“ method for estimating the confidence regions

In cases of non-normal distribution of measurement errors (merit function is not the Chi-square) the above-
mentioned relation between confidence regions and confidence levels are not valid. In that case only empirical
statistical methods can be applied to calculate the confidence regions, because analytical calculations are usually
not feasible. A suitable empirical method is the "bootstrap“ method, which is explained in the following: From
the available data set, new synthetic data sets are generated which have the same properties as the original data
set in a statistical sense. Such a synthetic data set can be generated by drawing N data pairs (with replacement!)
from the set of samples (N measured pairs). By drawing with replacement the data pairs of the synthetic data set
are not identical with the measured data set and an arbitrary number of different synthetic data sets can be
generated. Model fitting is then performed for each synthetic data set and the optimum set of parameters aopt is
determined which, of course, differs depending on the data set. Now the distribution of the determined sets of
parameters is calculated and those regions are stated as confidence regions in which a certain percentage of sets
of parameters are found.

Remark: The ”bootstrap“ method is only applicable, if the sequence of samples is not relevant and if the
measurement errors of all samples have the same (non-normal) distribution. If this is not the case, time-
consuming Monte-Carlo methods are to be applied, which are not treated here.

7.5 Goodness of fit

Having determined the optimum set of parameters and its confidence ranges, we have to investigate whether the
model sufficiently describes the data (“goodness of fit“). The probability Q that the remaining deviations of
samples from the model prediction result from random measurement errors is determined and it is investigated,
whether the deviations are of a systematic nature. Systematic deviations mean that the model does not cover
systematic effects in the measured data, i.e., its prediction capacity and quality are limited.

The goodness of fit is evaluated on the basis of the probability Q that the remaining
deviations of measured values from the model prediction result from random
measurement errors. The larger Q is, the larger is the goodness of fit.

If the Chi-square merit function is applied, results of mathematical statistics can be used in order to calculate Q.
For this, we make use of the fact that a K-fold sum of squares of independent N(0,1)-distributed random
variables (Gauss-distributed with mean 0 and variance 1) are χ2K –distributed (Chi-square-distributed with K
degrees of freedom). The Chi-square does not exactly correspond to this condition, because the model function
reduces the independence of random variables. Nevertheless, the Chi-square is assumed with good
approximation to be χ2(N-M) –distributed with N samples and M model parameters. Q can be derived from this as
follows (without proof):

60
 χ2 2
1
Q = 1− ⋅ ∫y
( N − M ) 2 −1
exp(− y )dy
Γ((N − M ) 2 ) 0

= 1 − gammainc (χ 2 2 , (N − M ) 2 )

The second addend represents the incomplete Gamma function. Q is directly calculated by inserting the
observed Chi-square values at the minimum of the merit function as well as the number (N-M) of degrees of
freedom into this formula.

Remark: The incomplete Gamma function is available in Matlab as the command gammainc.

On the basis of the value of Q the following statements on the model can be made:

Value of Q Conclusion
small (Q~10-18) The model is incorrect or
the variance of random measurement errors has been estimated too small or
the assumption of Gauss-distributed measuring errors is not correct.
Q~10-3 Model OK
Q~1 Too good to be true. Maybe the variance of random measurement errors has been
estimated too large.

Thus, it is obvious that the assumptions on the statistics of measurement errors must be correct for an
“applicable” evaluation of the value of Q.

The value of the Chi-square can also be assessed with a ”rule-of-thumb”. Assuming that the random deviation
of data from the predicted model values is approximately 1σ on average, then each addend in the Chi-square
function is approximately 1. The Chi-square thus takes approximately the value N for N addends (N samples).
Simultaneously, the adjustment becomes easier with increasing number of parameters. This reduces the expected
value of the Chi-square by about the number of parameters M. The rule-of-thumb thus reads, that the value
of the Chi-square should about equal the number of degrees of freedom N-M, if the remaining deviations
of samples from the model prediction are of random nature and to make the model acceptable. Hence, we
obtain the following evaluation table:

Value of χ2 Conclusion
χ2>>N-M The model is incorrect or
the variance of random measurement errors has been estimated too small or
the assumption of Gauss-distributed measuring errors is not correct.
χ2~N-M Model OK
χ2<<N-M Too good to be true. Maybe the variance of random measurement errors has been
estimated too large.

Remark: For the other merit functions which are not based on a Gauss-distribution of measuring errors, Q and
thus the goodness of fit are difficult to determine.

61
7.6 Summary of methods
The following table shows the different cases of model adjustment:

Case Merit function Minimization Confidence ranges Goodness of

fit
Random measurement Chi-square via SVD or normal Calculation of iso- Chi-square
errors Gauss-distributed equations χ2-curves via SVD distribution or
and model function rule-of-thumb
linear in the parameters
(”general least-squares“)
Measurement errors Chi-square Simplex method Calculation of iso- Chi-square
Gauss- distributed, but χ2-curves by distribution or
model function nonlinear sampling the rule-of-thumb
in the parameters parameter space
Measurement errors not Merit function Simplex method “Bootstrap“ method Only feasible if
Gauss-distributed according to the the probability
Maximum-Likeli- distribution of
hood principle based the values of
on the assumed the merit
distribution of function can be
measurement errors derived,
(e.g. Lorentz perhaps rule-
distribution of-thumb is
(“robust“ fit) or applicable
Poisson distribution)

62
7.7 Exercises

1. The supplied script linfit.m fits a straight line to test data. The χ2-function is used as merit function and the
Simplex method for function minimization (Matlab function FMINS).

i. Try to understand the script linfit.m and alter the test data as well as the initial values of the
parameters. Which influence do “outliers” have?
ii. Alter the merit function in such a way that it is based on Lorentz-distributed measurement
errors. How does the influence of outliers change?

2. In order to measure the half-life value of a radioactive isotope, the number of decays is counted every 10s
for an interval of ∆t=1s. The following values are measured:

Time / Aktivity Time / s Aktivity

s (Number of (Number of decays)
decays)
10 24 60 5
20 17 70 4
30 11 80 2
40 10 90 3
50 6 100 1

i. Estimate the half-life. Start from the model that the activity decreases exponentially:
A( t ) = A0e −αt
Moreover, assume that the samples are Poisson-distributed:
λk
p( k , λ ) = e− λ
k!
k being the number of decays and λ the expected number of decays (the expected value is the product
of activity and counting interval: λ=A(t)* ∆t)
(Hint: The merit function to be minimized is available on request as Matlab script f_pois.m)

ii. By which percentage does the estimated half-life value change, if Gauss-distributed samples are
assumed (least-squares fit)?

3. The following figure shows the measured values of CO2 concentration in Hawaii from 1981 on (registered
two-weekly). The standard deviation of the measurements is about σ = 0.16 ppm.
C o 2 -K o n ze n tra ti o n i n M a u n a L o a , H a w a i i , g e m e s s e n 1 4 -tä g i g a b 1 9 8 1
356

354

352

350

348
C O2 / ppm

346

344

342

340

338

336
0 50 100 150 200 250
Z e it a b 1 9 8 1 / 1 4 T a g e

Hint: The data are found in the file mauna.dat. Use the command „load mauna.dat –ascii“ for loading the
data into Matlab’s workspace.

63
a) Determine the rate of increase in the CO2 concentration in ppm/year. As a model use a linear increase in
concentration. When will it be 15% higher than that of 1981 according to this model?
Use the Chi-square (σ = 0.16 ppm for all samples) as merit function and the Simplex method for
minimization.
b) Like a), however, under the model assumption that the concentration increases quadratically.
Hint: The linear model function used in the script linfit.m can be easily extended to polynomials of the
n-th order by simply stating n+1 initial values of the parameters.
c) Assess the goodness of fit of the models used in a) and b) by applying the “rule-of-thumb for the
Chi-square. Can we rely on the results from a) and b) on the expected time of a 15%
increase of the concentration?
d) Which simple extension of the quadratic model function would improve the model? Repeat the
adjustment with this extended model function and assess the goodness of fit of this model using the
“rule-of-thumb”.
e) (*) Repeat the adjustment as performed in exercise d) using the general least-squares“ method (with SVD).

64
8 Analog-Digital Transform and Discrete Fourier Transform
Numerical methods play an important role in the analysis of sampled signals (e.g. the time course of temperature
and pressure or voltages and currents as functions of time). Among these methods, spectral analysis, i.e. the
analysis of the frequency contents of a signal by means of Discrete Fourier Transform (DFT) is particularly
worth mentioning and will be explained in the following. First, the process of Analog-Digital Transform (A/D-
Transform) is described which is used to discretize analog signals in such a way that they can be processed on
computers. On that basis the DFT as well as the Fast Fourier Transform (FFT) as a fast algorithm for calculating
the DFT are introduced. Finally, the filtering of signals is described as an application of DFT. The so-called
convolution theorem is of special importance in this context.

8.1 Analog-Digital Transform (A/D Transform)

Physical signals are generally continuous in time and value. For example, the time course of the pressure can
generally be written as a real-valued function of continuous time p = p(t). We speak of analog signals and
continuous functions, respectively11. Since computers can only process finite quantities of numbers, these
measured signals/functions have to be transformed into finite sequences of numbers:

x (t ) → {x (n )} ; n ∈ K ⊂ .

The sequence consists of sampled values of the function at multiples of a predetermined period T:

x (n ) := x (n ⋅ T ) ; n ∈K ⊂
x (n ) : Samples
T : Sampling period
fs = 1 T : Sampling period

This process of transformation into a sequence of numbers is called sampling or discretization. Its most
important quantity is the sampling frequency fs, the reciprocal value of which is the sampling period T. The
technical realization of sampling in order to discretize physical signals is called Analog-Digital Transform (A/D
Transform) or digitization.

0.6

0.4

0.2

0
Amplitude

-0.2

-0.4

-0.6

-0.8

-1
0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5
Zeit / s

Analog-Digital Transform: Analog signal x(t) (solid line), sampling points (dashed lines) and
samples (rhombs). The samples form the sequence of numbers {x(n)} = { x(0), x(1), x(2) , x(3)
, . . .} related to the function x(t).

11
In case a real physical quantity is given as a function of time, the term “signal” is used in the following. The
term “function“ is used as an abstract mathematical representation of a signal.

65
Remark: Sampling is mathematically also feasible for infinite signals, however, on the computer it can be
applied to finite signals only, otherwise the sequence of samples would become infinitely long.

Remark: Because of the finite precision of the number representation on the computer the values of the samples
do not exactly correspond to the values of the function, but have a random roundoff error depending on the
number representation. The deviation of the value of the sample represented on the computer from the true value
of the function at the related point of time is called quantization error. If, for example, a sine function is
represented by a series of samples, the quantization error just corresponds to the precision with which the values
of the sine function are calculated.

The course of the function between the sampling points is lost by sampling, however, the so-called sampling
theorem states that no information is lost under a certain condition:

Sampling theorem (without proof):

A function x(t) can be uniquely reconstructed from the sequence of its samples {x(n)}, if it
comprises only frequencies below half of the sampling frequency. In other words: The sequence
of samples is a unique representation of the function only under this condition.

To illustrate the latter, the signal cannot alter arbitrarily fast, if high frequencies are lacking. Therefore, it
becomes predictable in between two samples (can be interpolated).

Remark: In order not to infringe the sampling theorem, the analog signals must be low-pass filtered before
sampling, i.e. all frequencies above half of the sampling frequency are blocked and only the low frequencies are
passed.

Remark: If a signal contains frequencies higher than half of the sampling frequency, they will be reflected at
lower frequencies in the sequence of samples (“mirror frequencies”). This phenomenon is called aliasing (see
Figure below and the Exercises).
1

0.8

0.6

0.4

0.2
Amplitude

-0.2

-0.4

-0.6

-0.8

-1
0 10 20 30 40 50 60 70
Zeit / ms

The aliasing phenomenon: Sampling of a sine wave of high frequency (solid line) at a
sampling frequency which is too low. The samples (rhombs) also represent a sine wave of a
lower frequency (dotted line) so that an ambiguity occurs (aliasing). The lower frequency is
the ‘mirror frequency’ to the original frequency.

8.2 Diskrete Fourier Transform (DFT)

The Fourier Transform is defined for continuous functions in general. It describes an expansion of the function
into sine and cosine functions comprising any possible frequency (continuous frequency spectrum). It states at
which amplitude and phase each frequency is contained in the function. In order to obtain a version of the
Fourier Transform which uses finite sequences of numbers only and thus can be performed on the computer,
sampled periodic functions are considered in the following. Let us assume the sampling period to be T and the
period of the signal T0 = N⋅T, with N being a positive integer. Then such a function is unambiguously
represented by the N samples of a period

66
 {x (n )} ; n = 0, …, N − 1

x(N) is again identical with x(0) due to periodicity etc. This finite sequence of numbers can be processed on a
computer and nevertheless represents the entire infinitely long periodic function.

What does a frequency spectrum of a sampled periodic function look like? It can only be composed of sine
functions, which also have the period T0. Otherwise, a sine function would not be identical in two successive
periods of the function, which contradicts the assumed periodicity. All sine functions with the period T0 / k with
k as an arbitrary integer number also have the period T0. Thus, the spectrum comprises only frequencies fk = k /
T0; hence, it is discrete with the sampling period f0 = 1 / T0 (this corresponds to the known Fourier series).
Therefore, one prerequisite for representing the Fourier integral on the computer is already fulfilled: The
spectrum can be represented as a sequence of numbers. But is this sequence really finite? This question is treated
in the next paragraph.

Owing to the symmetry of the Fourier Transform (transform and reverse transform are only distinguished by a
negative sign, which corresponds to a time reversal) the argument used above is valid universally as well: If the
frequency spectrum is periodic with the period fs = 1 / T (periodic with the sampling frequency), the time
function is discrete with the sampling period T (which is just the starting point). In case of frequency samples at
fk = k / T0 and a period of the spectrum of 1 / T = N / T0, there are exactly N different frequency components.

Altogether, the frequency spectrum is periodic and discrete for discrete periodic functions. The function as well
as its frequency spectrum can uniquely be represented as a finite sequence of numbers of the length N on the
computer by the values of samples of one period each. The Discrete Fourier Transform relates the samples of the
time function to those of the spectrum (without proof):

Discrete Fourier Transform of Length N :

N −1
 kn 
X (k ) = ∑ x (n ) ⋅ exp −2πi  ; k = 0 …N − 1
n =0
 N .
1 N −1
 kn 
x (n ) =
N
∑ X (k ) ⋅ exp 2πi N 
k =0
; n = 0 …N − 1

X(k) and x(n), respectively, represent the samples of the spectrum and the time function, respectively. The first
equation denotes the DFT and the second equation denotes the inverse DFT. The relation between the actual
index k and n, respectively and the frequency and time, respectively, is given by:

n = 0,1,2, …, N − 1 t = 0,T ,2T , …, (N − 1)T

1 2 (N − 1)
k = 0,1,2, …, N − 1 f = 0, , , …,
T0 T0 T0
fs f f
= 0,1 ,2 s , …, (N − 1) s
N N N
These formulas are applied to calculate the frequency and time axes of the sequences of numbers calculated by
the DFT. The sample values x(n) may be real- or complex-valued, physical signals being always real-valued.
The X(k) are complex-valued in both cases. The value |X(k)| represents the amplitude and the argument
(arg(X(k))) represents the phase (shift), with which a sine/cosine function of the corresponding frequency is
represented in the function x(n).

67
 1.2

0.8

0.6

0.4

0.2

-0.2

-0.4

-0.6

-0.8
0 0.1 0.2 0.3 0.4 0.5

DFT: Periodic physical signal (solid line, only one period is represented) and five of its
different frequency components. Note the different frequencies, amplitudes, and relative
shift (”phases“) of the components.

Remark: It is important to know that the DFT observes periodic functions exclusively. Considering an arbitrary
finite function or signal, the DFT does not calculate the spectrum of that finite signal, but the spectrum of the
periodically continued version of that period (for the consequences of this fact see the Exercises).

8.2.1 Real-valued time functions

If the time function x(n) is real-valued, the following relation holds for the spectral samples:
N
X (N − k ) = X * (k ) ; k = 0, …,
2
, * denoting the complex conjugation. Thus, N/2 complex values X(k) are sufficient in order to represent N real
values x(n) and the spectrum is uniquely represented within the range of up to half the sampling frequency (cf.
formula presented above for conversion between frequency index and frequency). Therefore, the spectrum is
only represented in the range between 0 Hz and half the sampling frequency fs/2 for real signals. The remainder
is redundant.

Remark: The sampling theorem and the aliasing mentioned above result from the periodicity of the spectrum. If
the spectrum is periodic with a sampling frequency fs, only one of the periods is physically significant. For
example, the values of the spectrum at the frequencies f = 0.3 fs and f = (0.3 fs + fs) are always identical. Thus,
the sampled signal cannot have comprised sine functions of those frequencies independent of each other, so they
are indistinguishable. For real signals this applies also to the frequencies f = 0.3 fs and f = (0.3 fs + fs/2) owing to
the property mentioned above, so that the sampling theorem follows.

8.2.2 Intensity spectrum

Often we are not interested in the complex spectral values X(k), but the question is: With which intensity is a
sine of a certain frequency contained in the signal? The intensity is just the square of the spectral value |X(k)|2.
Within the range of typical physical signals, it may fluctuate by several orders of magnitude across the
frequency so that it is often useful to take the logarithm for plotting the data. As intensity spectrum in dB
(decibel) the following quantity is defined:

(
L (k ) dB = 10 ⋅ log10 X (k )
2
) ; k = 0, …, N − 1

Remark: When the amount |X(k)| just increases tenfold, the value in the intensity spectrum changes by 20 dB.

68
8.3 Fast Fourier Transform (FFT)
The DFT is a relatively time-consuming method when calculating the spectral values X(k) directly according to
the formula of the DFT. For each of the N spectral values N multiplications and additions are to be calculated so
that the total expenditure is proportional to N2. However, Gauss already found out that partial sums of different
frequencies (index k) are identical, if N is a power of 2, i.e. N = 21, 22, 23, 24, 25, 26.etc. The FFT makes use of
these symmetry properties and calculates partial sums only once. Thus, the expenditure is reduced from N2 to
N⋅log2(N), which means a considerable difference for large N. Therefore, the FFT is nothing but a fast algorithm
for calculating the DFT. It can be assumed that the DFT would not have gained such an importance within the
field of computer-controlled signal analysis, if this fast algorithm for its calculation had not been found.

One possible derivation of the FFT algorithm is explained in the following. With the definition

 2πi 
WN := exp −
 N 
the DFT reads

N −1
X (k ) = ∑ x (n ) ⋅WNkn ; k = 0, …, N − 1 .
n =0

First, this sum is divided into two partial sums containing the even and odd samples, respectively:

N −1
X (k ) = ∑ x (n ) ⋅WNkn
n =0
N −1 N −1
2 2
= ∑ x (2n ) ⋅W
n =0
k 2n
N
+ ∑ x (2n + 1) ⋅W
n =0
k (2n +1)
N

N −1 N −1
2 2
2 kn 2 kn
= ∑ x (2n ) ⋅ (W )
n =0
N
+W −k
N
⋅ ∑ x (2n + 1) ⋅ (W )
n =0
N

With the identity (see definition)

WN2 = WN
2

we obtain:

N −1 N 2−1 
2

X (k ) = ∑ x (2n ) ⋅W kn
N
+ W ⋅  ∑ x (2n + 1) ⋅WNkn 
k

 n =0
N
2
n =0 2


The two remaining sums just represent the formula for the DFT of the even and odd samples, respectively,
which are referred to as Xg(k) and Xu(k) in the following:

N −1
2
X g (k ) = ∑ x (2n ) ⋅W
n =0
kn
N
2
; k = 0, …, N 2 − 1
N −1
.
2
X u (k ) = ∑ x (2n + 1) ⋅W
n =0
kn
N
2
; k = 0, …, N 2 − 1

For these DFT‘s of length N/2 the index k ranges between 0 and N/2 –1, however, for the spectral values X(k) it
ranges from 0 to N-1. Owing to the identity

69
 WN(k + 2)n = WNkn
N

2 2

the partial sums for k = N/2 up to N-1, however, again represent the DFT of the even and odd samples,
respectively. Thus, the modulo function mod yields:

(
X (k ) = X g (k )modN
2
) +W k
N (
⋅ X u (k )modN
2
) ; k = 0, …, N − 1

The unknown spectral values X(k) thus result from the two spectral values Xg(k) and Xu(k) by multiplication and
addition. Each Xg(k) and Xu(k) is used twice due to the modulo function. This is the time-saving fact.

Altogether, these conversions divided the DFT of the length N into two DFT’s of the length N/2 each and a
subsequent arithmetic combination. The two DFT’s work on the samples with even and odd indices,
respectively, and the subsequent combination calculates the result of the entire DFT from the results of the
partial DFT’s. The subsequent combination requires a total of N operations (one operation per spectral value).
The two partial DFT’s can be split in the same way again in a further step, which results in two DFT’s of the
length N/4 and a subsequent combination with N/2 operations each. Altogether, there are four DFT’s of the
length N/4 and two subsequent combinations with a total of N operations in this next stage. This splitting can be
repeated until the DFT works on two samples only. Then N/2 DFT’s of the length 2 with an expenditure of 2
each are required, i.e., a total of N operations again.

Thus, N operations are required in each stage (for the DFT’s of the length 2 in the lowest stage and for the
subsequent combinations in the following stages, respectively) and there is a total of 10 log2(N) stages (e.g. N=8
can be divided by 2 twice, until a DFT of the length 2 is reached. Thus, there are three stages altogether (one
stage with DFT’s of the length 2 and two stages with subsequent combinations)). The total expenditure of the
FFT is therefore proportional to N⋅log2(N), as stated above.

Remark: The time required for calculating the DFT can be reduced to a certain degree also for the lengths N
that do not correspond to a power of two. For this purpose N is split into a product of as many powers of two as
possible, a separate FFT is performed for these partial sums of samples and the whole thing is assembled
skillfully afterwards. This is called mixed-radix FFT, which is implemented in Matlab. The command fft
performs a pure FFT or a mixed-radix FFT depending on the length of the given vector. The expenditure for a
mixed-radix FFT ranges between that for a DFT and that for an FFT depending on how a number can be split iin
powers of two. As it is not known beforehand in most cases how well a number can be split, powers of two
should be used if possible. In order to obtain a desired sequence length, zeros should be added to the signal
which does not alter the spectrum.

8.4 Filtering and convolution theorem

In signal theory the term filtering means the change of the frequency spectrum of a signal/function. For
example, if certain frequency ranges in a physical signal to be analyzed are of particular interest, they may be
increased by filtering, while the unimportant ones, e.g. the frequency ranges covered by noise can be reduced.
Basically, four types of filters can be distinguished:

1. Low-pass filter: Passes only frequencies below a certain cutoff frequency.

2. High-pass filter: Passes only frequencies above a certain cutoff frequency.
3. Band-pass filter: Passes only frequencies within a range between a given lower and a higher cutoff
frequency.
4. „Notch filter“: Passes only frequencies beyond a range between a given lower and higher cutoff
frequency.

In the frequency domain the filtering process can be mathematically formulated in a simple way. Since the
spectral values state the amplitude and phase of the frequency content, filtering can be described as a
multiplication of spectral values by a frequency-dependent factor:

Y (k ) = X (k ) ⋅ H (k ) ; k = 0, …, N − 1

70
X(k) and Y(k) are the spectral values of the original and the filtered signals. The sequence of (complex) factors
H(k) is a unique definition of the filter and can be freely predetermined for realizing the required type of
filtering. It is called transfer function of the filter.

Which operation within the time domain corresponds to the multiplication in the frequency range, i.e., how must
x(n) and h(n) be linked in order that y(n) results (capitals and small letters mean pairs of time signals and
spectral values that belong together)? For this purpose the convolution theorem of the DFT is applied (without
proof):

{x(n)} and {h(n)} be two sequences of length N. Be y(n)} calculated by means of a cyclic convolution:

zyklic convolution :
y (n ) = (x ⊗ h ) (n )
N −1 ; n = 0, …, N − 1
= ∑ x (m ) ⋅ h ((n − m )mod N )
m =0

The convolution theorem states that the spectral values Y(k) of the sequence {y(n)} are just the product of the
spectral values of {x(n)} and {h(n)}:

convolution theorem :
y (n ) = (x ⊗ h ) (n ) ; n = 0, …, N − 1
⇔ Y (k ) = X (k ) ⋅ H (k ) ; k = 0, …, N − 1
mit : Y = DFT (y ) ; X = DFT (x ) ; H = DFT (h )

Thus, a (cyclic) convolution in the time domain corresponds to a multiplication in the frequency domain. This is
valid for the inverse, too: Multiplication in the time domain corresponds to conversion in the frequency domain.

Thus, the realization of a filter by multiplication in the frequency domain with the transfer function H(k)
corresponds to the cyclic convolution with the sequence {h(n)} in the time domain, which is calculated from the
H(k) by inverse DFT. The sequence {h(n)} is called impulse response.

Hence, filtering can be realized in two ways:

1. Applying the convolution theorem: First, calculate the DFT of the sequence, multiply it by the selected
transfer function and then form the inverse DFT, in order to obtain the time function of the output signal:

y = IDFT (DFT (x ) ⋅ H )

2. Generate the output sequence directly by convolution in the time domain.

y (n ) = (x ⊗ h )(n ) ; n = 0, …, N − 1

8.4.1 Convolution and cyclic convolution

In general, convolution for signals of any length is defined as:

y (n ) = (x × h )(n )
∞ ; n∈
= ∑
m =−∞
x (m ) ⋅ h (n − m )

For indices beyond the domain of x or h, a zero is inserted as value for this general definition. Therefore, the
sum may always range from −∞ to ∞ .

71
Since convolution is a time-consuming algorithm, it is often realized by means of the convolution theorem, first
switching to the frequency domain, then multiplying and finally transforming back into time domain. If the
signals have the length of a power of two, the FFT can be applied. In that case the total expenditure for forward
and inverse FFT as well as for the multiplication is lower than for a direct application of the convolution sum for
signal lengths exceeding N=128. For realizing the convolution using the convolution theorem the following has
to be considered: If the signals x and h are periodic and equal in length, the general form of the convolution is
identical with the cyclic convolution. The cyclic convolution yields exactly one period of the periodic signal
generated by convolution of two periodic signals. Since the DFT deals with periodic signals exclusively, the
convolution theorem of the DFT applies to cyclic convolution and not to the general convolution sum, even if
only finite signals are considered! The cyclic convolution realized by applying the convolution theorem and the
general convolution sum thus yield different results in general (cf. Exercises).

Remark: The general form of the convolution sum is often applied as a mathematical definition of a filter. This
also explains, why h is called the impulse response: If a pulse (x(1)=1, x(n)=0 otherwise), is used as input signal,
the output function y is identical to h.

8.4.2 Deconvolution
The inversion of the filtering process is called deconvolution, which is applied, if the filtering effect of a filter
needs to be compensated for (e.g. frequency-dependent alteration of the sensitivity of a sensor). Thus, an
impulse response h* is required which exactly neutralizes the effect of a known filter (impulse response h,
transfer function H):

be : y (n ) = (x × h ) (n )
desired : h * (n )
with : x (n ) = (y × h * ) (n ) = (x × h × h * ) (n )

With the help of the convolution theorem, h* is easy to determine:

DFT (y ) = DFT (x ) ⋅ DFT (h )

1
⇒ DFT (x ) = DFT (y ) ⋅
DFT (h )
1
⇒ x = IDFT (DFT (y ) ⋅ )
DFT (h ) ,
 1 
= y × IDFT  
 DFT (h )
 1 
⇒ h * = IDFT  
 DFT (h )

This equation shows that deconvolution is formally very easy to perform: For deconvolution, we simply divide
by the transfer function H in the frequency range. Thus H is cancelled. In the time domain this corresponds to a
convolution with the inverse DFT of the inverse transfer function 1/H. Division by H is, however, numerically
very problematic, if – as frequently occurring in practice – the value of H approaches zero for one or several
frequencies. There are special deconvolution methods, which are not treated here, which solve this problem by
approximation.

72
8.5 Exercises

1. Generate the different signals stated below with the sampling frequency fs = 10 kHz and a duration of 10 ms
(corresponding to100 samples). Calculate the intensity spectrums using the FFT (Matlab function fft) and
plot four periods (frequency axis -2 fs to +2 fs). Plot also four periods of each of the time signals (time axis –
20 ms up to +20 ms). Where do the noticable differences in the spectra come from (qualitative explanation)?

(a) Sine of the frequency 1 kHz.

(b) Sine of the frequency 925 Hz
(c) Sine of the frequency 1050 Hz

2. Aliasing phenomenon: Generate a sine of the frequency 1 kHz with the sampling frequencies fs = 1.5, 1.7,
1.9, and 2.1 kHz and a duration of 50 periods (50 ms) each. Calculate the respective intensity spectrum and
plot four periods (frequency axis -2 fs up to +2 fs). Identify the aliasing components. According to which
formula can their frequency be calculated?

3. Convolution theorem: The following sampled functions be defined:

x 1 = {0, 0, 0,1, 0, 0, 0, 0, 0, 0}
x 2 = {1,2, 3, 4, 5, 0, 0, 0, 0, 0}
x 3 = {1,2, 3, 4, 5, 6, 7, 8, 9,10}

Calculate the convolution x 2 ∗ x 1 and x 3 ∗ x 1 by insertion into the convolution sum (non-cyclic
convolution) and by application of the convolution theorem of the DFT (Hint: Matlab provides the functions
fft and ifft for the Discrete Fourier Transform). What strikes us about the result and how can it be explained?

73
9 Partial Differential Equations
While ordinary differential equations (ODE) describe functions of one variable, e.g. x=x(t), partial differential
equations (PDE) deal with functions of several variables, e.g. the temperature as a function of space and time, T
= T(x,t). PDE are generally subdivided into three classes. Unfortunately, there are no universally applicable
numerical methods for solving PDE like, e.g., the Runge-Kutta method for solving ODE. Thus, the three classes
will be treated in the following by means of important examples and their specific solution methods.

9.1 Classification of partial differential equations

9.1.1 Parabolic equations

Parabolic PDE describe diffusion processes in a broadest sense, for example the heat equation

∂ ∂
T (x , t ) = κ 2 T (x , t ) .
∂t ∂x

(with: T, temperature, x and t, variables of position and time and κ, coefficient of thermal conduction) and the
time-dependent Schrödinger equation from quantum mechanics

∂
i Ψ (x , t ) = H Ψ (x , t )
∂t
with the Hamilton operator H.

9.1.2 Hyperbolic equations

Hyperbolic PDE describe transport processes, for example the advection equation

∂a (x , t ) ∂a (x , t )
= −c
∂t ∂x
and the wave equation

∂2A (x , t ) 2 ∂ A (x , t )
2
= c .
∂t 2 ∂x 2
A is the amplitude of the wave, x and t are variables of position and time, and c is the velocity of the wave.

9.1.3 Elliptic equations

The prototype of an elliptic PDE is the Poisson equation of electrostatics describing the potential of a charge
distribution:

∂2Φ (x , y ) ∂ 2Φ (x , y ) 1
2
+ 2
= − ρ (x , y )
∂x ∂y ε0

Φ is the potential, x and y are position variables, and ρ is the charge density distribution. If the charge density is
zero, the resulting equation is called Laplace equation.

Remark: All equations are written one- or two-dimensional here for clarity. Formulation in several dimensions
is possible (see textbooks).

Remark: The equations stated here are prototypes of the three classes. Real problems are often represented by
mixtures of these prototype equations, e.g., the wave equation with damping.

74
9.2 Initial value problems
Equations that are position- as well as time-dependent are treated as initial value problems. The state at time t =
0 must be known, just as in the case of ordinary DE’s. For the heat equation, for example, this means

T0 = T (x , t = 0)

and for the wave equation (2 integration constants, because it is an equation of 2nd order)

dA (x , t = 0)
A0 = A (x , t = 0) and A0' = .
dx
T(x,t) and A(x,t), respectively for t>0 are then to be calculated. To determine the initial state, however, is not
sufficient for solving the problem. In addition, the boundary conditions must be formulated. For this purpose,
the solution is restricted to the range of x-values
x ∈ − L 2 , L 2
 
with arbitrary and fixed L. There are several ways to restrict the solution to the boundary (stated for the heat
equation here):

1. Fixed values of the solution on the boundary (Dirichlet’s boundary condition):

T (x = − L 2 , t ) = Ta , T (x = L 2 , t ) = Tb

2. Periodic boundary conditions:

T (x = − L 2 , t ) = T (x = L 2 , t )
or
dT dT
=
dx x =−L 2 dx x =L 2

Thus, the solution is searched within an area limited by t=0 and x=± L/2. The question of why boundary
conditions are needed is solved later on.

9.2.1 Discretization
For numerical calculation, time variable as well as position variable are discretized by sampling. The sampling
points are

tn := n ⋅ τ ; n = 0,1, 2, 3, 4, …
with the sampling period/step size τ. The area x=± L/2 is divided into N intervals, i.e., the sampling period/step
size of the position variable is

L
h=
N
and the sampling points (number of points: N+1) are

x i := i ⋅ h − L 2 ; i = 0,1,2, …, N .

The solution is then searched for the variable pairs (xi , tn), i.e. the solution points

Ti n := T ( xi , tn )

75
are to be determined (here written for the heat equation). These are the definitions of sampling points and
sampling positions used to describe the different methods in the following.

9.2.2 Heat equation: FTCS scheme

The FTCS scheme (Forward-Time-Centered-Space) makes use of the right-hand derivative formula for
discretizing the time variable (Forward Time) and the centered derivative formula for the position variable. For
the heat equation this means that the derivatives are being replaced by the following discrete approximations:

∂T (x , t ) T (x i , tn + τ ) − T (x i , tn ) Tin +1 − Tin
→ =
∂t τ τ .
∂ T (x , t )
2
T (x i
+ h , tn ) + T (x i
− h , tn ) − 2 ⋅ T (x ,
i n
t ) T n
+ T n
− 2 ⋅ T n

2
→ 2
= i +1 i −1
2
i

∂x h h
Note that the index notation introduced above was used for time and position index, respectively. Insertion into
the heat equation yields

n n n
Tin +1 − Tin T + Ti −1 − 2 ⋅ Ti
= κ ⋅ i +1
τ h2 .
κτ n
⇒ Ti = Ti + 2 (Ti +1 + Ti −1 − 2 ⋅ Ti )
n +1 n n n

h
0 n n n
From the initial values Ti and the boundary values T0 and TN the solution values Ti can be calculated
iteratively with this formula.

Remark: It is necessary to predetermine the boundary values for the centered derivative!

The step sizes τ and h have to be chosen according to the length and time scales occurring in the individual
physical problem. The relation between both values is of special importance, because it contradicts physical
intuition to calculate the solution with a large time step τ from closely neighbouring points (h relatively small).
The maximum time step for given h is (without proof):

h2
τ max =
2κ
If the time step is chosen larger than this maximum, the solution becomes unstable. This equation can be
empirically derived (see Exercises) or it can be derived by means of the von-Neumann stability analysis (cf.
Section 9.2.5).

9.2.3 Time-dependent Schrödinger equation: Implicit Scheme (Crank-Nicholson)

The time-dependent Schrödinger equation

∂
i Ψ (x , t ) = H Ψ (x , t )
∂t

with the wave function Ψ has the Hamilton operator

2
∂
H =− + V (x ) .
2m ∂x 2
for a particle of the mass m and the potential V. The square of the wave function describes the sojourn
probability density

2
P (x , t ) = Ψ (x , t ) .

76
With this FTCS scheme the discretized Schrödinger equation reads:

Ψ nj +1 − Ψ nj 2 Ψ nj +1 + Ψ nj −1 − 2 Ψ nj
i =− 2
+ Vj ⋅ Ψ nj
τ 2m h
i τ  
n
Ψ 2 + Ψ nj −1 − 2Ψ nj
⇔ Ψ nj +1 = Ψ nj − −
j +1
+ Vj ⋅ Ψ nj 
 2m 
2
h

using the above-mentioned notation for time and position indices. Writing the wave function as a time-
dependent column vector

 Ψ n 
 1 
 
Ψ :=  
n
 
 Ψ n 
 N 

we can write the discretized equation as a matrix equation solving the scheme for all sampling points of the
position simultaneously:

iτ
Ψ n +1 = Ψ n − H ⋅ Ψn
(*).
 iτ 
⇔ Ψ n +1
= 1 − H  ⋅ Ψ n
 

H is the matrix of the discretized Hamilton operator with the components

2 δi, j +1 + δi, j −1 − 2δi, j

H i, j = − + Vi δi, j .
2m h2
The matrix equation (*) solves the Schrödinger equation according to the FTCS scheme. Since such matrix
formulations are often found in the literature, this form is emphasized here.

Proceeding from the equation (*) the Crank-Nicholson method was developed as the most important type of the
so-called implicit schemes. On the right side of the equation not only the old value of Ψ, but the mean of the old
and new values is inserted:

iτ
H ⋅ ( Ψ n +1 + Ψ n )
Ψ n +1 = Ψ n −
2
.
 i τ  n +1  i τ  n
⇔ 1 + H  Ψ = 1 − HΨ
 2   2 

Solving for Ψ n+1 then yields the Crank-Nicholson scheme:

 i τ −1  i τ  n
Ψ n +1
= 1 + H  1 − HΨ .
 2   2 

As compared to the FTCS scheme, this scheme has the advantage that it is always stable and that the applied
matrix operators are unitary (without proof).

77
9.2.3.1 Wave package of a free particle
As an example of solving the Schrödinger equation we observe a Gaussian wave package which is used in
quantum mechanics to describe a free particle (potential V=0). As initial value for the wave function

1
Ψ (x , t = 0) =
σ0 π
( 2
exp (ik 0x ) exp − (x − x 0 ) 2σ02 ).
is used, x 0 being the current average position of the particle. σ0 is the current width of the package and
p0 = ⋅ k0 is the current impulse of the particle. The time development of the wave function can be calculated
from the initial state by means of the Crank-Nicholson scheme. We know from quantum mechanics that there is
also an analytical solution for the wave package. It reads:

1 σ0
( ( )) ( ( ) )
2
Ψ (x , t ) = exp ik0 x − p0t 2m exp − x − x 0 − p0t 2m 2α2
σ0 π α .
2 2
with : α = σ +i t m 0

Thus, the form of the wave package is preserved in time. The mean value (expected value) moves at the velocity
p0 m :

p0
x = ∫ x ⋅ P (x , t )dx = x 0 +
m
t

and the standard deviation expands according to

4
 α  2 2
σ (t ) = σ0   = σ 1 + t .
 σ  0
m 2σ04
 0

This analytical solution can be used in order to test the Crank-Nicholson scheme (see Exercises).

9.2.4 Advection equation: Lax-Wendroff scheme

We proceed from the wave equation

∂2A (x , t ) 2 ∂ A (x , t )
2
=c ⋅ .
∂t 2 ∂x 2
As in the case of ODE we rewrite this equation of the 2nd order into two equations of the 1st order. For this, we
define the intermediate variables

∂A ∂A
p= ; q =c⋅ .
∂t ∂x
With that we can write the wave equation as a pair of equations

∂p ∂q ∂q ∂p
=c⋅ ; =c⋅
∂t ∂x ∂t ∂x
The first equation has resulted from the insertion of the variable into the wave equation. The second equation is
yielded by a comparison of the mixed derivatives from p to x und from q to t. These two equations can be
written in a vector form again:

78
 ∂a (x , t ) ∂a (x , t )
= c ⋅B ⋅
∂t ∂x
 
 p (x , t ) 0 1
with : a (x , t ) =   ; B = 
 
q (x , t ) 1 0

The wave equation is not the simplest equation of this kind. Setting

−1 0
B =  
 0 0

, we obtain the so-called advection equation12

∂a (x , t ) ∂a (x , t )
= −c
∂t ∂x

, with a being the unknown quantity which we can imagine as the amplitude of a wave. In the following we
observe this simplest form of a hyperbolic equation. By insertion it is easy to demonstrate that each function of
the form:

a (x , t ) = f0 (x − ct )

with an arbitrary function f0 is a solution of the advection equation, c representing a diffusion velocity. For
example, a cosine-modulated Gauss pulse

 2π   (x − c ⋅ t )2 
 
a (x , t ) = cos  (x − c ⋅ t ) ⋅ exp −
λ   2σ 2 

is a solution. The function represents a short wave package of the wave length λ and of the spatial extent σ,
which moves in the direction x at the velocity c.

Since the analytical solution is known in the case of the advection equation, it can be used to test numerical
methods for solving hyperbolic equations. In the following the Lax-Wendroff scheme is presented for this
purpose.

Remark: The advection equation represents the simplest form of a conservation equation.

∂a ∂a ∂F (a )
= −grad (F (a )) oder 1-D: =−
∂t ∂t ∂x
a can be the density of the mass, impulse or energy and F(a) can be related to the flow of the mass, impulse or
energy, for example. The continuity equation for the mass in a current be stated as an example (one-
dimensional):

∂ρ (x , t ) ∂
=− (ρ (x , t ) ⋅ v (x , t )) .
∂t ∂x

ρ is the mass density and v the velocity of the current.

12
Advection: Horizontal movement of masses of air or horizontal flow of masses of water in the sea (according
to some dictionary).

79
9.2.4.1 Lax-Wendroff scheme
Theoretically, the advection equation can be discretized according to the FTCS scheme, however, this does not
yield a stable solution (see Section 9.2.5 and Exercises). The Lax-Wendroff scheme is of a higher order and
represents a stable method for solving the advection equation, although it is not without numerical errors (see
Exercises). To derive the scheme, we proceed from the Taylor expansion of the 2nd order for the time
development:

∂a (x , t ) τ 2 ∂ 2a (x , t )
a (x , t + τ ) ≈ a (x , t ) + τ ⋅ + ⋅
∂t 2 ∂t 2
Then the second derivative with respect to time is determined from the advection equation:

∂a ∂F (a )
=− with : F (a ) = c ⋅ a
∂t ∂x
2
∂a ∂ ∂F ∂ ∂F
⇒ 2
=− =−
∂t ∂t ∂x ∂ x ∂t
With

∂F (a ) dF ∂a ∂a ∂F (a )
= ⋅ = F ′ (a ) ⋅ = −F ′ (a ) ⋅
∂t da ∂t ∂t ∂x
we obtain the following form of the second derivative by insertion

∂ 2a ∂ ∂F (a )
2
= F ′ (a )
∂t ∂x ∂x
and hence for the Taylor development:

∂  τ 2  ∂ ∂F (a (x , t ))
a (x , t + τ ) ≈ a (x , t ) − τ ⋅  F (a (x , t )) + ⋅  F ′ (a (x , t )) 
 ∂x  2  ∂x ∂x 

This formula again allows the state at time t+τ to be calculated from states at time t. It is discretized now
(nomenclature of indices as for the heat equation). For the first derivative of time the centered formula is
applied:

∂ Fin+1 − Fi n−1
F (a ) →
∂x xi ,tn 2h

The calculation of the multiple derivative in the last term is rather complex. The outer derivative is a right-hand
formula, just like the inner derivative. The derivative F′ is calculated from the mean of the amplitude values
found in both terms of the outer right-hand derivative:

80
  ( )  ∂F (a )
F ′ (a ) ∂F a  − F ′ (a ) 
∂  ∂F (a )  ∂x  xi +1 ,tn
 ∂x  x ,t
F ′ (a )  → i n

∂x  ∂x  xi ,tn h
 ∂F (a ) Fin+1 − Fin
F ′ (a ) ∂x  → F ′ (ain+1 + ain ) 
2 ⋅
x ,t
  h
i +1 n

 ∂F (a ) Fin − Fin−1

 ′( )
F a ∂x  → F ′ (ain + ain−1 ) 
2 ⋅
x ,t
  h
i n

Altogether, this yields:

Fi n+1 − Fin−1 τ 2 1  n Fi +1 − Fi F n − Fin−1 

n n

ain +1 = ain − τ ⋅ + ⋅ Fi ′+1 2 ⋅ − Fi ′−n1 2 ⋅ i 

2h 2 h  h h 

with

Fin ≡ F (ain ) and Fi ′±n1 2 ≡ F ′ (ain±1 + ain ) 2 .

 

For the advection equation with

F (a ) = c ⋅ a
⇒ Fin = c ⋅ ain and Fi ′±n1 2 = c

the considerably simplified form reads:

cτ n c 2τ 2
ain +1 = ain −
2h
(ai +1 − ain−1 ) + 2 (ain+1 + ain−1 − 2ain )
2h

Just as in the case of the FTCS scheme, a sample at time t+τ can be calculated from three neighbouring points at
a time t in this way. Contrary to the FTCS scheme, however, it is stable if the time step is maximally τ max, with

h
τ max =
c
A proof is feasible by means of the von-Neumann stability analysis (cf. Section 9.2.5) (see also Exercises).

9.2.5 von-Neumann stability analysis

The von-Neumann stability analysis allows to determine stability rules for the different schemes for solving
PDE. We start from a little perturbation ”injected“ into the iteration equations and then track, whether this
perturbation diverges (instable scheme) or converges (“damping” of the perturbation, stable scheme). If the
perturbation is chosen skillfully, this can be decided analytically, i.e. without actual numerical iteration of the
equation. For the von-Neumann analysis the following special perturbation is used, which represents a fixed
wave with time-dependent amplitude:

81
 a nj = ξ n ⋅ e i ⋅k ⋅h ⋅ j
with :
ξ : Amplification factor
2π
k= : Wave number (arbitrary )
λ
h : spatial step
n : time index (left hand side )/ exponent (right hand side )
j : spatial index
i : complex number -1

The dependence on position (complex e-function) and the dependence on time (power function) are separated
here, i.e. they are represented as a product (Remarks: The index n means the index of time on the left hand side
of the equation (discretization) and the power of the amplitude on the right hand side. In order to avoid a
confusion with the complex number i, j was used here as the index of time). The amplitude ξ is also called
amplification factor. If the amplification factor for an iterative equation is larger than 1, it is instable. It is
stable only for values smaller than 1. Calculation of ξ by means of the FCTS scheme is shown in the following
exemplary for the advection equation.

9.2.5.1 Stability of the FTCS scheme for the advection equation

The FTCS scheme for the advection equation (insertion of the right-hand and centered formulas, respectively,
for time and position) reads:

cτ n
a nj +1 = a nj −
2h
(a j +1 − a nj −1 ) .

Insertion of the perturbation yields:

c τ n ikh ( j +1)
ξ n +1 ⋅ e ikhj = ξ n ⋅ e ikhj −
2h
( ξ ⋅e − ξ n ⋅ e ikh ( j −1) )
.
 cτ 
= ξ n ⋅ e ikhj 1 − (e ikh − e −ikh )
 2h 

n ikhj
Dividing both sides by ξ ⋅ e we directly obtain an expression for the amplification factor:

c τ ikh
ξ = 1− (e − e−ikh )
2h
cτ
= 1−i sin (kh )
h
The absolute value of the amplification factor then is:

2
c τ 
ξ = 1 +   sin2 (kh )
h 

This expression is always larger than 1 independent of the time step chosen. By this we have analytically shown
that the scheme is instable without having to iterate the equation numerically. This is made possible by the
skillful choice of the perturbation. Please refer to exercises for the application of the von-Neumann stability
analysis to the Lax-Wendroff scheme.

82
9.3 Boundary value problems: Poisson and Laplace equations
Equations independent of time are described as boundary value problems. The Poisson and Laplace equations
from electrostatics are two important examples. As initial condition, the solution on the boundary of an area
needs to be predetermined. For two dimensions a rectangle with the side lengths Lx and Ly can be chosen as the
boundary13. The x-variable then runs from x=0 to x= Lx and the y-variable from 0 to Ly. The boundary conditions
read

Φ (x = 0, y ) = Φ1 Φ (x = Lx , y ) = Φ2
Φ (x , y = 0) = Φ3 Φ (x , y = Ly ) = Φ4

(Dirichlet’s boundary condition) and

Φ (x = 0, y ) = Φ (x = Lx , y )
Φ (x , y = 0) = Φ (x , y = Ly )

(periodic boundary conditions). The solution is then to be determined in the interior of the rectangle. For the
numerical solution, the x- and y-variables are again discretized by sampling. For this purpose, N and M intervals
are used for the x- and y-dimension, respectively:

Lx
hx = ; x i := i ⋅ hx ; i = 0,1,2, …, N
N
.
L
hy = y ; y j := j ⋅ hy ; j = 0,1,2, …, M
M
Then the solution is calculated for the variable pairs (xi , yj), i.e. the solution points

Φi , j := Φ (x i , x j )

are calculated. Basic methods for solving the Laplace and Poisson equations are investigated in the following.

9.3.1 Laplace equation: Relaxation methods (Jacobi method)

The Laplace equation

∂ 2Φ (x , y ) ∂ 2Φ (x , y )
+ =0
∂x 2 ∂y 2

is identical to the heat equation in two dimensions in its spatial coordinate except for the factor (coefficient of
thermal conduction)

∂T (x , y, t )  ∂ 2T (x , y, t ) ∂ 2T (x , y, t )
= κ ⋅  + 
∂t  ∂x 2 ∂y 2 

For large times (t infinite) the solution of the heat equation approaches a stationary state (stationary temperature
distribution):

lim T (x , y, t ) = Ts (x , y )
t →∞

∂ T (x , y ) ∂ 2Ts (x , y )
2
⇔ s
+ =0
∂x 2 ∂y 2

13
There are other special solutions for cyclindrical and spherical geometries.

83
Thus, the stationary state Ts is at the same time a solution of the Laplace equation. The general idea of the
relaxation methods is to introduce a virtual dependence on time into a time-independent equation and then to
iterate it until a stationary state is reached. The latter is the solution for the stationary, time-independent
equation. In the simplest case this means for the Laplace equation that we write the equation formally like the
heat equation, as shown above, making the potential time-dependent:

∂Φ (x , y, t )  ∂2Φ (x , y, t ) ∂ 2Φ (x , y, t )
= µ ⋅  + 
∂t  ∂x 2 ∂y 2 

µ being an arbitrary convergence factor (corresponding to the coefficient of heat conduction). According to the
FTCS scheme the equation can be discretized as follows:

µτ n µτ n
Φni,+j 1 = Φni, j +
hx 2
(Φi +1, j + Φni−1, j − 2 ⋅ Φni, j ) +
hy 2
(Φi, j +1 + Φni, j −1 − 2 ⋅ Φni, j ) .

The nomenclature for time and position indices is as stated above. Note, however, that the scheme has been
applied to both spatial dimensions here. Hence, the new state at the point (xi,yj) is calculated from the values in
this point and all directly neighbouring points.

Similar to the FTCS scheme of the one-dimensional heat equation, the stability condition for this scheme reads
(without proof):

µτ µτ
+ ≤1 2 .
hx 2 hy 2

Observing the equation in case of identical step sizes hx and hy, i.e.

hx = hy = h ,

the stability condition reads

µτ
≤1 4 .
h2
Since only the stationary state is of interest, the maximal time step is chosen. Insertion of

µτ
=1 4
h2
into the relaxation scheme yields

1 n
Φni ,+j 1 =
4
( Φi +1, j + Φni −1, j + Φni , j +1 + Φni , j −1 ) .

Please note that the central point (xi,yj) does not enter the calculation of the potential at all but only the
neighbouring points because of the maximum time step chosen. This scheme is a special relaxation scheme for
solving the Laplace equation and is called the Jacobi method. Proceeding from the Jacobi method the
convergence velocity can be increased by varying the iteration rules. The schemes resulting from that are called
superrelaxation schemes. Well-known methods are Gauss-Seidel and simultaneous superrelaxation. Since
they are derived from the Jacobi method rather empirically, they will not be treated in detail here (refer to the
textbooks).

84
9.3.1.1 Calculating initial values
The convergence of relaxation methods depends on the initial values chosen to initialize the unknown values in
the interior of the rectangle. If they are chosen skillfully, the system can converge with distinctly less steps. It is
known from the literature that the solution of a Laplace equation can be determined analytically within a
rectangle using the method of separation of variables. Typically, this results in infinite series, which are
exactly defined analytically, but converge slowly in general. Therefore, a numerical solution is still desirable,
although an analytical solution exists. If the infinite series is known, it is the obvious thing to use the first
elements of the series as initial value for the relaxation methods. As an example let us look at the solution of the
Laplace equation on a rectangle with the boundary conditions

Φ (x = 0, y ) = Φ (x = Lx , y ) = Φ (x , y = 0) = 0
Φ (x , y = Ly ) = Φ0

The ansatz of the separation of variables

Φ (x , y ) = X (x ) ⋅Y (y )

with these boundary conditions yields the solution (cf. textbooks)

∞  n πx   n πy L 
4    .
∑
x
Φ (x , y ) = Φ0 ⋅ sin 
 sinh 
N =1,3,5,… πn  L
 x    π y x 
 n L L

If the first element of the series is taken as initial value of the Jacobi method, the convergence velocity increases
considerably as compared to a suboptimal choice of initial values (cf. Exercises).

9.3.2 Poisson equation

The Poisson equation

∂ 2Φ (x , y ) ∂ 2Φ (x , y ) 1
2
+ 2
= ρ (x , y )
∂x ∂y ε0

generally has no analytical solutions for an arbitrary distribution of charges ρ by means of separation of
variables. Therefore, two numerical methods are presented in the following. The charge density is discretized
like the potential:

ρi , j := ρ (x i , x j )

9.3.2.1 Jacobi method

The Jacobi method can be directly extended to solve the Poisson equation. The formalism used in Section 9.3.1
leads to the relaxation equation

1  n 1 
Φni ,+j 1 = Φi +1, j + Φni −1, j + Φni , j +1 + Φni , j −1 + h 2ρi , j  .
4  ε0 

Since there are generally no analytical solutions, the initial values cannot be chosen optimally. Therefore, the
convergence is to be expected to be slow.

85
9.3.2.2 Method of multiple Fourier transform
Contrary to Dirichlet’s boundary conditions, the method of multiple Fourier transform enables the Poisson
equation to be solved on a square with periodic boundary conditions14. The number of samples N and M as well
as the step sizes hx and hy in x- and y-direction must be identical:

N =M
hx = hy = h

Proceeding from that the partial derivatives are discretized with respect to the spatial coordinates by the centered
formula:

∂ 2Φ Φ + Φi −1, j − 2 ⋅ Φi , j
2
→ i +1, j
∂x i , j h2
∂ 2Φ Φ + Φi , j −1 − 2 ⋅ Φi , j
2
→ i , j +1
∂y i , j h2

Insertion into the Poisson equation thus yields:

1 1
2 ( i +1, j
Φ + Φi −1, j + Φi , j +1 + Φi , j −1 − 4 ⋅ Φi , j ) = ρi , j (*)
h ε0

Now we define the two-dimensional Discrete Fourier transform (2-D DFT) of the potential and the charge
density as

N −1 N −1
Fm ,n = ∑ ∑ Φ j ,k ⋅WNjm ⋅WNkn ; m, n = 0, …, N − 1
j =0 k =0
N −1 N −1
Rm ,n = ∑ ∑ ρ j ,k ⋅WNjm ⋅WNkn ; m, n = 0, …, N − 1
j =0 k =0

 2πi 
mit : WN = exp −
 N 

(as compared to the one-dimensional DFT). The equation (*) is now Fourier-transformed on both sides. This
yields:

1 1
2 ( N
W −m + WNm + WN−n + WNn − 4) ⋅ Fm ,n = − Rm ,n
h ε0
(**)
  2
⇔ 2 cos  2πm  + 2 cos  2πn  − 4 ⋅ Fm ,n = − h Rm ,n
  N   N   ε0

using the theorem that:

displacement : Φ'j ,k := Φ j +l ,k
⇒ Fm' ,n = WN−lm ⋅ Fm ,n

(for proof refer to Remark below). Hence, displacement means a multiplication of the spectral components by a
complex factor.

14
Since the discrete Fourier transform considers periodic functions only, it is not surprising that this method
uses periodic boundary conditions.

86
Altogether, the spectral components of the potential can be determined from the (known) spectral components of
the given distribution of the charge by solving the equation (**) for the Fm,n:

Fm ,n = Pm ,n ⋅ Rm ,n
 

2 
h  1 
with : Pm ,n =− 
2ε0  cos  2πm  + cos  2πn  − 2 
  N   N  
and : Φ = IDFT (F )

Thus, the unknown potential is yielded by inverse Fourier transform of the spectral components Fm,n., which, in
turn, result from filtering the spectral components of the charge distribution with the filter P.

Remark: Determination of the spectral components Fm,n from the Rm,n corresponds to a filtering process with the
transfer function Pm,n.

87
Remark: Proof of the displacement theorem:

Φ'j ,k := Φ j +l ,k
N −1 N −1
⇒ Fm' ,n = ∑ ∑ Φ'j ,k ⋅WNjm ⋅WNkn
j =0 k =0
N −1 N −1
⇔ Fm' ,n = ∑ ∑ Φ j +l ,k ⋅WNjm ⋅WNkn
j =0 k =0
N −1+l N −1

∑ ∑Φ ⋅WN( j −l )m ⋅WNkn
'
⇔ Fm' ,n = '
j ,k
j ' =l k =0

j' = j + l
N −1+l N −1 
Fm' ,n = WN−lm ⋅  ∑ ∑ Φ ' ⋅WNj m ⋅WNkn 
'
⇔
 j ' =l k =0 j ,k 
⇔ Fm' ,n = WN−lm ⋅ Fm ,n
Φ periodic with N

88
9.4 Exercises

1. The script dftcs solves the heat equation according to the FTCS scheme for a delta-shaped temperature
distribution across the x-axis as initial value. This distribution “decays” in the course of time by diffusion.
Examine the stability of the solution for different sets of parameters. Use N=41 and several values of the
time step τ within the ranges 10-3 and 10-5. Try several different values of N for τ = 2.0x10-4.
2. Show by insertion that the Gaussian function

1  − ( x − x )2 
T ( x, t ) =
0
exp  
σ (t ) 2π  2σ (t ) 
2
 

with

σ ( t ) = 2κ t

is the solution of the heat equation. Derive a stability rule from the time development of the variance σ
(Hint: Calculate the time required for the variance to rise from 0 to h, h being the grid width of the spatial
coordinate). Does this correspond to the results of Exercise 1?
3. Calculate the time development of the wave function of a free particle (Gaussian wave package) according
to the Crank-Nicholson scheme with periodic boundary conditions. Use the following parameters: L=100;
N=30 and 80; m=1; τ=1; = 1 ; σ0 = 1 and the mean velocity p0 m = k 0 m = 0.5 .
Proceed as follows:
0
- Calculate the initial state Ψ . Ensure that the boundary values are periodic.
- Calculate the matrix of the Hamilton operator. Ensure that it meets the periodic boundary conditions
(which demands are to be made on the matrix in this respect?)
n +1 n
- Calculate the total matrix in order to calculate Ψ from Ψ .
- Iterate the scheme until the wave package has circularly passed the system once
(calculate the number of the steps required for this from L, τ, and the mean velocity).
0 n 2
– Plot Ψ (real and imaginary parts) as well as Ψ for all time points (3D plot).

To what extent does the numerical solution correspond to the analytical solution? Why does the numerical
solution for N=30 not come up to our expectation (cue: sampling theorem)?
4. The script aftcs solves the advection equation according to the FTCS scheme with periodic boundary
conditions and a cosine-modulated Gauss pulse as initial value.

a.) Apply the script with a time step τ=0.002 and N=50 grid points. How do we know that the
obtained numerical solution is not correct?
b.) Solve the equation according to the Lax-Wendroff scheme. Modify the script accordingly and
judge the numerical solution. For this purpose, use time steps just below the maximum stable
time step and the maximum time step itself.

89
5. Apply the von-Neumann stability analysis to the Lax-Wendroff scheme. For this, calculate the formula for
the amount of the amplification factor |ξ | (calculation in writing). Derive the stability condition

h
τ max =
c

by discussing the maximum of |ξ | for the following cases:

< 1

τc 
:= 1
h  
> 1

Hint: Observe the maximum values of the different terms of |ξ | in dependence on the argument x = k⋅h⋅j of
the sine and cosine functions. For illustration, plot |ξ | as a function of x for x at the interval 0 to 2π for these
three cases.
6. The script jacobi solves the Laplace equation on a square. The potential on the margin is zero on three sides
and 1 on one side (y=L). As initial value for the inner points the first member of the infinite series is used
which represents the analytical solution (solution via separation of variables). The time expenditure is
critical. It is basically determined by the number of grid points and the initial values chosen.

a.) Use the script for different spatial solutions (number of grid points N=10 to 30).
How much does the time expenditure increase with the spatial solution? Plot the
number of iterations required for convergence over N and adjust a power law to
the data. Determine the exponent.
b.) Repeat a.), however, use badly chosen initial values (e.g. potential Φ=0 for all
inner points.
7. Solve the Poisson equation

∂ 2 Φ( x , y ) ∂ 2 Φ( x , y ) 1
+ = − ρ( x , y )
∂x 2
∂y 2
ε0

on a square with the side length 1 for the following charge distributions:

a. One point charge of 1 at [x,y]=[0.5 0.5].

b. One point charge of 1 at [x,y]=[0.4 0.6] (What is striking us as compared to a.?)
c. Two point charges of 1 and –1 at [x,y] = [0.5 0.45] and [0.5 0.55](dipole).

For a solution use the method of multiple Fourier Transform (script fftpoi.m). Try to understand the script
as far as possible.

90