Numerical recipes
Download
1 / 30

Numerical Recipes - PowerPoint PPT Presentation


  • 119 Views
  • Uploaded on

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() {

loader
I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
capcha
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

Numerical Recipes

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


Computer architecture
Computer Architecture

CPU

Memory

External Storage


Program organization
Program Organization

int main() {

}

double func(double x) {

}


First example
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
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
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

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

}


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

}



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”



Pre post increment address of
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
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)

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
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 b e e not ieee
Floating Point, sMBe-E, not IEEE


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

b-1

b-2

b-23

e

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


Numerical recipes

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