Numerical Recipes

1 / 30

# Numerical Recipes - PowerPoint PPT Presentation

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 &lt;stdio.h&gt; main() {

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.

## PowerPoint Slideshow about ' Numerical Recipes' - keegan

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)

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

}

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]

#include “nr.h”

#include “nrutil.h”

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

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

• “Numerical Recipes”, Chap 1.

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