1 / 31

Numerical Recipes

Numerical Recipes. The Art of Scientific Computing ( with applications in computational physics ). Computer Architecture. CPU. Memory. von Neumann architecture: both data and instructions are stored in memory. External Storage. Program Organization. int main(void) { … }

tnocera
Download Presentation

Numerical Recipes

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

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

  2. Computer Architecture CPU Memory von Neumann architecture: both data and instructions are stored in memory. External Storage

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

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

  5. Data Types #include <stdio.h> int main(void) { 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); }

  6. 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

  7. 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… • Large number of operators exists in C, such as ++i, --j, a+=b, &, !, &&, ||, ^, ?a:b, etc

  8. 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)

  9. 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);

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

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

  12. 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

  13. 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 • double b[4], *bb; bb=b-1; then valid range of index for b is from 0 to 3, but bb is 1 to 4

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

  15. 2D Array in C double m[3][5]; defines fixed size array. Example below defines dynamic 2D array. double **m; m = (double **) malloc(3*sizeof(double *)); for(i=0; i<3; ++i) { m[i] = (double *)malloc(5*sizeof(double)); }

  16. Representation of 2D Array, fixed size vs dynamic

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

  18. Call by Value vs Call by Reference • In a function call, a copy of the value is passed to a function argument. E.g.: void func(int a, int b[]){ a = a + 1; b[0] = b[0] + 1; } int main(void) { int i, A[5]; i = 1; A[0] = 1; func(i, A); i = ?; A[0] = ?; }

  19. Header Files in NR #include ”nr.h” #include ”nrutil.h” Distinction of file names in <….> and in ”… ”

  20. Precedence and Association

  21. 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;

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

  23. 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

  24. 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.

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

  26. 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

  27. 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

  28. (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)

  29. 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

  30. Reading Materials • “Numerical Recipes”, Chap 1. • “What every computer scientist should know about floating-point arithmetic”. Can be downloaded from http://perso.ens-lyon.fr/jean-michel.muller/goldberg.pdf • “The C Programming Language”, Kernighan & Ritchie

  31. Problems for Lecture 1 (C programming, representation of numbers in computer, error, accuracy and stability, assignment to be handed in next week thursday) 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 above 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 fixed size: 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 largest positive 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).

More Related