AP30002 Computational Physics (2021)

AP30002 Computational Physics (2021)

Assignments / Solutions
Computer Laboratory and Assignment Schedule
Time Mon 4:30 - 6:20 Assignment due date Tutor
Lab 2 11 Oct 18 Oct Anakin
Lab 3 25 Oct 1 Nov Yujia
Lab 4 8 Nov 15 Nov Yujia
Lab 5 22 Nov - Yujia
Test at HJ302 2:30pm 22 Oct Submit on A4 paper Anakin
Lectures will be conducted in hybrid mode: face-to-face at PolyU Campus and broadcast via Blackboard - Collaborate Ultra.
You can also find video recordings of computational physics lectures in 2020 and programming in physics lectures in 2020.
Website: http://apricot.ap.polyu.edu.hk/cp   (login: ep password: ep)
Lecturer: Dr. C.H. Lam (Office: BC615, Tel: 2766-5681)
Tutors (Lab & Grading):
    Mr. Anakin Lee (email: cs.lee@connect.polyu.hk )
    Ms. Yujia Xu (email: yujia1905.xu@connect.polyu.hk   Office: BC521)
Learning Outcomes and Syllabus: AP30002,   Grading rubrics
Assessment: Assignments (20%), Mid-term tests (20%), Exam (60%)
Students can use other general purpose programming languages other than C++ or Python, but please let me know first.
Please report bugs to C.H.Lam@polyu.edu.hk.
Reference books on Computational Physics:
Reference books on C++ Programming Language:
Online C++ Programming Environment Websites: (1) CodingGround, (2) cpp.sh, (3) OnlineGDB
Reference books on Python Programming Language:
Python Websites / Software: (1) OnlineGDB, (2) CodingGround, (3) Rex Tester, (4) Google Colab, (5) Python's website

Contents

1  Introduction
    1.1  What is Computational Physics
    1.2  Examples of Applications
    1.3  How to conduct studies in computational physics?
    1.4  Analytic methods
    1.5  Numerical methods
    1.6  Computer simulation and computer modeling
2  Solution of Nonlinear Equations
    2.1  Fixed point iterations
    2.2  Bisection method
    2.3  Newton-Raphson method
3  Lagrange interpolation
4  Least square fit
    4.1  Least square fit to straight lines
    4.2  Least square fit to polynomials
    4.3  Data with known r.m.s. error
5  Solution of linear equations
    5.1  Gaussian elimination
6  Numerical Integration
    6.1  One-dimensional definite integration
    6.2  Rectangular rule
    6.3  Trapezoidal rule
    6.4  Simpson's rule
7  Numerical differentiation
    7.1  Taylor expansion
    7.2  First order derivative: forward/backward differencing
    7.3  First order derivative: central differencing
    7.4  Second order derivative: central differencing
8  Ordinary differential equation (ODE)
    8.1  First-order ODE
    8.2  Euler Method for 1st-order ODE
    8.3  Second-order ODE
    8.4  Euler Method for 2nd-order ODE
    8.5  System of second-order ODEs
9  ODE: higher order methods
    9.1  Improved Euler method / Euler predictor-corrector method
    9.2  4th-order Runge Kutta (RK4) method
10  ODE: boundary value problem
    10.1  Second-order ODE
    10.2  Shooting method
11  Steady-state heat flow and electric potential
    11.1  Heat Flow in 1D
    11.2  Heated plate and 2D Laplace's Equation
    11.3  Gauss-Seidel relaxation method
    11.4  Electric Potential
12  Time-dependent Heat flow and waves on a string
    12.1  Heated rod (parabolic PDE)
    12.2  Forward Euler's method
    12.3  Waves along an elastic string (hyperbolic PDE)
    12.4  Forward Euler's method

Chapter 1
Introduction

1.1  What is Computational Physics

Approaches in Physics Research:
What is "Computational physics"?
Example of branches of Computational Physics:

1.2  Examples of Applications

    Computational fluid dynamics (CFD): To calculate air flow in order to find better designs which reduce air resistance and increase stability.
    Computational fluid dynamics (CFD): To calculate pressure distribution, air flow (wind), water vapor condensation (cloud and rain formation), and sea water evaporation for weather prediction.
    Computational fluid dynamics (CFD): air flow around the space shuttle during re-entry.
    Computational mechanics: Finite-Element method for calculating the stress field of a structure to identify most stressed regions for improvement.
    Density functional theory (DFT): to calculate electron density in molecules, and electron energy levels.
    Classical molecular dynamics (MD): to understand the flow of nano-meter thick polymer films on a substrate under an applied force. (see movie)

1.3  How to conduct studies in computational physics?

           1. Define the physics problem
         (e.g. weather study)

  − − >  2. understand the basic underlining physics
         (e.g. theories of viscous fluid flow, water vaporization/condensation, heat flow, etc.)

  − − >  3. formulate the underlining physics using mathematical equations
         (e.g. Navier-Stokes equation for fluid flow, gas laws, etc.)

  − − >  4. identify numerical methods to solve the mathematical equations
         (e.g. finite difference method for discretization of differential equations)

  − − >  5. develop (or obtain) computer programs to implement the numerical methods
         (e.g. millions of line of computer codes)

  − − >  6. run the programs to generate numerical results
         (e.g. run in supercomputers)

  − − >  7. analyze and visualize (e.g. plot) numerical results
         (e.g. extract local values of pressure, wind speed, temperature, etc)

  − − >  8. interpret the results to understand the physics
         (e.g. weather forecast and understanding storm formation, global warming, etc)

Steps 2 and 8 are important steps (often neglected by students). They distinguish computational physics from pure numerical analysis (with numerical methods and programming).
Numerical methods and programming are tools for computational physics.
We will spend much time on numerical methods and programming, but the aim is computational physics.

1.4  Analytic methods

Example: Quadratic equation

ax2+bx+c=0
(1.1)
The solution is
x=
−b ±


b2 − 4 a c

2a
(1.2)
We call this an analytic solution, an exact solution, a closed-form solution, or an algebraic solution.
Advantage of analytic solution: a single solution for (nearly) all values of a, b, and c.

1.5  Numerical methods

Example: Analytically non-solvable equation

x=βcos(x)
(1.3)
There is no analytic solution in the form x = f(β) where f is some closed-form (i.e. analytically expressible) function.
However, for any specific value of β, there can be a numerical solution (i.e. solution expressed as a number, rather than an algebra).
e.g. when β = 1, x ≅ 0.73909.
e.g. when β = 0.5, x ≅ 0.45018.
Numerical methods

1.6  Computer simulation and computer modeling

Closely related terminologies:
Computational physics studies in general involve computer modeling, with a purpose to solve a problem in physics.
Computational physics studies usually involve computer simulations (e.g. CFD, MD). Exceptions may be e.g. stress-calculation of bridges and electron density calculation (which are static calculations involving no time-event).
(Other people may have different understandings of these terms.)

Chapter 2
Solution of Nonlinear Equations

2.1  Fixed point iterations

To solve
x = F(x)
(2.1)
Fixed point iterations:
Select an initial guess x0 for the root (can be an arbitrary number, in general).
Let
xn+1 = F( xn )       for n=0,1,2, ....
(2.2)
If x converges,
xn+1=xn=x
(2.3)
Equation (2.2) then gives
x = F(x)
(2.4)
Therefore, x is the root.
In practice, xn for some large n is approximated as the root.
Example:

x = cos(x)
(2.5)
Solution:
Take initial guess x0=1. (Try another number if it does not converge.)
Calculate
xn+1 = cos( xn )       for n=0,1,2, ....
(2.6)
Print out xn.
It is implemented in the following programming for n up to about 100:
fixedpoint.py
# solve x= cos(x) using fixed point interation
from math import *

x = 1  # initial guess
for i in range(100):
    print(i, x)
    x = cos(x)  # iterate
( fixedpoint.py )     
fixedpoint.cpp
// solve x= cos(x) using fixed point interation
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

int main()
{
  cout << setprecision(12);
  double x=1.; // initial guess
  for (int i=0; i<100; ++i){
    cout << i << ' ' << x << endl;
    x = cos(x); // iterate
  }
}
( fixedpoint.cpp )     
Here is the program output. It has converged to x ≅ 0.7390851, which is our numerical approximate of the root.
To observe how xn converges, a plot of xn against n is shown below:
Properties of Fixed point iterations
Example:

x3 − sin(x) = 0
(2.7)
Solution:(outline)
Rearranging the equation, we have
x = sin(x)/x2
(2.8)
Select an initial guess x0. Let
xn+1 = F(xn)      where     F(x)=sin(x)/x2
(2.9)
If it converges, the root is approximated by xn at some large n.

2.2  Bisection method

To solve
F(x)=0
(2.10)
knowing that one root exists in the interval [a1, b1].
Then, it must be true that F(a1)F(b1) < 0.
    First step:
Define midpoint c1 = (a1+b1)/2.
If F(a1)F(c1) < 0, root is in the left half, i.e. [a1,c1].
If F(c1)F(b1) < 0, root is in the right half, i.e. [c1,b1].
The interval to search for the root is halved in each step. After n steps, the root is narrowed down to a small interval [an, bn]. At large n, an (or bn or cn) can be taken as an approximate of the root.
See this java applet for an interactive demonstration.
Algorithm of bisection method (with a,b,c representing an, bn, cn for all n)
    Repeat n times
        c = (a+b)/2
        if F(a) F(c) < 0
            b = c 
        else
            a = c
    Output a

(Algorithm means method of calculation)
Example:
Find a root of
x−cos(x) = 0
(2.11)
in the interval [0,1].
Solution:
bisect.py
# Solve x-cos(x)=0 in [0,1] using bisection method
from math import *

a=0; b=1  # initial range
N=10      # no. of iterations

def F(x): 
    return x-cos(x)

for i in range(N):
    c=(a+b)/2   
    if ( F(a)*F(c)<0 ):
      b=c # take left half of the interval
    else:
      a=c # take right half of the interval
    print (a, b)

( bisection.py )     
bisect.cpp
// Solve x-cos(x)=0 in [0,1] using bisection method
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

double a=0, b=1;  // initial range
const int N=10;   // no. of iterations

double F(double x){ 
  return x-cos(x); 
}

int main(){
  cout << setprecision(12);
  for (int i=0; i<N; ++i){
    double c=(a+b)/2.;  
    if ( F(a)*F(c)<0 ) 
      b=c; // take left half of the interval
    else 
      a=c; // take right half of the interval
    cout << a << "  " << b << endl;
  }
}

( bisection.cpp )     
The output is
0.5  1
0.5  0.75
0.625  0.75
0.6875  0.75
0.71875  0.75
0.734375  0.75
0.734375  0.7421875
0.73828125  0.7421875
0.73828125  0.740234375
0.73828125  0.7392578125

We get x ≅ 0.7383, (or any value in [0.738281, 0.739258]).
The error is 0.739258−0.738281 = 0.000977.
Therefore, x = 0.7383 ±0.0010.
Initial interval:
Error Analysis
The error ∆x for bisection method after n steps is the width of the final interval. Thus,
∆x = bn−cn = bn−1 − cn−1

2
= b1 − c1

2n−1
∝ 2−n
(2.12)
Therefore, ∆x = O( 2−n).     (The O, usually written in caligraphical font , means "of the order".)
Tolerance
One often needs to solve a problem up to within a given error, also called tolerance.
    e.g. Find a root of x−cos(x)=0 up to an error of ϵ = 10−6 (i.e. ∆x < ϵ).
The tolerance may also be given as a fractional error.
    e.g. Find a root of x−cos(x)=0 up to a fractional error of ϵ = 10−6 (i.e. ∆x/x < ϵ).
Example:
Find a root of
x−cos(x) = 0
(2.13)
in the interval [0,1] up to an error of 10−6.
Solution:
bisect_tol.py
# Solve x-cos(x)=0 in [0,1] using bisection method
from math import *

a=0; b=1  # initial range
tol=1e-6  # tolerance

def F(x): 
    return x-cos(x)

while ( b-a > tol ):
    c=(a+b)/2   
    if ( F(a)*F(c)<0 ):
      b=c # take left half of the interval
    else:
      a=c # take right half of the interval
    print (a, b)

print( "x = ", a, " +/- ", b-a )
( bisect_tol.py )     
bisect_tol.cpp
// Solve x-cos(x)=0 in [0,1] within a tolerance using bisection
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

double a=0, b=1;  // initial range
const double tol=1e-6;   // tolerance

double F(double x){ 
  return x-cos(x); 
}

int main(){
  cout << setprecision(10);
  while ( b-a > tol ){
    double c=(a+b)/2.;  
    if ( F(a)*F(c)<0 ) 
      b=c; // take left half of the interval
    else 
      a=c; // take right half of the interval
    cout << a << ' ' << b << endl;
  }
  cout << "x = " << a <<  " +/- " << b-a << endl;
}

( bisect_tol.cpp )     
The output is:
0.5 1
0.5 0.75
0.625 0.75
0.6875 0.75
0.71875 0.75
0.734375 0.75
0.734375 0.7421875
0.73828125 0.7421875
0.73828125 0.740234375
0.73828125 0.7392578125
0.73876953125 0.7392578125
0.739013671875 0.7392578125
0.739013671875 0.7391357421875
0.73907470703125 0.7391357421875
0.73907470703125 0.739105224609375
0.73907470703125 0.7390899658203125
0.7390823364257812 0.7390899658203125
0.7390823364257812 0.7390861511230469
0.7390842437744141 0.7390861511230469
0.7390842437744141 0.7390851974487305
x =  0.7390842437744141  +/-  9.5367431640625e-07

Further extensions:

2.3  Newton-Raphson method

With an estimate of the root xn, a tangent of F(x) at xn is taken as an approximation of F(x). The x−intercept of the tangent, xn+1, is taken as the next (and usually better) estimate of the root.
Since
F′(xn) = F(xn)

xn−xn+1
(2.14)
we arrive at the main equation of Newton-Raphson method:
⇒ xn+1 = xn F(xn)

F′(xn)
(2.15)
The error can be approximated by ∆x = xn+1−xn. (This may slightly underestimate the error.)
See this java applet again (after clicking method  − − >  Newton's method) for an interactive demonstration.
Example:
Find a root of
x−cos(x) = 0
(2.16)
up to an error of 10−6.
Solution:
Using Newton-Rapson method with an initial guess x0=1,
Newton-Raphson.py
# Solve x-cos(x)=0 using Newton-Raphson method
from math import *

tol=1e-6   # tolerance
x=1        # initial guess

def F(x): 
  return x-cos(x)

def Fp(x):  # Derivative of F(x)
  return 1+sin(x)

Dx=tol*10 # error, set a large initial value
n = 1
while Dx > tol:
    xlast = x;           # xlast represents x_{n+1}
    x = x - F(x)/Fp(x)   # x represents x_{n}
    Dx = fabs( x-xlast ) # estimated error
    print (n, x) 
    n = n+1
( Newton-Raphson.py )     
Newton-Raphson.cpp
// Solve x-cos(x)=0 using Newton-Raphson method
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

const double tol=1e-6;   // tolerance
double x=1; // initial guess

double F(double x){ 
  return x-cos(x); 
}

double Fp(double x) { // Derivative of F(x)
  return 1+sin(x); 
}

int main(){
  cout << setprecision(10);
  double Dx=tol*10; // error, set a large initial value

  for (int n=0; Dx > tol; ++n){
    double xlast = x;   // xlast represents x_{n+1}
    x = x - F(x)/Fp(x); // x represents x_{n}
    Dx = fabs( x-xlast ); // estimated error
    cout << n << ' ' << x << endl;
  }
}
( Newton-Raphson.cpp )     
The output is:
1 0.7503638678402439
2 0.7391128909113617
3 0.739085133385284
4 0.7390851332151607

Convergence is much faster than when using the bisection method.
Initial guess x0:
Error Analysis
Let β be the root to be calculated, i.e. F(β)=0.
Then, the error at step n is, ϵn = β− xn.
Newton-Raphson formula is:
xn+1
=
xn F(xn)

F′(xn)
(2.17)
β− xn+1
=
β− xn + F(xn)

F′(xn)
(2.18)
ϵn+1
=
ϵn + F(xn)

F′(xn)
(2.19)
Taylor series expansion about xn:
F(x) = F(xn) + (x−xn) F′(xn) + (x−xn)2

2
F"(xn) + ...
(2.20)
Putting x=β and neglecting higher order terms give
F(β)
F(xn) + (β−xn) F′(xn) + (β−xn)2

2
F"(xn)
(2.21)
⇒ 0
F(xn) + ϵn F′(xn) + ϵn2

2
F"(xn)
(2.22)
⇒ F(xn)
− ϵn F′(xn) − ϵn2

2
F"(xn)
(2.23)
Substituting Eq. (2.23) into Eq. (2.19):
ϵn+1
=
ϵn +
− ϵn F′(xn) − ϵn2

2
F"(xn)

F′(xn)
(2.24)
=
ϵn2

2
F"(xn)

F′(xn)
(2.25)
Close to convergence, F′(xn) ≅ F′(β) and F"(xn) ≅ F"(β). Then, after taking absolute values,
⇒ | ϵn+1 |
| ϵn2 |

2
| F"(β)

F′(β)
|
(2.26)
This is called quadratic convergence: the number of significant figures approximately double in each step. This is very fast convergence.
Issue

Chapter 3
Lagrange interpolation

Interpolation: To find a smooth interpolation function p(x) going exactly through m+1 points (x0,y0), (x1,y1), ... (xm,ym).
Lagrange interpolation: To construct an interpolation function which is a m-th order polynomial using Lagrange polynomials.
Example:
Interpolate the data points (x,y) = (0, 2), (1, 8), (1.5, 12), (2, 10) and (3, 20). Furthermore, estimate y when x=1.8.
Solution:
At x=1.8, y ≅ 11.3.
(method to be explained below)
See this java applet for an interactive demonstration.
Two data points (m=1)
Two points can be interpolated by a straight line (i.e. linear interpolation):
p1(x) = a1 x + a0
(3.1)
Since p1(x) should go through (x0,y0) and (x1,y1), we have
y0 = p(x0) = a1 x0 + a0
(3.2)
y1 = p(x1) = a1 x1 + a0
(3.3)
Solving a1 and a0 (or simple inspection) gives
p1(x) = y0 x−x1

x0−x1
+ y1 x−x0

x1−x0
(3.4)
It can be rewritten as
p1(x) = y0 L10(x) + y1 L11(x)
(3.5)
where the Lagrange polynomials of degree 1 is defined by
L10(x) = x−x1

x0−x1
        and             L11(x) = x−x0

x1−x0
(3.6)
The Lagrange polynomials satisfy
L10(x0)
= 1         and             L11(x0)
= 0
(3.7)
L10(x1)
= 0         and             L11(x1)
= 1
(3.8)
Three data points (m=2)
The Lagrange interpolation function is
p2(x) = y0 L20(x) + y1 L21(x) + y2 L22(x)
(3.9)
where the Lagrange polynomials of degree 2 is defined by
L20(x) = (x−x1)(x−x2)

(x0−x1)(x0−x2)
 ,    L21(x) = (x−x0)(x−x2)

(x1−x0)(x1−x2)
 ,    L22(x) = (x−x0)(x−x1)

(x2−x0)(x2−x1)
(3.10)
The Lagrange polynomials satisfy
L20(x0)
=1  ,    L21(x0)
=0  ,    L22(x0)=0
(3.11)
L20(x1)
=0  ,    L21(x1)
=1  ,    L22(x1)=0
(3.12)
L20(x2)
=0  ,    L21(x2)
=0  ,    L22(x2)=1
(3.13)
Many data points (m ≥ 1)
The Lagrange interpolation function is
pm(x) = y0 Lm0(x) + y1 Lm1(x) + ... + ym Lmm(x)
(3.14)
where the Lagrange polynomials of degree m is defined by
Lmj(x)
=
(x−x0) ... (x−xj−1)   (x−xj+1) ... (x−xm)

(xj−x0) ... (xj−xj−1)  (xj−xj+1) ... (xj−xm)
(3.15)
=
m

k=0, k ≠ j 
x − xk

xj −xk
(3.16)
Note that (x−xj) is omitted from the numerator while (xj−xj) is omitted from the denominator.
It satisfies
Lmj(xi) = {
1
for j = i
0
for ji
(3.17)
Example: (continued)
Lagrange.py
# Lagrange interpolation of m+1 data points
import numpy as np
import matplotlib.pyplot as plt

m=4   #  number of data points - 1 
X = np.array( [0., 1., 1.5, 2,   3 ] )  # x co-ord's 
Y = np.array( [2., 8., 12., 10., 20] )  # y co-ord's 

def p(x):  # Lagrange interpolation function
    px = 0; 
    for j in range(m+1):
        # calculate Lagrange polymomial L_mj(x) 
        Lmj=1   
        for k in range(m+1):
            if k!=j: 
                # multiply k-th factor to Lmj
                Lmj *= (x-X[k])/(X[j]-X[k])
        px += Y[j] * Lmj; # add j-th term to px
    return px

# plot data points and interpolation function
plt.scatter(X, Y)
Xinterp = np.linspace(min(X), max(X), 100)
plt.plot(Xinterp, p(Xinterp))
( Lagrange.py )     
Lagrange.cpp
// Lagrange interpolation of m+1 data points
#include <iostream>
#include <iomanip>
#include <math.h>
#include <fstream>
using namespace std;

const int m=4; // number of data points - 1 
const double X[m+1] = { 0., 1., 1.5, 2.,  3. };  // x co-ord's 
const double Y[m+1] = { 2., 8., 12., 10., 20.};  // y co-ord's 

double p(double x) // Lagrange interpolation function
{
  double px=0; 
  for (int j=0; j<=m; ++j){
    // calculate Lagrange polymomial L_mj(x) 
    double Lmj=1.;   
    for (int k=0; k<=m; ++k){
      if (k!=j) 
        // multiply k-th factor to Lmj
        Lmj *= (x-X[k])/(X[j]-X[k]); 
    }
    px += Y[j] * Lmj; // add j-th term to px
  }
  return px;              
}

int main()                    
{
  // output data points to file
  ofstream datafile("Lagrange_data.txt");
  for (int j=0; j<=m; ++j)
    datafile << X[j] << ' ' << Y[j] << endl;
  datafile.close();

  // calculate and output 100 interpolated points to file
  ofstream interfile("Lagrange_inter.txt");
  for (double x=X[0]; x<=X[m]; x += (X[m]-X[0])/100.)  
    interfile << x << "  " << p(x) << endl;   
  interfile.close();
}

( Lagrange.cpp )     
data files: Lagrange_data.txt,    Lagrange_inter.txt
Issues:

Chapter 4
Least square fit

4.1  Least square fit to straight lines

Given N data points (x1, y1), (x2, y2), ... (xN, yN) and a fitted function which is a straight line
p(x)=c0 + c1 x
(4.1)
find the fitting parameters c0 and c1 which gives the best fit.
See this java applet for an interactive demonstration.

In least square fit, best fit occurs when the sum of squared errors S(c0,c1) is minimized, where
S(c0,c1)
=
N

j=1 
( p(xj)−yj )2 = Σ ( c0 + c1 xj −yj )2
(4.2)
Minimizing S(c0,c1) w.r.t. c0 and c1,
∂S(c0,c1)

∂c0
=
Σ 2 (c0+c1 xj −yj ) = 0
(4.3)
∂S(c0,c1)

∂c1
=
Σ 2 xj (c0+c1 xj −yj )=0
(4.4)
Simplification gives
c0 N
+
c1Σ xj
Σ yj
=
0
(4.5)
c0 Σ xj
+
c1 Σ xj2
Σ xj yj
=
0
(4.6)
Solution gives the fitted parameters:
c0
=
Σ yjΣ xj2 −Σ xj Σ xj yj

N Σ xj2 − ( Σ xj)2
(4.7)
c1
=
N Σ xj yj −Σ yj Σ xj

N Σ xj2 − ( Σ xj)2
(4.8)

4.2  Least square fit to polynomials

Fitting N data points using a polynomial as a fitting function, i.e.
p(x)=c0 + c1 x + c2 x2 + ... + cm xm
(4.9)
Sum of squared errors becomes
S(c0,c1, ... cm)
=
N

j=1 
( p(xj)−yj )2
(4.10)
=
Σ ( c0 + c1 xj + c2 xj2 + ..... + cm xjm −yj )2
(4.11)
Minimizing S(c0,c1, ... cm) w.r.t. ci  (i=0,1, .. m) implies ∂S(c0, c1, ... cm)/∂ci=0,
which implies the i-th equation in
c0 N
+
c1Σ xj
+
c2 Σ xj2
+ ..... +
cm Σ xjm
Σ yj
=
0
(4.12)
c0 Σ xj
+
c1Σ xj2
+
c2 Σ xj3
+ ..... +
cm Σ xjm+1
Σ xj yj
=
0
(4.13)
......
(4.14)
c0 Σ xjm
+
c1Σ xjm+1
+
c2 Σ xjm+2
+ ..... +
cm Σ xjm+m
Σ xjm yj
=
0
(4.15)
i.e. m+1 linear equations with m+1 unknowns.
Solution by Gaussian elimination (to be discussed later) gives c0... cm.

4.3  Data with known r.m.s. error

To fit data points
(x1, y1±σ1), (x2, y2±σ2), ..... (xN, yN±σm)
where σj is the r.m.s. error (i.e. standard deviate) of yj. Least square fit is performed by minimizing the weighted sum of squared errors:
S(c0,c1,..,cm)
=
N

j=1 

p(xj)−yj

σi

2

 
(4.16)

Chapter 5
Solution of linear equations

5.1  Gaussian elimination

To solve a system of n linear equations with n unknowns (x1,x2,... , xn),
A11 x1 + A12 x2 + ..... + A1n xn
= b1
A21 x1 + A22 x2 + ..... + A2n xn
= b2
...                
An1 x1 + An2 x2 + ..... + Ann xn
= bn
Using matrix representation:






A11
A12
.....
A1n
A21
A22
.....
A2n
...
...
...
An1
An2
.....
Ann












x1
x2
...
xn






=





b1
b2
...
bn






i.e.
A x = b
where A is a n×n matrix and x and b are column vectors.

Gaussian elimination for n=3
A11 x1 + A12 x2 + A13 x3
= b1
(5.1)
A21 x1 + A22 x2 + A23 x3
= b2
(5.2)
A31 x1 + A32 x2 + A33 x3
= b3
(5.3)
Step A: Construction of upper triangular system
(A1) Eliminating x1 from Eqs (5.2) and (5.3)
(5.1) : 
   A11 x1
+
A12 x2
+
A13 x3
=
b1
(5.4)
(5.2)− A21

A11
×(5.1) : 
    0
+
(A22 A21

A11
A12) x2
+
(A23 A21

A11
A13) x3
=
b2 A21

A11
b1
(5.5)
(5.3)− A31

A11
×(5.1) : 
    0
+
(A32 A31

A11
A12) x2
+
(A33 A31

A11
A13) x3
=
b3 A31

A11
b1
(5.6)
and rewrite as
A11 x1 + A12 x2 + A13 x3
= b1
(5.7)
0 + A22′x2 + A23′x3
= b2
(5.8)
0 + A32′x2 + A33′x3
= b3
(5.9)
(A2) Eliminate x2 from Eq (5.9)
            (5.7)⇒ 
A11 x1
+
A12 x2
+
A13 x3
=
b1             
(5.10)
(5.8)⇒ 
0
+
A22′x2
+
A23′x3
=
b2′            
(5.11)
(5.9)− A32

A22
×(5.8)⇒ 
0
+
0
+
(A33′− A32

A22
A23′) x3
=
b3′− A32

A22
b2′            
(5.12)
Rewrite as
            
A11 x1
+
A12 x2
+
A13 x3
=
b1                   
(5.13)
0
+
A22′x2
+
A23′x3
=
b2
(5.14)
0
+
0
+
A33"x3
=
b3"
(5.15)
The non-zero terms on the lhs follows an upper triangular form.
Solving Eq. (5.15)   gives x3
Step B: Back-substitution
(B1) Solving Eq. (5.14) using known x3   gives x2
(B2) Solving Eq. (5.13) using known x2 & x3   gives x1

Note that
Example:
Solve
−3  x1
+
2 x2
6 x3
=
6
(5.16)
5  x1
+
7 x2
5 x3
=
6
(5.17)
1  x1
+
4 x2
2 x3
=
8
(5.18)
Solution:
GaussianElimination.py
# Solving linear equations Ax=b by Gaussian elimination
# without checking for zero-coefficients
import numpy as np

n=3
A = np.array( [[ -3., 2., -6. ],
               [  5., 7., -5. ],
               [  1., 4., -2. ] ] ) 
b = np.array( [ 6., 
                6., 
                8. ] ) 
x = np.array( [0]*3 ) # answer: -2, 3, 1

# Construction of upper triangular system
for k in range(n-1): #  substep A-k
    for i in range(k+1, n): # i-th equation
        for j in range(k+1, n): # j-th term
            A[i][j] = A[i][j] - (A[i][k]/A[k][k])*A[k][j]  # Aij -> Aij' -> Aij'' -> ...
        b[i] = b[i] - (A[i][k]/A[k][k])*b[k]  # Bi -> Bi' -> Bi'' -> ... 
x[n-1]= b[n-1]/A[n-1][n-1]  # solution of last equation of upper triangle

# back-substitution 
for k in range(n-1): # substep B-k
    i=n-k-2  # i-th equation
    for j in range(i+1, n):  # j-th term
        b[i] -= A[i][j] * x[j]  # subtract j-th term in i-th equation
    x[i] = b[i] / A[i][i]

# output results
print(x)
( GaussianElimination.py )     
GaussianElimination.cpp
// Solving linear equations Ax=b by Gaussian elimination
// without checking for zero-coefficients
#include <iostream>
#include <math.h>
using namespace std;

const int n=3;
double A[n][n] = { { -3, 2, -6 },
                   {  5, 7, -5 },
                   {  1, 4, -2 } };
double b[n] = { 6, 
                6, 
                8 };
double x[n]; // answer: -2, 3, 1

int main()
{
  // Construction of upper triangular system
  for (int k=0; k<n-1; ++k) { // substep A-k
    for (int i=k+1; i<n; ++i) { // i-th equation
      for (int j=k+1; j<n; ++j) { // j-th term
        A[i][j] = A[i][j] - (A[i][k]/A[k][k])*A[k][j]; // Aij -> Aij' -> Aij'' -> ...
      }
      b[i] = b[i] - (A[i][k]/A[k][k])*b[k]; // Bi -> Bi' -> Bi'' -> ... 
    }
  }
  x[n-1]= b[n-1]/A[n-1][n-1]; // solution of last equation of upper triangle

  /* back-substitution */
  for (int k=0; k<n-1; ++k) { // substep B-k
    int i=n-k-2; // i-th equation
    for (int j=i+1; j<n; ++j) // j-th term
      b[i] -= A[i][j] * x[j]; // subtract j-th term in i-th equation
    x[i] = b[i] / A[i][i];
  }

  // output results
  for (int i=0; i<n; ++i) cout << x[i] << endl;
}
( GaussianElimination.cpp )     

Chapter 6
Numerical Integration

6.1  One-dimensional definite integration

To find the integral numerically
I =
b

a 
f(x) dx
(6.1)
Possible difficulties in analytic approach:

6.2  Rectangular rule

I equals the area under the curve within [a,b]. It is approximated by the sum of the areas of the rectangles.
See this java applet (which takes the right-end-point approach) for an interactive demonstration.

Divide [a,b] into n equal intervals of width h given by
h = b−a

n
(6.2)
The intervals are bounded by n+1 points
xi = a + i h       for i=0,1, ... n
(6.3)
so that x0=a and xn=b.
Area of the i−th rectangle is hf(xi), taking the left-end point approach. By summing up areas of all n rectangles, we get
I ≅ n−1

i=0 
hf(xi)
(6.4)
which gives the rectangular rule
I ≅ h n−1

i=0 
f(xi)
(6.5)

Error analysis
Local error: error in one interval
For 1st interval for example, Taylor expansion of f(x) about x0:
f(x)
=
f(x0) + (x−x0) f′(x0) + ...
(6.6)

x1

x0 
f(x) dx
=

x0+h

x0 
f(x0) + (x−x0) f′(x0) dx + ...
(6.7)
=
h f(x0)          + h2

2
f′(x0) + ...
(6.8)
=
approximate + local error
(6.9)
Therefore,
Local error = h2

2
f′(x0) +  ...  = O(h2)
(6.10)
It is similar for other intervals. Thus, the local error is of order h2.
Global error: error in I  (i.e. the sum over all intervals)
Eq. (6.4) implies
Global error
=
n−1

i=0 
local error at i−th interval
(6.11)
=
n−1

i=0 
O(h2) = n O(h2) = b−a

h
O(h2) = O(h)
(6.12)
Thus, the global error is of order h.

6.3  Trapezoidal rule

I is approximated by the sum of the areas of the trapezoids, i.e.
I
n−1

i=0 
h [ f(xi)+ f(xi+1) ]

2
(6.13)
=
h

2
[ f(x0) + 2f(x1) + 2f(x2) + ... + 2f(xn−1) + f(xn) ]
(6.14)
=
h [ 1

2
f(x0) + f(x1) + f(x2) + ... + f(xn−1) + 1

2
f(xn) ]
(6.15)
Using Taylor's expansion up to the f"(x0) term, i.e.
f(x) = f(x0) + (x−x0) f′(x0) + (x−x0)2

2
f"(x0) + ...
(6.16)
it can be similarly shown that local error = O(h3) and global error = O(h2).

6.4  Simpson's rule

Intervals are paired up. At the bi-interval [ xi, xi+2]   (i is even), the curve is approximated by a parabola
p(x) = a(x−xi)2 + b(x−xi) + c
(6.17)
For best approximation, we require that f(x)=p(x) at xi, xi+1 and xi+2, i.e.
f(xi)
=
p(xi)
=
c
(6.18)
f(xi+1)
=
p(xi+1)
=
p(xi+h)
=
ah2 + bh+c
(6.19)
f(xi+2)
=
p(xi+2)
=
p(xi+2h)
=
4ah2 + 2bh+c
(6.20)
Then,

xi+2

xi 
f(x) dx

xi+2

xi 
p(x) dx
(6.21)
=
h

3
[ 8h2a+6hb+6c ]                (using Eq. (6.17))
(6.22)
=
h

3
[ c + 4(ah2+bh+c) + (4ah2 + 2bh+c) ]
(6.23)
=
h

3
[ f(xi) + 4f(xi+1) + f(xi+2) ]            (using Eqs. (6.186.20))
(6.24)
Summing up all terms,
I
=
n−2

i=0,2, ... 

xi+2

xi 
f(x) dx
(6.25)
n−2

i=0,2, ... 
h

3
[ f(xi) + 4f(xi+1) + f(xi+2) ]
(6.26)
=
h

3
[ f(x0) + 4f(x1) + 2f(x2) + 4f(x3) + 2f(x4) +   ...   + 4f(xn−1) + f(xn) ]
(6.27)
Using Taylor's expansion up to the f(4)(x0) term, it can be similarly shown that local error = O(h5) and global error = O(h4).

Chapter 7
Numerical differentiation

7.1  Taylor expansion

Taylor's expansion of f(x) about point x
f(x+h) = f(x) + hf′(x) + h2

2!
f"(x) + ... + hn

n!
fn(x) + Rn
(7.1)
where Rn is called the n-th remainder and is given by
Rn
=
hn+1

(n+1)!
fn+1(ξ),         ξ ∈ (x,x+h)
(7.2)
=
O(hn+1)
(7.3)
Therefore, the Taylor's expansion can also be written as
f(x+h) = f(x) + hf′(x) + h2

2!
f"(x) + ... + hn

n!
fn(x) + O(hn+1)
(7.4)
i.e. the error term is of order hn+1
In 2D, Taylor's expansion is
f(x+h, y+k)
=
f(x,y) + (h

∂x
+ k

∂y
) f|x,y + 1

2!
(h

∂x
+ k

∂y
)2 f|x,y     + ... +
(7.5)
+ 1

n!
(h

∂x
+ k

∂y
)n f|x,y + 1

(n+1)!
(h

∂x
+ k

∂y
)n+1 f|x+αh,y+βk
(7.6)
where α, β ∈ (0,1).
Assuming k ∝ h,
f(x+h, y+k)
=
f(x,y) + (h

∂x
+ k

∂y
) f|x,y + 1

2!
(h

∂x
+ k

∂y
)2 f|x,y     + ... +
(7.7)
+ 1

n!
(h

∂x
+ k

∂y
)n f|x,y + O(hn+1)
(7.8)
For n=1, we get
f(x+h, y+k)
=
f(x,y) + h fx(x,y) + k fy(x,y) + O(h2)
(7.9)

7.2  First order derivative: forward/backward differencing

A derivative is defined based on the infinitesimal limit of h as
f′(x) =
lim
h→ 0 
f(x+h) − f(x)

h
(7.10)
Now, assuming that h is a finite value without taking the infinitesimal limit, we get the finite difference approximation of the derivative
f′(x) ≅ f(x+h) − f(x)

h
(7.11)
Usually, h > 0 is assumed. The above formula is further classified as a forward difference formula.
Similarly, a backward difference formula is
f′(x) ≅ f(x) − f(x−h)

h
(7.12)
Error Analysis
Taylor expansion (Eq. (7.4)) gives
f(x+h)
=
f(x) + hf′(x) + O(h2)
(7.13)
f′(x)
=
f(x+h) − f(x)

h
O(h2)

h
(7.14)
=
f(x+h) − f(x)

h
+ O(h)
(7.15)
It is similar for backward difference.
Therefore, forward/backward difference formula has an error O(h).

7.3  First order derivative: central differencing

Taking the mean of Eqs (7.11) and (7.12), we get the central difference formula
f′(x) ≅ f(x+h) − f(x−h)

2h
(7.16)
Error Analysis
Taylor's expansions (with +h and -h) :
f(x+h)
=
f(x) + hf′(x) + h2

2!
f"(x) + O(h3)
(7.17)
f(x−h)
=
f(x) − hf′(x) + h2

2!
f"(x) + O(h3)
(7.18)
Subtracting the two expansions gives
f(x+h) − f(x−h)
=
2 hf′(x) + O(h3)
(7.19)
⇒ f′(x)
=
f(x+h) − f(x−h)

2h
+ O(h3)

2h
(7.20)
=
f(x+h) − f(x−h)

2h
+ O(h2)
(7.21)
Therefore, central difference formula has an error O(h2).
Note that

7.4  Second order derivative: central differencing

Central difference formula:

f"(x) ≅ f(x+h) −2f(x) + f(x−h)

h2
(7.22)
Error Analysis
Taylor's expansions (with +h and -h) :
f(x+h)
=
f(x) + hf′(x) + h2

2!
f"(x) + h3

3!
f"′(x) + O(h4)
(7.23)
f(x−h)
=
f(x) − hf′(x) + h2

2!
f"(x) − h3

3!
f"′(x) + O(h4)
(7.24)
Adding the two expansions gives
f(x+h) + f(x−h)
=
2 f(x) + h2f"(x) + O(h4)
(7.25)
⇒ f"(x)
=
f(x+h) −2f(x) + f(x−h)

h2
+ O(h4)

h2
(7.26)
=
f(x+h) −2f(x) + f(x−h)

h2
+ O(h2)
(7.27)
Therefore, central difference formula has an error O(h2).

Chapter 8
Ordinary differential equation (ODE)

8.1  First-order ODE

Example: Falling sphere under viscous drag
    Sphere under gravity Fg and viscous drag Fd.
Equation of motion:
Mx" = Mg − λx′
(8.1)
where M=mass and λ = viscosity coefficient.
Velocity, v=x′. Thus,
Mv′ = Mg − λv
(8.2)
Analytic solution:
v = M

λ
( g − exp( − λ

M
(t+c) )
(8.3)
where c is an integration constant. Assuming an initial condition v(0)=0,
v = Mg

λ

1 − exp( − λt

M
)
(8.4)
Example: Falling sphere under viscous drag and external force
Equation of motion:
Mv′(t) = Mg − λv(t) +F(t)
(8.5)
where F(t) is a time-dependent external force, e.g. F(t) = cos(ωt).
In general, there is no analytic solution for v(t).
Example: Classical particle in 1D potential
A particle of mass M and total energy E moving in the +x direction has a potential energy V(x) at position x. The velocity x′ is given by
x′(t) =

 

2 (E−V( x(t) )) /M
 
(8.6)
Find x(t).
In general, consider the initial value problem of 1st-order ODE
x′(t)
=
f( t, x(t) )       
differential equation (DE)
(8.7)
x(t0)
=
x0
initial condition (I.C.)
(8.8)
The aim is to find x(t) for t ≥ t0 numerically.
Note that

8.2  Euler Method for 1st-order ODE

One step
Taylor's expansion:
x(t0+∆t) = x(t0) + ∆t  x′(t0) + O(∆t2)
(8.9)
It gives the Euler method (also called the forward Euler method):
x(t0+∆t) ≅ x(t0) + ∆t  f( t0, x(t0) )
(8.10)
which has a local error O(∆t2). Note that all values on the R.H.S. are known or readily computable.
N steps
To calculate x( t0+T ), divide the interval T into N steps of size ∆t = T/N. Let
ti
=
t0 + i ∆t
(8.11)
xi
=
x(ti)
(8.12)
and hence tN=t0+T.
Similar to Eq. (8.10), Euler method at the i-th iteration step can be written as
xi+1 = xi + ∆t  f( ti, xi )
(8.13)
Algorithm:
       Using x0 (I.C.) and applying Eq. (8.13) with i=0, we get x1.
       Then, using x1 and applying Eq. (8.13) with i=1, we get x2.
       Then, using x2 and applying Eq. (8.13) with i=2, we get x3.
                       ...
       Then, using xi and applying Eq. (8.13) with i ≥ 0, we get xi+1.
Error:
       Calculating x(t0+T) requires N = T/∆t steps.
       Therefore, global error ≅ N × local error = (T/∆t) ×O(∆t2) = O(∆t).
Example: Classical particle in 1D potential (revisit)
A particle of mass M=1 g and total energy E=1.001 J moving in the +x direction has a potential energy V(x) = c/x at position x, where c=0.1 Jm. The velocity x′ is given by
x′(t) =

 

2 (E−V( x(t) )) /M
 
(8.14)
At t=0, x=0.1 m. Find x(t) for 0 ≤ t ≤ 0.01 s.
Solution: Putting
f(t,x) =

 

2 (E−V( x )) /M
 
(8.15)
Euler-ptcle.py
# Euler method for particle in a 1D potential 
from math import *
import numpy as np

M = 1e-3      # mass in kg
E = 1.001     # total energy in J
c = 0.1       # potential energy parameter in Jm
T = 1e-2      # final time in s
dt= 0.001     # time step of Euler method
x = 0.1       # position in m   

def f(t, x): 
    V = c/x                       # 1D potential
    return sqrt( 2 * (E - V) /M ) # velocity

for t in np.arange(0, T, dt):
    print(t, x)            # output t, x(t)
    x = x + dt * f(t, x)    # Euler method iteration     

( Euler-ptcle.py )     
Euler-ptcle.cpp
// Euler method for particle in a 1D potential 
#include <iostream>
#include <math.h>
using namespace std;

const double M = 1e-3;      // mass in kg
const double E = 1.001;     // total energy in J
const double c = 0.1;       // potential energy parameter in Jm
const double T = 1e-2;      // final time in s
const double dt= 0.001;     // time step of Euler method
double x = 0.1;             // position in m

double f(double t, double x){ 
  double V = c/x;                // 1D potential
  return sqrt( 2 * (E - V) /M ); // velocity
}

int main(){
  for (double t=0; t<=T; t+=dt){
    cout << t << " " << x << endl;   // output t, x(t)
    x = x + dt * f(t, x);            // Euler method iteration     
  }
}

( Euler-ptcle.cpp )     
Plot of x(t) against t for ∆t = 0.001, 0.0001, and 0.00001:
i.e. particle is under repulsion from the origin with a small initial velocity.

8.3  Second-order ODE

Example: Forced and damped harmonic oscillator
    Equation of motion:
Mx" = − k x + Mg − λx′ + F(t)
(8.16)
where
        x=extension of spring,
        k=spring constant,
        M=mass,
        λ = viscosity coefficient,
        F(t)=external force.
With the −kx term, it cannot be reduced to a 1st-order ODE
In general, consider the initial value problem of 2nd-order ODE
x"(t)
=
f( t, x(t), x′(t) )       
differential equation (DE)
(8.17)
x(t0)
=
x0
initial conditions (I.C.)
(8.18)
x′(t0)
=
v0
(8.19)
The aim is to find x(t) for t ≥ t0 numerically.

8.4  Euler Method for 2nd-order ODE

Introducing a new variable v(t)=x′(t), Eq. (8.17) is transformed into a system of 1st-order ODEs:
x′(t)
=
v(t)       
differential equation (DE)
(8.20)
v′(t)
=
f( t, x(t), v(t) )
(8.21)
x(t0)
=
x0
initial condition (I.C.)
(8.22)
v(t0)
=
v0
(8.23)
To calculate x( t0+T ), divide the interval T into N steps of size ∆t = T/N. Let
ti
=
t0 + i ∆t
(8.24)
xi
=
x(ti)
(8.25)
vi
=
v(ti)
(8.26)
The Euler method at the i-th iteration step is
xi+1
=
xi + ∆t  vi
(8.27)
vi+1
=
vi + ∆t  f( ti, xi, vi )
(8.28)
Avoid a common mistake in programming which implements the above equations incorrectly as
xi+1
=
xi + ∆t  vi
(8.29)
vi+1
=
vi + ∆t  f( ti, xi+1, vi )
(8.30)
Example: Forced and damped harmonic oscillator (continued)
Now, assume that M=0.1 kg, k=10 N/m, λ=0.7 Ns/m, and F(t)=0. Initially, x=0 and v=0. Find x(t).
Solution: The equation of motion can be rewritten as :
x" = f( t, x, v) = − k

M
x + g − λ

M
v
(8.31)
Euler-oscillator.py
# Euler method for damped oscillator 
import numpy as np

M = 0.1      # mass
k = 10.      # spring constant
g = 9.8      # gravity
Lambda = 0.7 # damping coefficient
T = 1.       # final time
dt= 0.01     # time step of Euler method

x = 0.       # position  m
v = 0.       # velocity in m/s
eps=1e-12    # tolerance for machine error

def f(t, x, v): 
    return -k/M*x +g - Lambda/M*v   # particle acceleration

for t in np.arange(0, T+eps, dt):
    print( t, x)             # output t, x(t)
    dx = dt * v              # Euler method iteration     
    dv = dt * f(t, x, v)
    x += dx
    v += dv
( Euler-oscillator.py )     
Euler-oscillator.cpp
// Euler method for damped oscillator 
#include <iostream>
#include <math.h>
using namespace std;

const double M = 0.1;      // mass
const double k = 10;       // spring constant
const double g = 9.8;      // gravity
const double lambda = 0.7; // damping coefficient
const double T = 1.;       // final time
const double dt= 0.01;     // time step of Euler method
double x = 0.;             // position 
double v = 0.;             // velocity 
const double eps=1e-12;    // tolerance for machine error

double f(double t, double x, double v){ 
  return -k/M*x +g - lambda/M*v;   // particle acceleration
}

int main(){
  for (double t=0; t<=T+eps; t+=dt){
    cout << t << " " << x << endl;   // output t, x(t)
    double dx = dt * v;              // Euler method iteration     
    double dv = dt * f(t, x, v);
    x += dx;
    v += dv;
  }
}

( Euler-oscillator.cpp )     

8.5  System of second-order ODEs

Consider N particles in 3D interacting with each other.
Their dynamics follow a system of 3N 2nd-order ODE for the variables xk, yk, zk      (k=1 .. N).
They can be transformed into 6N 1st-order ODE for the variables xk, yk, zk, vxk, vyk, vzk.
They can then be solved by Euler method.

Chapter 9
ODE: higher order methods

9.1  Improved Euler method / Euler predictor-corrector method

Consider the 1st order ODE
x′(t)
=
f( t, x(t) )
(9.1)
    Prediction based on average slope is more accurate.
First, obtain a predicted value by applying the Euler method (which is a forward difference method)
xpi+1
=
xi + ∆t  x′(ti)
(9.2)
=
xi + ∆t  f( ti, xi )
(9.3)
Second, obtain a more accurate corrected value, using a central difference method:
x(ti+∆t)
xi + ∆t  x′(ti+∆t/2)
(9.4)
xi + ∆t  x′(ti) + x′(ti+1)

2
(9.5)
xi + 1

2
∆t  ( f( ti, xi ) + f( ti+1, xpi+1 ) )
(9.6)
Therefore, the improved Euler method is a two-step method defined by
predicted value:
       xpi+1
=
xi + ∆t  f( ti, xi )
(9.7)
corrected value:
       xi+1
=
xi + 1

2
∆t  ( f( ti, xi ) + f( ti+1, xpi+1 ) )
(9.8)
Error:
Denote exact value by x(t0+∆t) and estimates by xpi+1 and xi+1.
From error analysis of Euler method, error of predictor is given by
x(ti+1) = xpi+1 + O(∆t2)
(9.9)
Then,
x′(ti+1) = f( ti+1, x(ti+1) )
=
f( ti+1, xPi+1 + O(∆t2) )
(9.10)
=
f( ti+1, xPi+1 )  +  ∂f/∂x   O(∆t2)           (Using Taylor′s expansion)
(9.11)
=
f( ti+1, xPi+1 )  +  O(∆t2)
(9.12)
Taylor's expansion:
x(ti+1)
=
x(ti) + ∆t  x′(ti) + 1

2
∆t2 x"(ti) + O(∆t3)
(9.13)
=
xi + ∆t  x′(ti) + 1

2
∆t2 ( x′(ti+1) − x′(ti)

∆t
+ O(∆t) ) + O(∆t3)       (forward differencing)
(9.14)
=
xi + 1

2
∆t  ( x′(ti) + x′(ti+1) ) + O(∆t3)
(9.15)
=
xi + 1

2
∆t  ( f( ti, xi ) + f( ti+1, xPi+1 )  +  O(∆t2) ) + O(∆t3)
(9.16)
=
xi + 1

2
∆t  ( f( ti, xi ) + f( ti+1, xPi+1 ) ) + O(∆t3)
(9.17)
=
xi+1 + O(∆t3)
(9.18)
i.e. local error = O(∆t3), and global error = O(∆t2).
 
Note: Improved Euler method can be generalized to 2nd order ODEs, and systems of ODEs.

9.2  4th-order Runge Kutta (RK4) method

Most widely used method for ODEs.
xi+1 = xi + ∆t

6
( f0 + 2f1 + 2f2 + f3)
(9.19)
where
f0
=
f(ti, xi)
(9.20)
f1
=
f(ti + ∆t

2
, xi + ∆t

2
f0 )
(9.21)
f2
=
f(ti + ∆t

2
, xi + ∆t

2
f1 )
(9.22)
f3
=
f(ti + ∆t, xi + ∆t f2 )
(9.23)
By comparing with Taylor's expansion, it can be shown that RK4 has a local error = O(∆t5) and global error = O(∆t4).
For f(t,x) = f(t), i.e. independent of x, RK4 reduces to Simpsion's rule.

Chapter 10
ODE: boundary value problem

10.1  Second-order ODE

    Equation of motion:
Mx" = − k x + Mg − λx′ + F(t)
(10.1)
Now, assume that M=0.1 kg, k=10 N/m, λ=0.7 Ns/m, and F(t)=0.
Also assume boundary conditions (B.C.):
x(t)
=
0 cm     
at t=0 s
(10.2)
x(t)
=
2 cm
at t=0.1 s
(10.3)
Find x(t).
In general, consider the boundary value problem of 2nd-order ODE
x"(t)
=
f( t, x(t), x′(t) )       
differential equation (DE)
(10.4)
x(t0)
=
x0
boundary conditions (B.C.)
(10.5)
x(t1)
=
x1
(10.6)
The aim is to find x(t) for t0 ≤ t ≤ t1 numerically.

10.2  Shooting method

    For any estimate v0 of x′(t0), we can solve the initial value problem (using e.g. RK4):
x"(t)
=
f( t, x(t), x′(t) )
(10.7)
x(t0)
=
x0
(10.8)
x′(t0)
=
v0
(10.9)
Thus, we can define g(v0) = x(t1), because x(t1) is a function of v0.
To solve the boundary value problem, solve g(v0) − x1 = 0 using e.g. bisection method.
Each function evaluation of f in the bisection method involves performing a complete RK4 integration.
Example: (Advanced!!) Forced and damped harmonic oscillator (continued: boundary value problem)
Now, assume that M=0.1 kg, k=10 N/m, λ=0.7 Ns/m, and F(t)=0. At time t=0, x=0 cm. At t=0.1 s, x=2 cm. Find x(t).
Solution: We basically combine the Euler method solution (for 2nd order ODE) with the bisection method solution.
shoot-oscillator.py
# Shooting method for damped oscillator, with Euler method
# based on bisect.cpp & Euler-oscillator.cpp
import numpy as np

M = 0.1      # mass
k = 10       # spring constant
g = 9.8      # gravity
Lambda = 0.7 # damping coefficient
T = 0.1      # final time
dt= 0.001    # time step of Euler method
x0 = 0.            # intial position 
x1 = 2             # final position 
eps=1e-12    # tolerance for machine error

def f(t, x, v):
    return -k/M*x +g - Lambda/M*v   # particle acceleration

def xfinal(v0):   # return x(T) as calculated by Euler method
    x=x0
    v=v0
    for t in np.arange(0, T+eps, dt):
        # print(t, x)     # output t, x(t)
        dx = dt * v              # Euler method iteration     
        dv = dt * f(t, x, v)
        x += dx
        v += dv
    return x

def F(v0):  # solving F(v0)=0 (i.e. x(T)=x1) with bisection 
    return xfinal(v0)-x1 

a=0; b=100  # initial range of v0
N=10        # no. of bisection iterations
for i in range(N):
    c=(a+b)/2
    if F(a)*F(c)<0: 
      b=c # take left half of the interval
    else:
      a=c # take right half of the interval
    print (a, b)
( shoot-oscillator.py )     
shoot-oscillator.cpp
// Shooting method for damped oscillator, with Euler method
// based on bisect.cpp & Euler-oscillator.cpp
#include <iostream>
#include <iomanip>
#include <math.h>
using namespace std;

const double M = 0.1;      // mass
const double k = 10;       // spring constant
const double g = 9.8;      // gravity
const double lambda = 0.7; // damping coefficient
const double T = 0.1;      // final time
const double dt= 0.001;    // time step of Euler method
double x0 = 0.;            // intial position 
double x1 = 2;             // final position 
const double eps=1e-12;    // tolerance for machine error

double f(double t, double x, double v){ 
  return -k/M*x +g - lambda/M*v;   // particle acceleration
}

double xfinal(double v0){ // return x(T) as calculated by Euler method
  double x=x0;
  double v=v0;
  for (double t=0; t<=T+eps; t+=dt){
    // cout << t << " " << x << endl;   // output t, x(t)
    double dx = dt * v;              // Euler method iteration     
    double dv = dt * f(t, x, v);
    x += dx;
    v += dv;
  }
  return x;
}

double F(double v0){  // solving F(v0)=0 (i.e. x(T)=x1) with bisection 
   return xfinal(v0)-x1; 
}

int main(){
  double a=0, b=100;  // initial range of v0
  const int N=10;     // no. of bisection iterations
  cout << setprecision(12);
  for (int i=0; i<N; ++i){
    double c=(a+b)/2.;  
    if ( F(a)*F(c)<0 ) 
      b=c; // take left half of the interval
    else 
      a=c; // take right half of the interval
    cout << a << "  " << b << endl;
  }
}
( shoot-oscillator.cpp )     

Chapter 11
Steady-state heat flow and electric potential

11.1  Heat Flow in 1D

A metal rod of length L is heated at one end. Wait until temperature stabilizes, i.e. steady state.
Let u(x) be the temperature at x ( 0 ≤ x ≤ L ), i.e. non-uniform temperature.
Heat current j(x) at position x follows, ( after letting k=thermal conductivity)
j(x) = − k du

dx
(11.1)
Conservation of energy implies
d j

dx
=0
(11.2)
Combining them gives the 1D steady-state heat flow equation, which is the 1D Laplace's Equation:
d2 u(x)

dx2
= 0
(11.3)
Boundary conditions (e.g.):
u(0)
=
700° C
(11.4)
u(L)
=
200° C
(11.5)
This is a 2nd-order ODE, boundary value problem. Solve by shooting method with the central difference approximation:
⇒ u"(x)
=
u(x+∆x) + u(x−∆x) − 2 u(x)

∆x2
+ O(∆x2)
(11.6)

11.2  Heated plate and 2D Laplace's Equation

    A metal plate of size Lx ×Ly is heated from one edge.
Let u(x,y) be the temperature at (x,y)     ( 0 ≤ x ≤ Lx    0 ≤ y ≤ Ly ).
Then, u(x,y) follows the 2D steady-state heat flow equation, which is the 2D Laplace's equation:
2 u(x,y)

∂x2
+ 2 u(x,y)

∂y2
= 0
(11.7)
The Laplace's equation in vector form is, with u(r ) ≡ u(x,y)
2 u = 0
(11.8)
It is a 2nd-order partial differential equation (PDE).
Applying central difference approximation to both derivatives, the PDE becomes
u(x+∆x,y) + u(x−∆x,y ) − 2 u(x,y)

∆x2
+ u(x, y+∆y) + u(x, y−∆y) − 2 u(x,y)

∆y2
= 0
(11.9)
Putting ∆x = ∆y for convenience, we have
u(x,y) = 1

4
[ u(x+∆x,y) + u(x−∆x,y ) + u(x, y+∆y) + u(x,y−∆y) ]
(11.10)
    Applying a spatial discretization,
xi
=
i ∆x
(11.11)
yj
=
j ∆y
(11.12)
ui,j
=
u(xi, yj)
(11.13)
where ∆x = Lx/Nx and ∆y = Ly/Ny.                                      
Finite difference equation for 2D heat flow becomes:
ui,j = 1

4
[ ui+1,j + ui−1,j + ui,j+1 + ui,j−1 ]
(11.14)
Boundary conditions:
    e.g.:
u(x,y)
= 700° C         
for x=0
(11.15)
(left edge)
(11.16)
u(x,y)
= 200° C
for x=Lx or y=0 or y=Ly
(11.17)
(right, bottom or top edges)  
(11.18)
Discretized form:
ui,j
= 700° C         
for i=0
(11.19)
ui,j
= 200° C
for i=Nx or j=0 or j=Ny
(11.20)
                                                                                                               

11.3  Gauss-Seidel relaxation method

Iterative method to obtain ui,j for all internal points (i=1,2, .., Nx−1;    y=1,2, .., Ny−1)
Initial guess (0th iteration), e.g.:
u0i,j = 200° C             for i=1,2, .., Nx−1;    y=1,2, .., Ny−1
(11.21)
k-th iteration with Jacobi relaxation method (slow):
uk+1i,j = 1

4
[ uki+1,j + uki−1,j + uki,j+1 + uki,j−1 ]
(11.22)
k-th iteration with Gauss-Seidel relaxation method (fast, and easy to program):
uk+1i,j = 1

4
[ uai+1,j + ubi−1,j + uci,j+1 + udi,j−1 ]
(11.23)
where a, b, c, d = k or k+1 specifies the most updated value of u,   e.g.
uk+1i,j = 1

4
[ uk i+1,j + uk+1i−1,j + uki,j+1 + uk+1i,j−1 ]
(11.24)
At convergence (i.e. uk+1i,j = uki,j ), solution is ui,j = uki,j
Programming tips:
Use only one array u[Nx+1][Ny+1] to store only the most updated value of uki,j
C++ code for one Gauss-Seidel step:
u[i][j]=1./4.*(u[i+1][j]+u[i-1][j]+u[i][j+1]+u[i][j-1]);
Algorithm
Accuracy
To obtain good accuracy:

11.4  Electric Potential

Assume a charge density ρ(x,y) and an electric potential V(x,y) at (x,y) in 2D. Then V(x,y) follows the 2D Poisson's equation
2 V(x,y)

∂x2
+ 2 V(x,y)

∂y2
= − ρ(x,y)

ϵ
(11.25)
where ϵ = electric permittivity of the medium.
It can similarly be discretized into finite difference equations and solved by Gauss-Seidel relaxation method.
It reduces to Laplace equation when ρ(x,y) ≡ 0.
See this java applet for a demonstration about a parallel plate capacitor.

Chapter 12
Time-dependent Heat flow and waves on a string

12.1  Heated rod (parabolic PDE)

    A metal rod of length L is heated at one end.
Let u(x,t) be the temperature at x ( 0 ≤ x ≤ L ) at time t.
Then, u(x,t) follows the 1D time-dependent heat flow equation, also called the diffusion equation,
∂u(x,t)

∂t
= D 2 u(x,t)

∂x2
      (parabolic PDE)
(12.1)
where D=k/C is the heat diffusion coefficient which depends on the thermal conductivity k and the heat capacity C.
Initial conditions (e.g.):
u(x,0)
= 200° C       for 0 < x < L
(12.2)
Boundary conditions (e.g.):
u(0,t)
= 700° C       for t > 0
(12.3)
u(L,t)
= 200° C       for t > 0
(12.4)

12.2  Forward Euler's method

Discretizing space: xi = i∆x
Discretizing time: tk = k∆t
Then, uik = u(xi, tk).
Finite difference equation for time-dependent heat flow:
u(x, t+∆t) − u(x,t )

∆t
= D u(x+∆x,t) + u(x−∆x,t ) − 2 u(x,t)

∆x2
(12.5)
⇒   uik+1 − uik

∆t
= D ui+1k + ui−1k − 2 uik

∆x2
(12.6)
This gives the Forward Euler's equation
uik+1 = D ∆t

∆x2
( ui+1k + ui−1k ) + (1 − 2 D ∆t

∆x2
) uik       for k=0, 1, 2, 3 ....
(12.7)
Algorithm
Accuracy To obtain good accuracy:

12.3  Waves along an elastic string (hyperbolic PDE)

    Let y(x,t) be the transverse displacement at position x and time t.
It follows the wave equation (which is a hyperbolic PDE)
2 y(x,t)

∂t2
= v2 2 y(x,t)

∂x2
(12.8)
where v = (T/μ)1/2 is the wave velocity, T is the tension, and μ is the mass per unit length of the string.
Initial conditions (e.g.):
 
y(x,0) = A sin(k x)        
for 0 < x < L, k = π/L
(12.9)
 
∂y

∂t
(x,0) = 0
for 0 < x < L
(12.10)
Boundary conditions (e.g.):
y(0,t) = y(L,t) = 0       for t > 0
(12.11)

12.4  Forward Euler's method

Discretizing space: xi = i∆x
Discretizing time: tk = k∆t
Then, yik = y(xi, tk).
Finite difference equation for the wave equation:
y(x, t+∆t) + y(x, t−∆t) − 2 y(x,t )

∆t2
=
v2 y(x+∆x,t) + y(x−∆x,t ) − 2 y(x,t)

∆x2
(12.12)
⇒   yik+1 + yik−1 − 2 yik

∆t2
=
v2 yi+1k + yi−1k − 2 yik

∆x2
(12.13)
This gives the forward Euler's equation
yik+1 = v2 ∆t2

∆x2
( yi+1k + yi−1k ) + 2 (1 − v2 ∆t2

∆x2
) yik − yik−1       for k=1, 2, 3 ....
(12.14)
For k=0, initial conditions gives
yi1 − yi0

∆t
=
0
(12.15)
⇒ yi1
=
yi0
(12.16)
Algorithm
wave.py
# wave.py: Wave motion on an elastic string
from math import *
import numpy as np

v=25.    # velocity of wave in m/s
L=1.     # length in m
A=0.02   # initial ampiltude in m

N=20     # number of spatial segments
dx=L/N   # spatial discretization
kmax=80  # number of time steps to simulate
dt=0.001 # temporal discretization

yk_1 = np.array([0.]*(N+1))
yk = np.array([0.]*(N+1))
yk1 = np.array([0.]*(N+1))
#yk_1[N+1], yk[N+1], yk1[N+1] # yi at time step k-1, k and k+1

for i in range(N+1):         # all points
    x = i*dx
    yk_1[i] = A*sin(pi*x/L)  # initial displacment and boundary conditions
    yk[i] = yk_1[i]          # becasue initial velocity = 0
a = (v*dt/dx)**2 

for k in range(1,kmax+1):
    for i in range(N+1):
        print (i*dx, yk[i])
    print ()

    # forward Euler
    for i in range(1,N):    # internal points
      yk1[i] = a * (yk[i+1] + yk[i-1]) + 2*(1-a)*yk[i] - yk_1[i] 

    yk_1 = np.copy(yk)      # copy y(k) to y(k-1)
    yk   = np.copy(yk1)     # copy y(k+1) to y(k)
( wave.py )     
wave.cpp
// wave.cpp: Wave motion on an elastic string
#include <iostream>
#include <math.h>
using namespace std;

const double v=25.;    // velocity of wave in m/s
const double L=1.;    // length in m
const double A=0.02;  // initial ampiltude in m

const int N=20; // number of spatial segments
const double dx=L/N;   // spatial discretization
const int kmax=80; // number of time steps to simulate
const double dt=0.001; // temporal discretization

double yk_1[N+1], yk[N+1], yk1[N+1]; // yi at time step k-1, k and k+1

int main(){

  for (int i=0; i<=N; ++i){   // all points
    double x = i*dx;
    yk_1[i]=A*sin(M_PI*x/L);  // initial displacment and boundary conditions
    yk[i]=yk_1[i];            // becasue initial velocity = 0
  }
  double a = pow( v*dt/dx, 2. ); 

  for (int k=1; k<=kmax; ++k){
    for (int i=0; i<=N; ++i) 
      cout << i*dx << " " << yk[i] << endl;     // output result
    cout << endl;

    // forward Euler
    for (int i=1; i<N; ++i)   // internal points
      yk1[i] = a * (yk[i+1] + yk[i-1]) + 2*(1-a)*yk[i] - yk_1[i]; 

    for (int i=1; i<N; ++i){
      yk_1[i] = yk[i];     // copy y(k) to y(k-1)
      yk[i] = yk1[i];      // copy y(k+1) to y(k), 
    }
  }
}

( wave.cpp )     



File translated from TEX by TTH, version 4.15.
On 8 Oct 2021, 13:55.