Computational Astrophysics

E. Larrañaga

Observatorio Astronómico Nacional

Universidad Nacional de Colombia

March 29, 2019

Outline

1 Sources of Error
Round-off Error
Truncation Error

2 Finite Differences
Differentiation of a Discrete Function
Differentiation of an Analytic Function
Round-off Error

Round-off error arises from the error inherent in representing a

floating point number with a finite number of bits in the computer.
Round-off Error

"""
Round-off Error
We find the value of epsilon for which 1. + epsilon = 1.

This gives the machine epsilon value, representing the error inherent to
representing a floating point number
"""

epsilon = 1. # Initial value for epsilon

while (1. + epsilon != 1.):

epsilon = epsilon /2.

print (epsilon)

#Iterates, halving epsilon until 1 + epsilon = 1

# Prints the value of the machine epsilon
Truncation Error

Truncation error is a feature of an algorithm.

Typically we expand the expressions about some small quantity.

When throwing away higher-order terms, there is a truncation in
the expression.

This introduces an error in the representation !!.

If the quantity we expand about truly is small, then the error is

small.
Convergence Test

Example:
Function:
f (x) = sin(x) (1)
Taylor series representation:
∞
X x 2n−1
f (x) = (−1)n−1 (2)
(2n − 1)!
n=1

Truncation for |x|  1:

x3
f (x) = x − + O(x 5 ) (3)
6
Convergence Test

"""
Truncation Error - Convergence Test
We implement a 5th order accurate approximation of the function Sin(x)

We define a procedure that calculates the difference between the value of

the Sin(x) function and the truncated approximation.

Calculating the value of this epsilon for x<1 and then taking half of
this value of x we show that the epsilon reduces by 2**5 = 32,
demonstrating 5th-order accuracy
"""

import math as m

# Definition of the function giving the truncation error

def epsilon(x):
return m.sin(x) - (x - (x**3)/6)

# Value to calculate the function

x = 0.1

# Results. We use the formating in the print function to show the results
print("\nFor x = %f the value of the truncation error is epsilon = %e" %( x, epsilon(x)))
print("For x = %f the value of the truncation error is epsilon = %e" %( x/2, epsilon(x/2)))
print("")
print("The ratio of these values is %f " %(epsilon(x)/epsilon(x/2)))
Differentiation of a Discrete Function

Collection of equally spaced points xi (such that ∆x = xi+1 − xi )

and the corresponding values fi . The derivative of this discrete
function can be calculated using
Differentiation of a Discrete Function

Collection of equally spaced points xi (such that ∆x = xi+1 − xi )

and the corresponding values fi . The derivative of this discrete
function can be calculated using

df  fi − fi−1
= (4)
dx i ∆x
or using 
df  fi+1 − fi
= . (5)
dx i ∆x
These finite difference derivatives are first order accurate.
Differentiation of a Discrete Function

Collection of equally spaced points xi (such that ∆x = xi+1 − xi )

and the corresponding values fi . The derivative of this discrete
function can be calculated using

df  fi − fi−1
= (4)
dx i ∆x
or using 
df  fi+1 − fi
= . (5)
dx i ∆x
These finite difference derivatives are first order accurate.
A second order accurate method to calculate this derivative uses
the centered difference
Differentiation of a Discrete Function

Collection of equally spaced points xi (such that ∆x = xi+1 − xi )

and the corresponding values fi . The derivative of this discrete
function can be calculated using

df  fi − fi−1
= (4)
dx i ∆x
or using 
df  fi+1 − fi
= . (5)
dx i ∆x
These finite difference derivatives are first order accurate.
A second order accurate method to calculate this derivative uses
the centered difference

df  fi+1 − fi−1
= . (6)
dx i 2∆x
Differentiation of an Analytic Function

Numerically find the derivative of an analytic function f (x)

Differentiation of an Analytic Function

Numerically find the derivative of an analytic function f (x)


0 df  ∆f
f (x0 ) = = lim (x0 ) (7)
dx x0 ∆x→0 ∆x
Differentiation of an Analytic Function

Numerically find the derivative of an analytic function f (x)


0 df  ∆f
f (x0 ) = = lim (x0 ) (7)
dx x0 ∆x→0 ∆x

Taylor expansion
∞
X f (n) (x0 )
f (x0 + δx) = δx n (8)
n!
n=0
= f (x0 ) + f 0 (x0 )δx + O(δx 2 ) (9)
Differentiation


df  f (x0 + δx) − f (x0 )
= + O(δx) (10)
dx x0 δx
Differentiation


df  f (x0 + δx) − f (x0 )
= + O(δx) (10)
dx x0 δx

Forward difference estimate for f 0 (x0 )

Differentiation of an Analytic Function

exact
1.0 left-sided
right-sided
0.8 centered

0.6

0.4

0.2
∆x
0.0
0.0 0.5 1.0 1.5 2.0 2.5 3.0
Differentiation of an Analytic Function

f (x0 + δx) − f (x0 )

f 0 (x0 ) = + O(δx) (11)
δx
First order forward difference estimate for f 0 (x0 )

f (x0 ) − f (x0 − δx)

f 0 (x0 ) = + O(δx) (12)
δx
First order backward difference estimate for f 0 (x0 )

f (x0 + δx) − f (x0 − δx)

f 0 (x0 ) = + O(δx 2 ) (13)
2δx
Second order central difference estimate for f 0 (x0 )
Differentiation of an Analytic Function

An optimal value for δx requires a balance of truncation error

(which needs a small δx) and the round-off error (which becomes
large when δx is close to the machine  ).
Differentiation of an Analytic Function

An optimal value for δx requires a balance of truncation error

(which needs a small δx) and the round-off error (which becomes
large when δx is close to the machine  ).
√
A rule-of-thumb for the election is δ ≈ .
Differentiation of an Analytic Function in an Unevenly
Spaced Grid

Equidistant grids are the conceptually simplest way of discretizing

a problem.
Differentiation of an Analytic Function in an Unevenly
Spaced Grid

Equidistant grids are the conceptually simplest way of discretizing

a problem.

They are, however, in many cases by far not the most efficient way,
since many problems need more resolution in some parts of their
domain than in others.
Differentiation of an Analytic Function in an Unevenly
Spaced Grid

Equidistant grids are the conceptually simplest way of discretizing

a problem.

They are, however, in many cases by far not the most efficient way,
since many problems need more resolution in some parts of their
domain than in others.

For example, when modeling the stellar collapse of an iron core to

a neutron star it is useful to resolve the steep density gradients
near the neutron star. In this process, a resolution of order 100 m
is required within ∼ 30 km of the origin, while a cell size of order
10 km is sufficient at radii above ∼1000 km.
Differentiation of an Analytic Function in an Unevenly
Spaced Grid

We want to evaluate f 0 (x) at x = xi .

Consider the steps δx1 = xi − xi−1 and δx2 = xi+1 − xi .
Differentiation of an Analytic Function in an Unevenly
Spaced Grid

We want to evaluate f 0 (x) at x = xi .

Consider the steps δx1 = xi − xi−1 and δx2 = xi+1 − xi .
Then,

δx22 00
f (xi + δx2 ) = f (xi ) + δx2 f 0 (xi ) + f (xi ) + O(δx23 )
2 (14)
δx 2
f (xi − δx1 ) = f (xi ) − δx1 f 0 (xi ) + 1 f 00 (xi ) + O(δx13 )
2
Differentiation of an Analytic Function in an Unevenly
Spaced Grid

We want to evaluate f 0 (x) at x = xi .

Consider the steps δx1 = xi − xi−1 and δx2 = xi+1 − xi .
Then,

δx22 00
f (xi + δx2 ) = f (xi ) + δx2 f 0 (xi ) + f (xi ) + O(δx23 )
2 (14)
δx 2
f (xi − δx1 ) = f (xi ) − δx1 f 0 (xi ) + 1 f 00 (xi ) + O(δx13 )
2
Eliminating f 00 (xi ) and solving for f 0 (xi ),

δx1 δx1 − δx2 δx2

f 0 (xi ) = f (xi+1 ) − f (xi ) − f (xi−1 ) .
δx2 (δx1 + δx2 ) δx2 δx1 δx1 (δx1 + δx2 )
(15)
Differentiation of an Analytic Function in an Unevenly
Spaced Grid

We want to evaluate f 0 (x) at x = xi .

Consider the steps δx1 = xi − xi−1 and δx2 = xi+1 − xi .
Then,

δx22 00
f (xi + δx2 ) = f (xi ) + δx2 f 0 (xi ) + f (xi ) + O(δx23 )
2 (14)
δx 2
f (xi − δx1 ) = f (xi ) − δx1 f 0 (xi ) + 1 f 00 (xi ) + O(δx13 )
2
Eliminating f 00 (xi ) and solving for f 0 (xi ),

δx1 δx1 − δx2 δx2

f 0 (xi ) = f (xi+1 ) − f (xi ) − f (xi−1 ) .
δx2 (δx1 + δx2 ) δx2 δx1 δx1 (δx1 + δx2 )
(15)
If δx1 = δx2 = δx this reduces to the standard difference equation (13).