Calculation of Pi using polygons inscribed in a circle

Pi is the ratio of the perimeter of the polygon to twice the radius of the circle

> restart;

> with(plottools):
c := circle([1,1], 1, color=red):l1:=line([1,1],[2,1],color=blue):
l2:=line([1,1],[1,2],color=blue):l3:=line([1,2],[2,1],color=blue,linestyle=3):plots[display]({c,l1,l2,l3});

In the above, the dotted line has length r*sqrt(2) and hence for n=4 the approximation to Pi is 4*sqrt(2)/2 ~ 2.828

This can be generalized to n*sin(Pi/n ) for a polygon of n sides

If we let x=1/n , then the approximation for Pi becomes g(x)= sin(Pi*x)/x which approaches Pi in the limit x->0

> n:=16;dataset:=evalf(seq([1+cos(2*Pi*m/n),1+sin(2*Pi*m/n)],m=0..n));

> with(plots): d:=plot([dataset],color=blue,title="n=16"):plots[display]({c,d});

Suppose we are given the measured values for n=8,16,32,64. Can we predict the value of Pi?

This is an extrapolation problem. Let's make a simple assumption that the measured values approach Pi in powers of 1/n

Assume Pi ~ a0 + a1*x + a2*x^2 +a3*x^3 and then solve for a0,a1,a2,a3

> f:=(x,a0,a1,a2,a3)->a0+a1*x+a2*x^2+a3*x^3;

> sln:=solve({f(1/8,a0,a1,a2,a3)=3.061467,f(1/16,a0,a1,a2,a3)=3.121445,f(1/32,a0,a1,a2,a3)=3.136548,f(1/64,a0,a1,a2,a3)=3.140331},{a0,a1,a2,a3});

> assign(sln);

> with(plots):

> p1:=plot(f(x,a0,a1,a2,a3),x=0..0.15):

> data:=[1/8,3.061467],[1/16,3.121445],[1/32,3.136548],[1/64,3.140331]:

> p2:=plot([data],style=point,symbol=cross,color=blue,symbolsize=20):

> display({p1,p2});

This is the approximate value for Pi

> a0;

The error in a0 is

> evalf(a0-Pi);

However the exact expression for Pi is an even function of x

> g:=x->sin(Pi*x)/x;evalf(taylor(g(x),x,6));

> restart;

Assume Pi ~ a0 + a1*x ^2+ a2*x^4 +a3*x^6 and then solve for

a0,a1,a2,a3

> f:=(x,a0,a1,a2,a3)->a0+a1*x^2+a2*x^4+a3*x^6;

> sln:=solve({f(1/8,a0,a1,a2,a3)=3.061467,f(1/16,a0,a1,a2,a3)=3.121445,f

> (1/32,a0,a1,a2,a3)=3.136548,f(1/64,a0,a1,a2,a3)=3.140331},{a0,a1,a2,a3

> });

> assign(sln);

> with(plots):

Warning, the name changecoords has been redefined

> p1:=plot(f(x,a0,a1,a2,a3),x=0..0.15):

> data:=[1/8,3.061467],[1/16,3.121445],[1/32,3.136548],[1/64,3.140331]:

> p2:=plot([data],style=point,symbol=cross,color=blue,symbolsize=20):

> display({p1,p2});

This is a better approximate value for Pi

> a0;

The exact value is

> evalf(Pi);

>

The error in a0 is now

> evalf(a0-Pi);

>