%% Composite numerical integration rule

%% 計算の手順：
%% 1) f.m　をダウンロードして、numerical_integration.mと同じフォルダの下に保存する。
%% 2) Freematの命令Windowで、numerical_integrationを打って、実行する。

clc; % 画面をクリアする。

index=1;

h_list = 1 ./ [10,20,40];
E_h_1 = zeros(3,1);

for n = [10,20,40]
    
    index
    
    h=1/n;    
    x=0:h:1;
    
    I_h_1 = 0;
    for k=1:n
        I_h_1 = I_h_1 + f( x(k) + h/2 )*h ;
    end
    
    E_h_1(index) =  abs( pi - I_h_1);  
    
    fprintf('Composite midpoint rule: I_h_%d = %17.15f \n', index, I_h_1 )
    fprintf('Composite midpoint rule: E_h_%d = %17.15f \n', index, E_h_1(index) )

    index=index+1;    
end    

plot( log(h_list), log(E_h_1),'+-')
grind on


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Gaussian Integration 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I_h = 0.5*( 5/9 * f(-sqrt(3/5)/2+0.5) + 8/9*f(0.5) + 5/9*f(sqrt(3/5)/2+0.5) ) ;
E_h = abs(pi - I_h);

fprintf('\n\nGaussian Integration: I_h = %17.15f \n', I_h )
fprintf('Gaussian Integration: E_h = %17.15f \n', E_h )
