===================================================
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.
~
~
~
~