Numerical Recipes

Numerical Recipes The Art of Scientific Computing (with some applications in computational physics)

Computer Architecture CPU Memory External Storage

Program Organization int main() { … } double func(double x) { … }

First Example #include <stdio.h> main() { printf(“hello, world\n”); } gcc hello.c (to get a.out) [Or other way depending on your OS]

Data Types #include <stdio.h> main() { int i,j,k; double a,b,f; char c, str[100]; j = 3; a = 1.05; c = ‘a’; str[0] = ‘p’; str[1] = ‘c’; str[2] = ‘\0’; printf(“j=%d, a=%10.6f, c=%c, str=%s\n”, j, a, c, str); }

Equal is not equal • x = 10 is not 10 = x • x = x + 1 made no sense if it is math • x = a + b OK, but a+b = x is not C. • In general, left side of ‘=’ refers to memory location, right side expression can be evaluated to numerical values

Expressions • Expressions can be formed with +, -, *, / with the usual meaning • Use parentheses ( …) if meaning is not clear, e.g., (a+b)*c • Be careful 2/3 is 0, not 0.666…. • Other large class of operators exists in C, such as ++i, --j, a+=b, &, !, &&, ||, ^, ?a:b, etc

Use a Clear Style (A) k=(2-j)*(1+3*j)/2; k=j+1; if(k == 3) k=0; switch(j) { case 0: k=1; break; case 1: k=2; break; case 2: k=0; break; default: { fprintf(stderr, “unexpected value for j”); exit(1); } } (D) k=(j+1)%3; (B) (C)

Control Structures in C - loop for (j=0; j < 10; ++j) { a[j] = j; } while (n < 1000) { n *= 2; } do { n *= 2; } while (n < 1000);

Control Structure - conditional if (b > 3) { a = 1; } if (n < 1000) { n *= 2; } else { n = 0; }

Control Structure - break for( ; ; ) { ... if(. . .) break; }

Pointers • Pointer is a variable in C that stores address of certain type • Int *p; double *ap; char *str; • You make it pointing to something by (1) address operator &, e.g. p = &j, (2) malloc() function, (3) or assignment, str = “abcd”. • Get the value the pointer is pointing to by dereferencing, *p

1D Array in C • int a[4]; defines elements a[0],a[1],a[2], and a[3] • a[j] is same as *(a+j),a has a pointer value • float b[4], *bb; bb=b-1; then valid range of index for b is from 0 to 3, but bb is 1 to 4.

1D Array Argument Passing void routine(float bb[], int n) // bb[1..n] (range is 1 to n) • We can use as float a[4]; routine(a-1, 4);

2D Array in C int m[13][4]; defines fixed size array. Example below defines dynamic 2D array. float **a; a = (float **) malloc(13*sizeof(float *)); for(i=0; i<13; ++i) { a[i] = (float *)malloc(4*sizeof(float)); }

Representation of 2D Array

Special Treatment of Array in NR • float *vector(long nl, long nh) allocate a float vector with index [nl..nh] • float **matrix(long nrl, long nrh, long ncl, long nch) allocate a 2D matrix with range [nrl..nrh] by [ncl..nch]

Header File in NR #include “nr.h” #include “nrutil.h”

Precedence and Association

Pre/post Increment, Address of • Consider f(++i) vs f(i++), what is the difference? • &a vs *a • Conditional expression x = (a < b) ? c : d;

Macros in C #define DEBUG #define PI 3.141592653 #define SQR(x) ((x)*(x))

Computer Representation of Numbers • Unsigned or two’s complement integers (e.g., char) 1000 0000 = 128 or -128 1000 0001 = 129 or -127 1000 0010 = 130 or -126 1000 0011 = 131 or -125 . . . 1111 1100 = 252 or -4 1111 1101 = 253 or -3 1111 1110 = 254 or -2 1111 1111 = 255 or -1 0000 0000 = 0 0000 0001 = 1 0000 0010 = 2 0000 0011 = 3 0000 0100 = 4 0000 0101 = 5 0000 0110 = 6 . . . 0111 1111 = 127

Real Numbers on Computer Example for β=2, p=3, emin= -1, emax=2 ε 0 0.5 1 2 3 4 5 6 7 ε is called machine epsilon.

Floating Point, sMBe-E, not IEEE

IEEE 754 Standard (32-bit) • The bit pattern represents If e = 0: (-1)sf 2-126 If 0<e<255: (-1)s (1+f) 2e-127 If e=255 and f = 0: +∞ or -∞ and f ≠ 0: NaN … … s b-1 b-2 b-23 e f = b-12-1 + b-22-2 + … + b-232-23

Error in Numerical Computation • Integer overflow • Round off error • E.g., adding a big number with a small number, subtracting two nearby numbers, etc • How does round off error accumulate? • Truncation error (i.e. discretization error) • The field of numerical analysis is to control truncation error

(Machine) Accuracy • Limited range for integers (char, int, long int, and long long int) • Limited precision in floating point. We define machine ε as such that the next representable floating point number is (1 + ε) after 1. ε 10-7 for float (32-bit) and  10-15 for double (64-bit)

Stability • An example of computing Φnwhere • We can compute either by Φn+1 = ΦnΦ or Φn+1 = Φn-1 – Φn Results are shown in stability.c program

Reading Materials • “Numerical Recipes”, Chap 1. • “What every computer scientist should know about floating-point arithmetic”. Can be downloaded from http://www.validlab.com/goldberg/paper.ps • “The C Programming Language”, Kernighan & Ritchie

Problems for Lecture 1 (C programming, representation of numbers in computer, error, accuracy and stability, assignment to be handed in next week) 1. (a) An array is declared as char *s[4] = {“this”, “that”, “we”, “!”}; What is the value of s[0][0], s[0][4], and s[2][1]? (b) If the array is to be passed to a function, how should it be used, i.e., the declaration of the function and use of the function in calling program? If the array is declared as char t[4][5] ; instead, then how should it be passed to a function? 2. (a) Study the IEEE 754 standard floating point representation for 32-bit single precision numbers (float in C) and write out the bit-pattern for the float numbers 0.0, 1.0, 0.1, and 1/3. (b) For the single precision floating point representation (32-bit number), what is the exact value of machine epsilon? What is the smallest possible number and largest possible number? 3. For the recursion relation: F n+1 =Fn-1 – Fn with F0 and F1 arbitrary, find the general solution Fn. Based on its solution, discuss why is it unstable for computing the power of golden mean Φ? (Hint: consider trial solution of the form Fn = Arn).

Numerical Recipes