=================================================== Math 451 Laboratory 2: Tools and Pi's Date: 2013-02-15 =================================================== So far in our course, we've focused on using Matlab (or octave) for our numerical calculations, primarily because of the convenience and the need for a common language. However, there are many tools available, some good, some bad. In this lab, we'll see some simple examples of other tools that are freely available and useful for calculations. Goals ---------------------------- - See a simple Monte-Carlo approach to calculating pi in Matlab. - See an infinite series calculation in Python - Learn how to create and run a python script - See a recursive calculation in bc - Learn how to do arbitrary-precision calculations in Linux with bc - See an optimal implementation in C - Learn how to create and compile a C program Methods of calculating Pi =============================================================================== Dartboard estimates for pi in Matlab ---------------------------------------- This method is what's called a Monte Carlo method -- we make a bunch of guesses, test the results and look for the average outcome as our estimate. (1-a) Run the following script in Matlab seq = [1,1,1]; for L=1:8 N = 10^L; a = rand(2,N); pi_est = sum(sum(a.*a)<1)*4/N; seq(end+1,1:3) = [N,pi_est,abs(pi_est-pi)/pi]; end figure(1); clf; loglog(seq(:,1),seq(:,3),'o-'); seq (1-b) Save your figure and hand it in with this assignment. (1-c) How many decimal places of accuracy do you predict with 10^9 points? Series with Python =============================================================================== Python is an interpreted like Matlab, but much more powerful and widely used for non-mathematical problems. The main website is "python.org" which contains lots of documentation. The dictator is Guido. If you recall from trigonometry, tan(pi/4) = 1. When we take the inverse, then, we can show pi = 4 arctan(1). Thus, if we can approximate the arctan function, we can approximate pi. This method is referred to as the Gregory-Leibniz series, in honor of it's discoverers in European history. The series is ... 3 5 7 9 x x x x / 10\ atan(x) = x - --- + --- - --- + --- + O\x / 3 5 7 9 (2-a) Open a text file on the desktop named "gregory_leibniz.py". Enter the following code. WARNING: python indentation is fundamental to the language. So make sure you do it exactly as shown below. #!/usr/bin/env python max_term = 200 pi = 3.1415926535897932384626433832795 def term(n,x): return (2 - n%4) * x**n / n pi_est = 0 for n in range(1,max_term,2): pi_est += term(n,1.) print n, 4*pi_est, abs(4*pi_est-pi) (2-b) Open a terminal window. Type the following 2 commands. cd Desktop python gregory_leibniz.py (2-c) What are the best upper and lower bounds on pi achieved by this method? (2-d) Is this method better or worse than the dart-board method above? Explain. (2-e) What does "(2-n%4)" do? (2-f) Would this script work better if we used tan(pi/3) or tan(pi/6)? Bryson's method with square roots in bc =============================================================================== Two of the first and most widely available calculator tools are "dc" standing for "desktop calculator" and "bc", it's "better" brother. "dc" is a reverse Polish notation calculator of great power, but since reverse Polish notation is less familiar notation to us, we will make use of "bc", which uses infix notation like we are used to. The big advance of these tools is that they can do ARBITRARY PRECISION calculations. Want 20 decimal places -- easy. Want 100 decimal places -- no problem! The only limitation is the size of you computer's memory. Unfortunately, "bc" and "dc" are baron of may common features of today's languages. For instance, in many versions, variables must be represented by a single letter. Every once and a while, you run into a problem where you need there raw power and you don't want to pay for a big, fancy tool. Here, we'll use bc to calculate pi using a method based on Bryson's method and the half-angle formulas for sine and cosine in trigonometry. Basically, we make smaller and smaller triangles inside of a circle with sides x by y, and then add the area's of all these triangles up. But to use half-angle formula's, we have to be able to calculate square roots efficiently, which was not needed in the previous methods. (3-a) Open a text file on the desktop named "bryson.bc". Enter the following code. #!/usr/bin/bc -lq # http://www.piday.org/million.php pi = 3.1415926535897932384626 scale = 16 # number of digits b = 1. n = 2 t = 25; # max number of iterations print "Iteration\tEstimate\t\tError\n" for (i=0;i a[52514],b,c=52514,d,e,f=1e4,g,h; main(){ for(;b=c-=14;h=printf("%04d",e+d/f)) for (e=d%=f;g=--b*2;d/=g)d=d*b+f*(h?a[b]:f/5),a[b]=d%--g; } (2-b) Now, we need to compile and run our code. Open a terminal window. Type the following 2 commands. cd Desktop gcc -O3 -o ramanujan ramanujan.c ramanujan (2-c) What's the 33rd digit of pi? (2-d) Anybody who can reconstruct the formula this code uses, and explain the implementation get's an A on this lab. Ancient tools and Latest fads ================================================================================ Fortran is king still in many circles. See http://www.netlib.org/ A newly released programming language that is stirring interest is "Julia". Might be a fad, might be the next great thing in numerical methods programming. http://julialang.org/downloads/ http://ee380.stanford.edu/cgi-bin/videologger.php?target=120229-ee380-300.asx Check it out! Maybe next year, we'll use this to teach the course instead of Matlab. ~ ~ ~ ~