nrlec1

Numerical Recipes: 

Numerical Recipes The Art of Scientific Computing

Program Organization: 

Program Organization

Use a Clear Style: 

Use a Clear Style 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); } } (A) (B) (C)

Control Structure in C - loop: 

Control Structure 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: 

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

Control Structure - break: 

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

1D Array in C: 

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 b is from 0 to 3, but bb is 1 to 4.

1D Array Argument Passing: 

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: 

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: 

Representation of 2D Array

Special Treatment of Array in NR: 

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: 

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

Precedence and Association: 

Precedence and Association

Macros in C: 

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

Computer Representation of Numbers: 

Computer Representation of Numbers Unsigned or two’s complement integers (e.g., char) 0000 0000 = 0 0000 0001 = 1 0000 0010 = 2 0000 0011 = 3 0000 0100 = 4 0000 0101 = 5 0000 0110 = 6 . . . 0111 1111 = 127 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

Real Numbers on Computer: 

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

Floating Point, sMBe-E : 

Floating Point, sMBe-E

IEEE 754 Standard (32-bit): 

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 e f = b-12-1 + b-22-2 + … + b-232-23 b-1 b-2 b-23

Error in Numerical Computation: 

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: 

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

Stability: 

Stability An example of computing Φn where We can compute either by Φn+1 = Φn Φ or Φn+1 = Φn-1 – Φn Results are shown in a simple program

Reading Materials: 

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

Slide23: 

Problems for Lecture 1 (C programming, representation of numbers in computer, error, accuracy and stability) 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][4] ; 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 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: use generating function, or rewrite as a matrix equation.)