Full HTML for

Scripted foilset CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm

Given by Geoffrey C. Fox at Delivered Lectures of CPS615 Basic Simulation Track for Computational Science on 10 October 96. Foils prepared 28 December 1996
Outside Index Summary of Material Secs 44.6

This discusses solution of systems of ordinary differential equations (ODE) in the context of N squared particle dynamics problems
We start with motivation with brief discussion of relevant problems
We go through basic theory of ODE's including Euler, Runge-Kutta, Predictor-Corrector and Multistep Methods
We begin the discussion of solving N body problems using classic Connection Machine elegant but inefficient algorithm
Note -- Some foils expanded to two after talk and second one is given without audio in script

Table of Contents for full HTML of CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm

 Denote Foils where Image Critical
Denote Foils where HTML is sufficient
Indicates Available audio which is lightpurpleed out if missing
1 Delivered Lectures for CPS615 -- Base Course for the Simulation Track of Computational Science
Fall Semester 1996 --
Lecture of October 10 - 1996
2 Abstract of Oct 10 1996 CPS615 Lecture
3 Particle Applications - Ordinary Differential Equations (ODE's)
4 Particle Applications - the N-body problem
5 Newton's First Law -- The Gravitational Force on a Particle
6 Equations of Motion -- Newton's Second Law
7 Numerical techniques for solving ODE's
8 Second and Higher Order Equations
9 Basic Discretization of Single First Order Equation
10 Errors in numerical approximations
11 Runge-Kutta Methods: Euler's method
12 Estimate of Error in Euler's method
13 Relationship of Error to Computation
14 Example using Euler's method from the CSEP book
15 Approximate solutions at t=1,using Euler's method with different values of h
16 Runge-Kutta Methods: Modified Euler's method
17 Approximate solutions of the ODE for et at t=1, using modified Euler's method with different values of h
18 The Classical Runge-Kutta -- In Words
19 The Classical Runge-Kutta -- Formally
20 The Classical Runge-Kutta Pictorially
21 Predictor / Corrector Methods
22 Definition of Multi-step methods
23 Features of Multi-Step Methods
24 Comparison of Explicit and Implicit Methods
25 Solving the N-body equations of motion
26 Representing the N-Body problem
27 Form of the Computation
28 N-body Runge Kutta Routine in Fortran90 - I
29 Runge Kutta Routine in Fortran90 - II
30 Computation of accelerations - a simple parallel array algorithm
31 Simple Data Parallel Version of N Body Force Computation -- Grav -- I
32 The Grav Function in Data Parallel Algorithm - II

Outside Index Summary of Material

HTML version of Scripted Foils prepared 28 December 1996

Foil 1 Delivered Lectures for CPS615 -- Base Course for the Simulation Track of Computational Science
Fall Semester 1996 --
Lecture of October 10 - 1996

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index
Geoffrey Fox
NPAC
Room 3-131 CST
111 College Place
Syracuse NY 13244-4100
HTML version of Scripted Foils prepared 28 December 1996

Foil 2 Abstract of Oct 10 1996 CPS615 Lecture

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 44.6
This discusses solution of systems of ordinary differential equations (ODE) in the context of N squared particle dynamics problems
We start with motivation with brief discussion of relevant problems
We go through basic theory of ODE's including Euler, Runge-Kutta, Predictor-Corrector and Multistep Methods
We begin the discussion of solving N body problems using classic Connection Machine elegant but inefficient algorithm
Note -- Some foils expanded to two after talk and second one is given without audio in script
HTML version of Scripted Foils prepared 28 December 1996

Foil 3 Particle Applications - Ordinary Differential Equations (ODE's)

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 211.6
Consider Models of physical systems represented as sets of particles rather than densities (fields) evolving over time
Examples:
  • the solar system
  • an electrical or electonic network
  • a globular star cluster
  • the atoms in a crystal vibrating under interatomic forces
  • the molecule rotating and flexing under interatomic force
Laws of Motion are typically ordinary differential Equations
Ordinary means differentiate wrt One Variable -- typically time
HTML version of Scripted Foils prepared 28 December 1996

Foil 4 Particle Applications - the N-body problem

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 40.3
N particles, each with a mass mi , moving with velocity Vi through 3-dimensional space
Are governed by Newton's equations of motion
Basic Kinematics
HTML version of Scripted Foils prepared 28 December 1996

Foil 5 Newton's First Law -- The Gravitational Force on a Particle

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 36
HTML version of Scripted Foils prepared 28 December 1996

Foil 6 Equations of Motion -- Newton's Second Law

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 72
Newton's Second Law
Incorporate laws into equations of motion
Example of force law for molecular dynamics
HTML version of Scripted Foils prepared 28 December 1996

Foil 7 Numerical techniques for solving ODE's

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 195.8
ODE's give an equation for the derivative of X with respect to time t
They can be classified by the boundary conditions used
Initial value problems
  • Give initial value Xo for X and first derivative X'o
  • of X, X', at A
Boundary value
  • Give endpoint
  • values X o X1
  • for X
  • at A and B
HTML version of Scripted Foils prepared 28 December 1996

Foil 8 Second and Higher Order Equations

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 92.1
Second order equations such as
can always be rewritten as a system of 2 first-order equations involving X and a new variable Y representing the first order derivative:
HTML version of Scripted Foils prepared 28 December 1996

Foil 9 Basic Discretization of Single First Order Equation

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 135.3
For simplicity, we assume just one first-order equation, where f is the function on right hand side which depends on X and t
We can always solve this by setting up a grid of equidistant points with grid size h=(B-A)/n where n is an integer.
Starting from the initial value, we calculate positions one step at a time
HTML version of Scripted Foils prepared 28 December 1996

Foil 10 Errors in numerical approximations

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 300.9
Two sources of error:
Computational error includes such things as roundoff error, etc. and is generally controlled by having enough significant digits in the computer arithmetic
Discretization error is the accuracy of the numerical method and has two measures:
  • Global error - difference between solution and numerical approximation at any point. Difficult or expensive to estimate.
  • Local error - difference between solution and numerical approximation over only one time step. Easier to calculate and indicates global error indirectly.
HTML version of Scripted Foils prepared 28 December 1996

Foil 11 Runge-Kutta Methods: Euler's method

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 148.3
Euler's method is not practical, but illustrates the technique.
It involves Linear approximation to get next point
HTML version of Scripted Foils prepared 28 December 1996

Foil 12 Estimate of Error in Euler's method

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 260.6
Use Taylor's theorem to represent the exact solution:
HTML version of Scripted Foils prepared 28 December 1996

Foil 13 Relationship of Error to Computation

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 60.4
Whenever f satisfies certain smoothness conditions, there is always a sufficiently small step size h such that the difference between the real function value at ti and the approximation Xi+1 is less than some required error magnitude e. [Burden and Faires]
Euler's method: one computation of the derivative function f at each step.
Other methods require less computation in order to produce the specified error e.
HTML version of Scripted Foils prepared 28 December 1996

Foil 14 Example using Euler's method from the CSEP book

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 84.9
Initial Value Problem with known analytical solution:
The approximation with Euler's method and h=0.25:
Calculate a few values of the approximation:
HTML version of Scripted Foils prepared 28 December 1996

Foil 15 Approximate solutions at t=1,using Euler's method with different values of h

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 155.5
Note it will take about one million iterations to get an error of order O(10-6)
The last column shows global error is of order O(h) as expected
HTML version of Scripted Foils prepared 28 December 1996

Foil 16 Runge-Kutta Methods: Modified Euler's method

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 269.2
Use the derivative at one time step to extrapolate the midpoint value - use midpoint derivative to extrapolate the function value at the next time step
Evaluates the derivative function twice at each time step. Global error - O(h2), second order method
Sometimes called the midpoint method
HTML version of Scripted Foils prepared 28 December 1996

Foil 17 Approximate solutions of the ODE for et at t=1, using modified Euler's method with different values of h

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 252
Note global error is now O(h2) and we get an error of O(10-5) after 128 iterations which would take about 1000 more iterations for Euler's method to achieve
So Euler has roughly half the computational effort per iteration but requires the square of the number of iterations
HTML version of Scripted Foils prepared 28 December 1996

Foil 18 The Classical Runge-Kutta -- In Words

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 253.4
Runge Kutta methods achieve better results than Euler by using intermediate computations at intermediate time values
The fourth-order rule is the favorite method as it achieves good accuracy with modest computational complexity -- the algorithm is in words:
Use derivative of first time step to get trial midpoint
Use its derivative at first time step to get second trial midpoint
Use its derivative to get a trial end point
Integrate by Simpon's Rule, using average of two midpoint estimates
Global error is fourth order
HTML version of Scripted Foils prepared 28 December 1996

Foil 19 The Classical Runge-Kutta -- Formally

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index
HTML version of Scripted Foils prepared 28 December 1996

Foil 20 The Classical Runge-Kutta Pictorially

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 142.5
Compared with Euler, Runge-Kutta has 4 times more calculation per time step, but should use fourth root as many time steps
HTML version of Scripted Foils prepared 28 December 1996

Foil 21 Predictor / Corrector Methods

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 169.9
First, predict Xi+1 using an explicit equation, with O(hn) error, and known values.
Then correct this value by using it in an implicit equation, with O(hn+1) error.
Simple example:
HTML version of Scripted Foils prepared 28 December 1996

Foil 22 Definition of Multi-step methods

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 303.8
The predictor/corrector methods use previous values Xi-1, Xi-2, ...... to increase accuracy p - not extra values between Xi and Xi+1 as in Runge Kutta
General form of multi-step difference equation:
Xi+1 = am-1Xi + am-2 Xi-1 + ..... + a0Xi+1-m + h(bmf(ti+1, Xi+1) + bm-1 f (ti, Xi) + .....
  • + b0f(ti+1-m, Xi+1-m)))
if coefficient of f evaluation dependent on ti+1 bm=0, this is an explicit equation
HTML version of Scripted Foils prepared 28 December 1996

Foil 23 Features of Multi-Step Methods

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index
Note these are essentially interpolation formulae for if one uses information from m t values, you can fit a degree m-1 polynomial
Taylor expansion and Polynomial fitting are essentially the same thing!
Implicit Multistep Methods are gotten using backwards difference interpolating polynomials starting at ti+1. But wherever we we need f(ti+1,y(ti+1)), we use f(ti+1, X*i+1) where X*i+1is derived from the explicit predictor equation
Note that implicit formulae should be best as explicit method involves extrapolation from ti to ti+1 whereas in implicit case ti+1 is endpoint of region in which interpolation done
Extrapolation is always unreliable and to be avoided!
HTML version of Scripted Foils prepared 28 December 1996

Foil 24 Comparison of Explicit and Implicit Methods

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 331.2
Adams-Bashforth fourth-order explicit method (4 step)
Adams-Moulton fourth order method (3 step) is an implicit Multistep method
Note coefficient of implicit method is a factor of 251/19 smaller than explicit method of same order -- this is why extrapolation is not so good!
HTML version of Scripted Foils prepared 28 December 1996

Foil 25 Solving the N-body equations of motion

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 221.7
Introduce 3 vectors X V A for position,velocity and acceleration of particles
Numerical techniques iterate equations over time using
Euler's Equation gives:
X(t+dt) = X(t) + h * V(t)
V(t+dt) = V(t) + h * Grav(X(t))
Note i labels particles not time steps
HTML version of Scripted Foils prepared 28 December 1996

Foil 26 Representing the N-Body problem

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 79.2
Positions and velocities are 3 N dimensional arrays, X and V
other variables
  • M - spread masses to another 3xN array
  • h - step size
  • ns - number of time steps
subroutine for numerical method will take these arguments and update X and V.
A subroutine called Grav is assumed to compute new accelerations
HTML version of Scripted Foils prepared 28 December 1996

Foil 27 Form of the Computation

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 51.8
Computation of numerical method is inherently iterative: at each time step, the solution depends on the immediately preceding one.
At each time step, Grav is computed:
  • For each particle i, sum over the forces due to all other particles j.
  • Computation is O(N2) - potential for parallel computation.
We will use 4th order Runge Kutta to integrate in time and the program is designed an overall routine looping over time with parallelism hidden in Grav routine
HTML version of Scripted Foils prepared 28 December 1996

Foil 28 N-body Runge Kutta Routine in Fortran90 - I

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 92.1
C Solves Newton's equations of motion using Runge-Kutta method
C which is globally 4th order. X and V are initial positions and
C velocities. The system is evolved over a time interval h*ns.
C X and V contain the updated state at that time.
  • real, dimension (1:3, 1:N) :: X,V, Xdelta1, Vdelta1, Xdelta2, Vdelta2, Xdelta3, Vdelta3, Xdelta4, Vdelta4
  • real h, h2
  • integer ns, k
  • INTERFACE
  • function Grav(X,M)
  • real array(1:3, 1:N) :: X, M , Grav
  • end function
  • end INTERFACE
C Grav is hard parallel algorithm and will be given later!
HTML version of Scripted Foils prepared 28 December 1996

Foil 29 Runge Kutta Routine in Fortran90 - II

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index
HTML version of Scripted Foils prepared 28 December 1996

Foil 30 Computation of accelerations - a simple parallel array algorithm

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 126.7
Spread positions of particles into 2 3D arrays so that extra dimension is labelled by index in sum over particles that interact with a given particle
Xj is essentially transpose of Xi in second and third dimension
HTML version of Scripted Foils prepared 28 December 1996

Foil 31 Simple Data Parallel Version of N Body Force Computation -- Grav -- I

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index Secs 279.3
function Grav(X,M)
C accepts positions of particles X and masses of particles M
C returns accelerations in Grav
C Uses completely parallel calculation, ignoring anti-symmetry of force
  • integer, parameter :: N =number of particles
  • real, parameter :: G = gravitational constant
  • real array (1:3, 1:N) :: X,M,Grav
! calculates acceleration on body i due to body j in entries ( :,i,j )
  • real array (1:3, 1:N,1:N) :: A, Xi, Xj, Ms, D, R
  • logical, array (1:3, 1:N,1:N) ::diag
  • integer i,j,k
! diag is true for diagonal of N by N slices
  • forall ( k=1:3, i=1:N, j=1:N ) diag = (i.eq.j)
HTML version of Scripted Foils prepared 28 December 1996

Foil 32 The Grav Function in Data Parallel Algorithm - II

From CPS615-Discussion of Ordinary Differential Equations and Start of Parallel N-Body Algorithm Delivered Lectures of CPS615 Basic Simulation Track for Computational Science -- 10 October 96. *
Full HTML Index
! set up arrays of particles Xi and particles Xj
  • Xi = spread (X, dim =3, ncopies = N )
  • Xj = spread (X, dim = 2, ncopies = N )
  • Ms = spread (M, dim=2, ncopies = N )
! displacements and Euclidean distance
  • D = Xj - Xi
  • R = sqrt (spread ( sum ( D*D, dim=1), dim =1, ncopies=3)))
! calculate accelerations for all pairs except on main diagonal
  • where (diag)
    • A=0.0
  • elsewhere
    • A=G*Ms*D / R**3
  • end where
  • Grav = sum (A, dim=3)
end function Grav
© Northeast Parallel Architectures Center, Syracuse University, npac@npac.syr.edu

If you have any comments about this server, send e-mail to webmaster@npac.syr.edu.

Page produced by wwwfoil on Fri Aug 15 1997