- 117 Views
- Uploaded on

Download Presentation
## PowerPoint Slideshow about ' Numerical Recipes' - keegan

**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.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Presentation Transcript

### Numerical Recipes

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

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

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

}

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]

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;

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.

IEEE 754 Standard (32-bit)

- The bit pattern

represents

If e = 0: (-1)sf 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).

Download Presentation

Connecting to Server..