Numerical Methods

EE1103
Instructor: Anil Prabhakar
Office: ESB 246E
anilpr@[Link]
Installing Linux

Install Virtualbox or WSL

● Ubuntu LTS

Check some basics

● cat
Survey1

Linux commands
C language: Pointers, Arrays, Structures

● Good tutorial on Pointers:

[Link]
● Structures: C-Lecture slides
● Writing pseudo code
● Compiling your code: gcc file.c -o executable -lm
○ #include<math.h>
● Profiling your code
○ Use gprof
○ Use breakpoints
Outline
● Linux commands
● Pseudo code of Random Walk - Assignment 1
● Understanding Noise
○ Uniform Distribution
○ Normal distribution - Box Muller Transformation
○ Generating random numbers - Assignment 2
● Truncation Errors
○ Accuracy versus precision
○ Taylor series approximation - Assignment 3
● Bracketing methods to find the root of f(x) = 0
● Iterative Methods - Assignment 6 (Mandelbrot Set)
Outline

● Interpolation - Newton, Lagrange, Quadratic and Cubic Splines

● Interpolation : Splines
○ Needs matrix inversion
○ Gaussian elimination
● LU decomposition

● Processing data files

○ Simulating an experiment - Assignment 7 (Lorentzian pulse train)
○ Histogram analysis - Assignment 8 (BER and Visibility of photons)
Outline

● Solving ODEs
○ ODE , Heun, Midpoint
○ Predictor-Corrector
○ Runge-Kutta
● Regression
Tips - Pareto Principle

● 80-20 principle
○ 80% of bugs from 20% of errors in code … invest time wisely
○ 20% of learning concepts are applicable to 80% of programming tasks

Source: The Art of Clean Code, Christian Mayer

Bash shell
● ls, (cat, less, more, tail), wc, grep, awk, head
○ head -n -7 [Link] > temp && mv temp [Link]
● pwd, environment variables ($PATH)
○ echo $PATH
○ set PATH=$PATH:/home/anil/ee1103
● Compiler gcc (has flags e.g. -o)
○ gcc myprogram.c -o myprogram
● Shell scripts
#!/bin/bash
# NOTE : Quote it else use array to avoid problems #
FILES="/path/to/*"
for f in $FILES
do
echo "Processing $f file..."
# take action on each file. $f store current file name
cat "$f"
done

● Use ping to generate a data file, redirect output to a file

● Use awk to process the 7th column and calculate the mean and stdev of the times
Accuracy versus Precision

● Additive Gaussian Noise

● The power of averaging
● Differentiate between two
distributions

● How many measurements are good enough in an expt?

 Normal Distribution

Exercise

Given an array of values, use the

empirical rules to determine if the
values fall within a normal
distribution.

Create your own array of values

and check if they follow the
empirical rules of a normal
distribution.

● What is a 6-sigma manufacturing process?

● [Link]
pective-ab89fcfd29b7
Box Muller transform

● Use samples U1 and U2 from uniform distributions

○ rand()/MAXRAND
● Calculate Z0 and Z1 using U1 and U2

● Z0 and Z1 will lie on a normal distribution

● Exercise: Generate 10k samples of Z0 and Z1 and check the histogram for its
adherence to a normal distribution
○ Confirm the number of samples within 2, 4 and 6 𝛔
Assign2: Generating numbers from a distn

● Uniform distribution of
random integers (between
0 and RAND_MAX) using
rand()
● Normal distribution using
the Box-Muller transform,
starting with two
uniformly distributed
random integers
● Counting the number of
values that lie within ∓σ, ∓
2σ and ∓ 3σ of the mean
Central Limit Theorem

● Law of Large Numbers: As you increase sample size (or the number
of samples), then the sample mean will approach the population
mean.
● With multiple large number of samples, the sampling distribution of
the mean is normally distributed, even if your original variable is not
normally distributed.
● Questions
a. How do we define this large number?
b. How many bins can we use to create a histogram of N values?

[Link]
Generating random numbers

● Functions like rand() and srand() are available in the math library to
generate streams of random numbers. These can be converted into
binary streams.
○ Extra: you can also use a linear shift register
○ Extra: can you use bits instead of integers? Bitwise operations
● Hamming distance is a measure of how different two streams of bits are.
○ Exercise: Find a codeword (of length M) within a block (of length N)
■ Write this code using pointers to elements in the arrays.
■ Use malloc() to allocate memory and free() to release it
Assign3: Truncation Errors

● Taylor Series expansions

○ Truncate it after some terms
○ Sum of residual terms in the series is the error

● Calculate the sin(x) up to N terms

● Use a function to calculate the factorial
● Extra: Attempt recursion
Extra: Generating a random bit stream

● Functions like rand() and srand() are available in the math library to
generate streams of random numbers. These can be converted into bit
streams.
○ You can also use a linear shift register
● Hamming distance is a measure of how different two streams of bits are.
○ Try to write this code using pointers to elements in the arrays.
● Good to write programs that take in multiple inputs with flags.
○ The getopt function is useful. Can also use strcmp()
● The NIST test suite checks how random a bit stream is. Download, make
and run it. Read its documentation.
Bracketing Methods

● Find the root to the equation f(x) = 0 in the interval {xl, xu}
○ Bisection Method
○ False Position Method
● Determine the root xr within the error
Newton - Raphson Method

● Use the slope to estimate the next point

● Can sometimes have a problem
Secant Method

● Use the slope from the previous point

Summary
 Hailstone Numbers
• A Hailstone Sequence is generated by a simple algorithm:

Start with an integer N. If N is even, the next number in the sequence is N/2. If
N is odd, the next number in the sequence is (3*N)+1

• 7, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1, 4, 2, 1, ...
repeats
• 12, 6, 3, 10, 5, 16, 8, 4, 2, 1, 4, 2, 1 ….
• 909, 2726, 1364, 682, 341, 1024, 512, 256, 128, 64, 32, 16, 8, 4, 2,
1, 4, 2, 1…
 Hailstone Numbers
[Link]

Exercise : Write a program to accept an input and count the number of iterations needed to get to
1, and the highest number reached. Generate a table of results…
Fixed point iteration

● Logistic map
● Montonone, oscillating or spiral patterns
● x[n+1] = r x[n](1-x[n])
● Convergence requires |g’(x)| < 1
● [Link]
Gnuplot

● Simple plotting program

○ p “myfi[Link]” using 1:3 with points … plots columns 3 vs 1
○ help plot … gives you help on different options in plot
○ Save “[Link]” … saves the gnuplot commands used in this session
○ An empty line in myfi[Link] will treat subsequent points as a new plot
○ set xrange [0:20] … sets the range of x-axis

● Incorporate gnuplot into your c-code

FILE *pipe_gp = popen("gnuplot -p", "w");
fputs("set terminal png \n",pipe_gp);
fputs("set output '[Link]' \n",pipe_gp);
fputs("set xlabel 'f' \n",pipe_gp);
fputs("set xrange [0:100] \n",pipe_gp);
fputs("set yrange [0:100] \n",pipe_gp);
fputs("plot '[Link]' u 1:2 w circles lc rgb 'pink' notitle \n",pipe_gp);
pclose(pipe_gp);
AWGN

● Additive White Gaussian Noise

○ Additive - is added to the signal, x(t) = s(t) + n(t)
○ White - uniform spectral content
○ Gaussian - amplitude distribution is a Gaussian

● Exercise: Generate a Lorentzian pulse train (time series) and add

noise to it
○ Amplitude noise
○ Uncertainty in location of the Lorentzian pulses
○ Uncertainty in the width of the Lorentzian pulses
Filtering Noisy Data

● Moving average filter

○ How does the window size affect the output?
● Exponential filter
○ Single exponential
■ y(k) = a * y(k-1) + (1-a) * x(k)
■ y(k) = y(k-1) + (1-a)*( x(k) - y(k-1) ) … Predictor-Corrector form
○ Double exponential
■ y[k] = a ( y[k-1] + b[k-1] ) + (1-a) x[k]
■ b[k] = g b[k-1] + (1-g) ( y[k] - y[k-1] )

● Alpha - beta filter

Time series data analysis

● Exercise: Generate and analyze a noisy pulse train (time series)

○ Inputs: number of pulses, time between pulses, width of pulses, type
of pulses (Gaussian or Lorentzian)
○ Add baseline amplitude noise, location noise, width noise
○ Extract the average time between pulses, and the average width of
the pulses. You can use
■ a window filter to reduce the noise
■ a threshold to identify the FWHM and location of the peak
○ Compare the extracted values of <T> and <a> from the inputs to
your code
○ How does the difference between extracted and input values
change with the different types of noise? What are <T> and <a>
most sensitive to?
Interpolation: Newton’s Polynomial

Bracketed terms are

finite difference
Interpolation: Lagrange Polynomial

Second order
polynomial

● Careful with division by small numbers

● Will sometimes give numerical errors for higher order polynomials
Newton and Lagrange interpolation

● Convert to C-code
● Can you use recursion?
● Check the pseudo code for bugs
before using it
Interpolation: Splines

● Interpolating polynomial cause “ringing”

● Splines are better
○ Lower order polynomials to subset of points
Quadratic Splines

● Adjacent polynomials are

equal at interior knots
● First and last functions
pass through end points
● First derivative at interior
knots must be equal

● n intervals, need 3n
conditions
Cubic Splines

● Fit a cubic polynomial between points

● Identify the number of intervals
● Identify the number of unknowns
● Find the corresponding number of equations by equating the
interior points and the derivatives
● Leave the second derivative at the end points unknown
Gauss Elimination

● Numerically solve a system of equations

Gauss Elimination - Pseudo code

Forward elimination Backward substitution

● Assignment: Convert this to C code (check for bugs)

● Reading: LU decomposition
Gauss Elimination - Pitfalls

● Division by zero

● Round off errors

○ Large numbers of equations

● Ill conditioned systems

○ Small changes in coefficients result in large changes in values
Gauss Elimination - Improvements

● Normalization
● Partial Pivoting
○ Find the largest coefficient in
the row and interchange
columns
Gauss Elimination

● Understand
○ check for bugs
● Convert the pseudo
code to C-code
LU decomposition

Multiply row 1
Subtract result from row 2

Multiply row 1
Subtract result from row 3

Multiply row 2
Subtract result from row 3

Sec. 10.1 Chapra

Regression

● Simple examples
● Which is a better fit?

[Link]
Linear Fits

● Goodness of a Fit : R^2

● Can R^2 be biased?

[Link]
Linear Fits - Sum of Residuals

Doesn’t work well

● Minimize the sum of squares of the residual

Linear Fit - Least Squares Minimization

● Set the partial derivatives of S to zero and solve for a0 and a1

Standard Error of the Estimate

Spread about
regression line

St is the sum of squares

around the mean of y

r is the correlation coeffi`cient

● Always check that the errors are normally distributed about your fit
● Else, remove an outlier
Curve fits: Lorentzian versus Gaussian

● Start with N points of additive white Gaussian noise n(x) with given standard
deviation σ (recall the Box-Muller algorithm that creates Gaussian noise and
normalize it to ±3σ)
● Create a Lorentzian L(x) that approximately fills N points, add it to the noise
f(x) = L(x)+ n(x)
● Fit f(x) to a Gaussian g(x) and extract R2
● Plot log(R2) versus σ and extract the power-law dependence

Learning Outcomes:
1. Reuse code on the Box-Muller transform
2. Work with gnuplot or the Gnu Scientific Library
3. Understand the limitations of curve fitting on noisy data
Floating point operations or FLOPS

● Usually focus on multiplication and division

Using hexadecimals

● Hexadecimal notation is easily converted to-from binary

● gcc understands that 0xff = 255
Bitwise Operations

& bitwise AND x=x & ~077(zeros last 6 bits of x)

! bitwise (inclusive) OR x=x ! mask (sets ON the 1 bits of the mask)
^ bitwise (exclusive) OR x^y (sets to 1 all positions where
x and y differ and 0 where they agree)
<< left shift x << 2 shifts bits in x to left, padding with zeros
>> right shift x >> 2 shifts bits in x to right, padding with 0s.
~ one’s complement of the argument ~x converts all 1 bits to 0 and 0 bits to 1

These operators have higher precedence than &&, ||, ?:, = etc but lower precedence than ==,
!=, <, <=, etc

These bits operations are useful to take apart the bits of a word. For instance if we know that
the 4th bit of a variable x is 1 if a certain action is to be taken (very common in operating
system situations), we “test” the variable by
if( x & 8 ){ … }
where 8 is 1000 in binary. x&8 returns 1 if the 4th bit of x is 1, and zero otherwise.
Profiling your code

● Compile your code using gcc with the -pg flag

● Running the code generates [Link]
● gprof will analyze [Link] and generate a profile

Try it out: [Link]

Note: The sample code in the gprof-tutorial uses a hexadecimal notation for 32 bits

0xffffffff = 11111111111111111111111111111111 (binary)

= 4294967295 (decimal)
Sample Size
Z = 1.95 for 95% confidence interval

[Link]
ODEs - first order method

● Euler’s formula:
● Errors
○ Round-off
■ Local
■ Global - local errors add up

○ Truncation
 ODEs : Modular Approach

Can easily change from

Euler to something else
ODEs

● Use the higher order terms in the Taylor series

● Chain differentiation if derivatives are functions of both x and y

ODEs: Heun’s method

● Predictor
● Corrector

● Use the value at i+1 to correct for slope at i

ODEs: Midpoint Rule

● Take a half step

● Calculate the slope
● Extrapolate
ODEs: Heun with iteration
ODEs: Runge Kutta

● Try to reach the accuracy of a Taylor series with an increment function ɸ

● Must calculate the unknown ai , pi and qi coefficients

● k1 is the increment based on the slope at the

beginning of the interval, using yi (Euler's method),

● k2 is the increment based on the slope at the

midpoint of the interval, using yi + h/2*k1 ,

● k3 is again the increment based on the slope at the

midpoint, but now using yi + h/2*k2
● We have infinite choice in a2
○ Heun: a2 = 0.5 ● k4 is the increment based on the slope at the end of
○ Midpoint: a2 = 1 the interval, using yi + h*k3
○ Ralston: a2 = 1/3
ODEs: 4th order Runge Kutta

● Most commonly used with 4th order corrector

 ODEs: Comparisons

● Start with a magnetization pointing almost at the N pole and apply a

field along the S-pole
● Estimate the trajectory of the magnetization on the unit sphere using
○ Euler
○ Heun

● Assume that RK45 gives you the most accurate results

○ Calculate the R2 error with the Euler and Heun trajectories
○ Use RK45 to estimate the switching time 𝛕
■ plot 𝛕 as a function of 𝜶
In the limit of H >> Hk, the switching time can be estimated analytically.
E.g. see (11) of Mallinson, “Damped Gyromagnetic Switching”, doi:10.1109/20.875251
Second order differential equations

● Swinging pendulum
● [Link]
Phase Portrait

● Nonlinear pendulum - chaotic

● [Link]