Scientific Computing Numerical Differentiation Dr Guy TelZur Clouds

Scientific Computing Numerical Differentiation Dr. Guy Tel-Zur Clouds. Picture by Peter Griffin, publicdomainpictures. net

Some Recent New Articles and Blog Posts

“MY SLICE OF PIZZA” Blog • ” - DIMACS Parallelism 2020: John Gustafson • “Throw out old Numerical Analysis textbooks! Algorithm designers have historically "measured" algorithm run times by counting the number of floating point operations / additions / multiplications. This made sense decades ago, when floating point arithmetic took 100 times as long as reading a word from memory. Now, one multiplication takes about 1. 3 nanoseconds (to go through the entire pipeline; this underestimates throughput), compared to 50 -100 nanoseconds for the memory access. Why do our algorithms measure the wrong thing? We should be counting memory accesses; it isn't reasonable to ignore the constant factor of 50. ” • Slides are here

Computing in Science and Engineering • March/April 2011 – A few papers about Python! • Python for Scientists and Engineers • Python: An Ecosystem for Scientific Computing • Mayavi: 3 D Visualization of Scientific Data • Cython: The Best of Both Worlds • The Num. Py Array: A Structure for Efficient Numerical Computation

Code examples from “Python for Scientists and Engineers“ Another example (Guy): Demo: symbolicmath. py from sympy import exp f=lambda x: exp(-x**2) integrate(f(x), (x, -inf, inf)) See also: “From Equations to Code: Automated Scientific Computing” By Andy R. Terrel”, Computing in Science and Engineering, March-April 2011

Next Slides from Michael T. Heath – Scientific Computing • Source: http: //www. cse. illinois. edu/heath/scicomp/n otes/chap 01. pdf

MHJ Chapter 3: Numerical Differentiation Should be f’c “ 2” stands for two points

forward/backward 1 st derivative: ±h f(x)=a+bx^2 Exact f’ f’=b f’=2 bx f 2’=((a+b(x+h))-(a+bx))/h=b f 2’=((a+b(x+h)^2)-(a+bx^2))/h = (a+bx^2+2 bxh+bh^2)-(a+bx^2))/h= 2 bx+bh Bad!!! Good! (1/2)(f 2 L’+f 2 R’)= f 2’=(f(x+h)-f(x-h))/(2 h) ((a+b(x+h)^2)-(a+b(x-h)^2))/(2 h)= (a+bx^2+2 bxh+bh^2 -a-bx^2+2 bxhbh^2)/(2 h)=4 bxh/2 h=2 bx Good!

N-points stencil

Example code • 2 nd derivative of exp(x) • Code in C++, we will learn more of the language features: – Pointers – Call by Value/Reference – Dynamic memory allocation – Files (I/O)

Call by value vs. call by reference • printf(“speed= %dn”, v); // this is a call by value as the value of v won’t be changed by the function (printf) – which is desired • scanf(“%dn”, &v); // this is a call by reference, we want to supply the address of v in order to set it’s value(s)

Analyzing an example code // This program module demonstrates memory allocation and data transfer // in between functions in C++ #include<stdio. h> #include<stdlib. h> int main(int argc, char *argv[]) int a; // line 1 int *b; // line 2 a = 10; // line 3 b = new int[10]; // line 4 for(i = 0; i < 10; i++) { b[i] = i; // line 5 } func( a, b); // line 6 return 0; } // End: function main() { void func( int x, int *y) // line 7 { x += 7; // line 8 *y += 10; // line 9 y[6] += 10; // line 10 return; // line 11 } // End: function func()

// // This program module // demonstrates memory allocation and data transfer in between functions in C++ #include<stdio. h> #include<stdlib. h> void func( int x, int *y); int main(int argc, char *argv[]) { int a; // line 1 int *b; // line 2 a = 10; // line 3 b = new int[10]; // line 4 for(int i = 0; i < 10; i++) { b[i] = i; // line 5 } func( a, b); // line 6 printf("a=%dn", a); printf("b[0]=%dn", b[0]); printf("b[1]=%dn", b[1]); printf("b[2]=%dn", b[2]); printf("b[3]=%dn", b[3]); printf("b[4]=%dn", b[4]); printf("b[5]=%dn", b[5]); printf("b[6]=%dn", b[6]); printf("b[7]=%dn", b[7]); printf("b[8]=%dn", b[8]); printf("b[9]=%dn", b[9]); return 0; } // End: function main() void func( int x, int *y) // line 7 { x += 7; // line 8 *y += 10; // line 9 y[6] += 10; // line 10 return; // line 11 } // End: function func() תכנית משופרת Check program: demo 1. cpp Under /lectures/02/code

Topics from MHJ Chapter 3 • Program 1: 2 nd derivative of exp(x) in C++ • Program 2: Working with files • Program 3: Same as program #1 plus working with files • Program 4: The same in Fortran 90 • Error Estimation • Plotting the error using gnuplot

! • A reminder for myself: open Dev. C++ for the demos!) • I slightly modified “program 1. cpp” from MHJ section 3. 2. 1 (2009 Fall edition): http: //www. fys. uio. no/compphys/cp/program s/FYS 3150/chapter 03/cpp/program 1. cpp • Install Dev. C++ on you personal laptop/computer and try it!

Explain program 1. cpp Working directory: c: UserstelzurDocumentsWeizmannScientific. ComputingSC 2011 BLectures2code> Open Dev. C++ IDE for the demo Program description: // Program to compute the second derivative of exp(x). // Three calling functions are included // in this version. In one function we read in the data from screen, // the next function computes the second derivative // while the last function prints out data to screen. // The variable h is the step size. We also fix the total number // of divisions by 2 of h. The total number of steps is read from // screen Usage: > program 1 <user input> 0. 01 10 100 Examine the output: > type out. dat

program 2. cpp • The book mentions program 2. cpp which is in cpp and the URL is indeed a cpp code, but the listing below the URL is in C. • This demonstrates the I/O differences between C and C++

using namespace std; #include <iostream> int main(int argc, char *argv[]) { FILE *in_file, *out_file; Working with files in C++ program 2. cpp if( argc < 3) { printf("The programs has the following structure : n"); printf("write in the name of the input and output files n"); exit(0); } in_file = fopen( argv[1], "r"); // returns pointer to the input file if( in_file == NULL ) { // NULL means that the file is missing printf("Can't find the input file %sn", argv[1]); exit(0); } out_file = fopen( argv[2], "w"); // returns a pointer to the output file if( out_file == NULL ) { // can't find the file printf("Can't find the output file%sn", argv[2]); exit(0); } fclose(in_file); fclose(out_file); return 0; }

program 3. cpp • Usage: > program 3 outfile_name • All the rest is like in program 1. cpp

Now lets check the f 90 version • Open in Sci. TE program 1. f 90 • In the image below: compilation and execution demo:

MHJ section 3. 2. 2 Error analysis

Content of exp(10)’’ computation • See MHJ section 3. 2. 2 and Fig. 3. 2 (Fall 2009 Edition) • Text output with 4 columns: h, computed_derivative, log(h), ε >program 1 Input: 0. 1 10. 10 >more out. dat 1. 00000 E-001 2. 72055 E+000 -1. 00000 E+000 -3. 07904 E+000 5. 00000 E-002 2. 71885 E+000 -1. 30103 E+000 -3. 68121 E+000 2. 50000 E-002 2. 71842 E+000 -1. 60206 E+000 -4. 28329 E+000 1. 25000 E-002 2. 71832 E+000 -1. 90309 E+000 -4. 88536 E+000 6. 25000 E-003 2. 71829 E+000 -2. 20412 E+000 -5. 48742 E+000 3. 12500 E-003 2. 71828 E+000 -2. 50515 E+000 -6. 08948 E+000 1. 56250 E-003 2. 71828 E+000 -2. 80618 E+000 -6. 69162 E+000 7. 81250 E-004 2. 71828 E+000 -3. 10721 E+000 -7. 29433 E+000 3. 90625 E-004 2. 71828 E+000 -3. 40824 E+000 -7. 89329 E+000 1. 95313 E-004 2. 71828 E+000 -3. 70927 E+000 -8. 44284 E+000

Download and install Gnuplot http: //www. gnuplot. info/

Visualization: Gnuplot • Reconstruct result from MHJ - Figure 3. 2 using gnuplot • Gnuplot is included in Python(x, y) package! • Gnuplot tutorial: http: //www. duke. edu/~hpgavin/gnuplot. html • Example: http: //www. physics. ohiostate. edu/~ntg/780/handouts/gnuplot_quadeq_example. pdf • 3 D Examples: http: //www. physics. ohio- state. edu/~ntg/780/handouts/gnuplot_3 d_example_v 2. pdf

Using gnuplot

Can we explain this behavior? Computed for x=10

Error Analysis εro = Round-Off error The approximation error: Recall Eq. 3. 4: The leading term in the error (j=1) is therefore:

The Round-Off Error (εro ) • εro depends on the precision level of the chosen variables (single or double precision) Single precision Double precision If the terms are very close to each other, the difference is at the level of the round off error

hmin = 10 -4 is therefore the step size that gives the minimal error in our case. If h>hmin the round-off error term will dominate

Summary • Next 3 slides are from: Michael T. Heath Scientific Computing • http: //www. cse. illinois. edu/heath/scicomp/n otes/chap 08. pdf

Let’s upgrade our visualization skills! • Mayavi • Included in the Python(x, y) package • 2 D/3 D • User guide: http: //code. enthought. com/projects/mayavi/ docs/development/html/mayavi/index. html

http: //code. enthought. com/projects/mayavi/

Mayavi is a general purpose, open source 3 D scientific visualization package that is tightly integrated with the rich ecosystem of Python scientific packages. Mayavi provides a continuum of tools for developing scientific applications, ranging from interactive and script-based data visualization in Python to fullblown custom end-user applications.

Computing in Science and Engineering, March-April 2011

Mayavi demos

from enthought. mayavi import mlab from numpy import ogrid x, y, z = ogrid[-10: 10: 100 j, -10: 100 j] ctr = mlab. contour 3 d(0. 5*x**2 + y**2 + 2*z**2) mlab. show() To execute: run cont. py

vector. py import numpy from enthought. mayavi. mlab import * x, y, z = numpy. mgrid[-2: 3, -2: 3] r = numpy. sqrt(x**2 + y**2 + z**4) u = y*numpy. sin(r)/(r+0. 001) v = -x*numpy. sin(r)/(r+0. 001) w = numpy. zeros_like(z) quiver 3 d(x, y, z, u, v, w, line_width=3, scale_factor=1)

import numpy from enthought. mayavi. mlab import * surface. py """Test surf on regularly spaced co-ordinates like Maya. Vi. """ def f(x, y): sin, cos = numpy. sin, numpy. cos return sin(x+y) + sin(2*x - y) + cos(3*x+4*y) x, y = numpy. mgrid[-7. : 7. 05: 0. 1, -5. : 5. 05: 0. 05] surf(x, y, f)

# Create the data. from numpy import pi, sin, cos, mgrid dphi, dtheta = pi/250. 0, pi/250. 0 [phi, theta] = mgrid[0: pi+dphi*1. 5: dphi, 0: 2*pi+dtheta *1. 5: dtheta] m 0 = 4; m 1 = 3; m 2 = 2; m 3 = 3; m 4 = 6; m 5 = 2; m 6 = 6; m 7 = 4; r = sin(m 0*phi)**m 1 + cos(m 2*phi)**m 3 + sin(m 4*theta)**m 5 + cos(m 6*theta)**m 7 x = r*sin(phi)*cos(theta) y = r*cos(phi) z = r*sin(phi)*sin(theta) # View it. from enthought. mayavi import mlab s = mlab. mesh(x, y, z) mlab. show() nice_mesh. py • http: //code. enthought. com/projects/mayavi/docs/d evelopment/html/mayavi/mlab. html#id 5

Mayavi environment

Move the object and zoom with the mouse

In one of the next lectures we will also learn visualization with Vis. It