% Approx av pi - praktiserar Richardson

NN=3*2.^[0:7];
%NN=[3 5 10 20];
%NN=[3 10 25 50 100 200 500 10.^[3:8]];


% Berakanr forst ett antal approximationer med olika N.

p=0;
for i=1:length(NN);
  N=NN(i);
  p(i)=N*sind(180/N);
end;
p

fel=p-pi;

% Tittar hur felet ser ut. Plottar pa olika satt och tittar.

%plot(NN,p,'-o')
%plot(NN,fel,'-o')
loglog(NN,abs(fel),'-o')

% Inser att felet blir rak linje i loglog, dvs ett polynom.

% Da hjalper Richardson!
% Raknar nu om med stegdubblering (dvs bytte NN ovan), 
% for att kunna Richardson-extrapolera:


format short

dp=diff(p);

kvotkoll= dp (1:(end-1)) ./ dp(2:end)


phatt=p(2:end)+dp/3;

d2=diff(phatt);

kvot2= d2(1:(end-1)) ./ d2(2:end)

pdubbhatt=phatt(2:end)+d2/15;

d3=diff(pdubbhatt);

kvot3= d3(1:(end-1)) ./ d3(2:end)

ptriphatt=pdubbhatt(2:end)+d3/63;

format long

rtab=[p', [0;phatt'],  [0;0;pdubbhatt'],  [0;0;0;ptriphatt']]



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1a. PLotta p vs N (3:16)
%1b. PLotta p vs N (3:50) - ännu fel redan i tredje decimalen!
%1c. PLotta p vs N (3:500) - ännu fel redan i femte decimalen!

%2 plotta loglog - blir rät linje! visar att felet är polynom i N. - Rich!!!

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Första beräkning, utan Richardson-extrapolation


%         N       p                   fel
%          3   2.598076211353316    -5.44 e-01
%         10   3.090169943749474    -5.14 e-02
%         25   3.133330839107606    -8.26 e-03
%         50   3.139525976465669    -2.07 e-03
%        100   3.141075907812829    -5.17 e-04
%        200   3.141463462364135    -1.29 e-04
%        500   3.141571982779475    -2.07 e-05
%       1000   3.141587485879564    -5.17 e-06
%      10000   3.141592601912665    -5.17 e-08
%     100000   3.141592653073022    -5.17 e-10
%    1000000   3.141592653584626    -5.17 e-12
%   10000000   3.141592653589741    -5.20 e-14
%  100000000   3.141592653589792    -8.88 e-16
  

  
% 
%          N       ptriphat            fel
%          24  3.141592650573243    -3.02 e-09
%          48  3.141592653577892    -1.19 e-11
%          96  3.141592653589746    -4.71 e-14
%         192  3.141592653589792    -8.88 e-16
%         384  3.141592653589793        0

% rtab =
%
% 2.598076211353316   0.000000000000000   0.000000000000000   0.000000000000000
% 3.000000000000000   3.133974596215561   0.000000000000000   0.000000000000000
% 3.105828541230249   3.141104721640332   3.141580063335316   0.000000000000000
% 3.132628613281238   3.141561970631567   3.141592453897650   3.141592650573243
% 3.139350203046867   3.141590732968743   3.141592650457888   3.141592653577892
% 3.141031950890509   3.141592533505057   3.141592653540811   3.141592653589746
% 3.141452472285462   3.141592646083779   3.141592653589027   3.141592653589792
% 3.141557607911857   3.141592653120656   3.141592653589781   3.141592653589793
