#Header Begins--Do not remove################################################## # # Maple code to compute pi a multitude of different ways. # # (c) Eugene Fiume, University of Toronto, 1995 (modified 1996). # # You are permitted to use this code for noncommercial, educational # purposes only. This header must not be removed from this file. # #Header Ends--Do not remove#################################################### # # Leibnitz' rule to get for pi/4. # L := proc (n) local i; sum( (-1)^i / (2*i+1), i=0..n ); end: # # Sequence to compute Pi by Leibnitz' rule. # P1 := proc (n) local i; 4*sum( (-1)^i / (2*i+1), i=0..n ); end: # # Sequence to compute Pi by Euler's rule. # P2 := proc (n) local i; sqrt(6*sum( 1 / i^2, i=1..n )); end: # # Compute pi by perimeter of regular 2^n-gon inscribing a circle of # unit diameter. Thus P(2) is 2*sqrt(2). (Archimedes' rule.) # P3 := proc (n) option remember; if (n<2) then 0 elif (n=2) then 2*sqrt(2) else 2^(n-1)*sqrt(2*(1-sqrt(1-(P3(n-1)/2^(n-1))^2))); fi; end: # # Compute pi by perimeter of regular n-gon inscribing a circle of # unit diameter. This is the nonrecursive form that uses Pi to # compute it! # P4 := proc (n) local i; n * sqrt( 1/2*(1-cos(2*Pi/n)) ); end: # # The recursive form of P4 to avoid use of Pi. See the notes. # P4r := proc (n) local i; 2^n*sqrt(1/2*(1-C(n-1))); end: C := proc (n) option remember; if (n<2) then 0 elif (n=2) then sqrt(2)/2 else sqrt( 1/2*(C(n-1)+1) ); fi; end: # # This evaluates any of the above schemes to a floating point # representation. It is useful to observe how sensitive the # above schemes are to floating point errors. # Pe := proc (P,n) evalf(P(n)); end: # # Compute a sequence of values of one of the sequences above suitable # for plotting or computing errors. # compP := proc(P,n) local i; [seq( [i, P(i)], i=4..n)]; end: # # Plot a sequence of approximations to pi using one of the above. # plotP := proc(P,n) plot(compP(P,n)); end: # # The routines below compute various forms of error between the # sequence approximation and Pi. # errP := proc(P,n) local v,p; p := compP(P,n); v := Pi; plot([ [i, v-p[i][2]]$i=1..nops(p)]); end: # # Absolute error for approximation of procedure P to value v. # errPabs := proc(P,v,n) local p; p := compP(P,n); plot([ [i, abs(v-p[i][2])]$i=1..nops(p)]); end: # # Relative error for approximation of procedure P to value v. # errPrel := proc(P,v,n) local p; p := compP(P,n); plot([ [i, abs(v-p[i][2])/v]$i=1..nops(p)]); end: # # Do composite absolute error for some of the above routines. # (This is mainly just an example of how to superimpose several # plots into one display.) You can also use the "display" # procedure to do this. # errPabsAll := proc(N) local p1,p2,p3,p4,p5, e1,e2,e3,e4,e5, n; option remember; p1 := compP(P4r,N); p4 := compP(P4,N); e1 := [ [i, abs(Pi-p1[i][2])]$i=1..nops(p1)]; e4 := [ [i, abs(Pi-p4[i][2])]$i=1..nops(p4)]; plot({e1,e4},n=0..nops(p1)); end: # # Sample executions: # # P1(10); # (or P2, or P3, or ...) # # Pe(P1,100); # # Compare the following: # # evalf(P4r(10)); # P3(P4r,10); # # Why are these results different? # # plotP(P4r,10); # errPrel(P4r,15); # explain the results!!