数値積分のいろいろ

http://www.takuichi.net/hobby/edu/em/mom/numerical/index-j.htm
を計算する

リーマン積分（区分求積法）

function riemann

clear all;
close all;

divlist = [10e2,10e3,10e4,10e5,10e6,10e7];
errorlist = [];
for cycle=1:6
    
    div = divlist(cycle);
    x = linspace(0,1,div);
    fx = 4./(1+x.^2);

    dx = 1/div;
    result = 0;
    for idx=1:div
        result = result + fx(idx)*dx;
    end
    error = result - pi;
    errorlist(end+1) = error;
end

loglog(divlist,errorlist,'bo');
title('Riemann integration');
xlabel('number of division');
ylabel('error');

end

台形公式（直線補間）

function trapezoid

clear all;
close all;

divlist = [10e2, 10e3,10e4,10e5,10e6,10e7];
errorlist = [];
errorlist2 = [];
for cycle=1:length(divlist)
    
    div = divlist(cycle);
    x = linspace(0,1,div);
    fx = 4./(1+x.^2);

    dx = 1/div;

    result = 0;
    for idx=1:div
        result = result + fx(idx)*dx;
    end
    error = result - pi;
    errorlist(end+1) = error;
       
    result2 = 0;
    for idx=1:div-1
        result2 = result + 0.5*(fx(idx)+fx(idx+1))*dx;
    end
    error2 = result2 - pi;
    errorlist2(end+1) = error2;
end

loglog(divlist,errorlist,'r')
hold on,
legend({'Riemann'});
ax2=axes('xaxislocation','top','yaxislocation','right','color','none', 'xscale', 'log', 'yscale', 'log')
hold on
loglog(divlist,errorlist,'g')
legend({'Trapezoid'});

title('Trapezoid formula integration');
xlabel('number of division');
ylabel('error');

end

台形公式の方が、短冊より一桁精度がよくなる

シンプソンの公式

この問題に適用するのは間違ってる→通常3次以下の方程式に用いる
繰り返し計算を必要としない利点がある

function simpson

clear all;
close all;

divlist = [10e2, 10e3,10e4,10e5,10e6,10e7];
errorlist = [];
errorlist2 = [];
errorlist3 = [];

for cycle=1:length(divlist)
    
    div = divlist(cycle);
    x = linspace(0,1,div);
    fx = 4./(1+x.^2);

    dx = 1/div;

    result = 0;
    for idx=1:div
        result = result + fx(idx)*dx;
    end
    error = result - pi;
    errorlist(end+1) = error;
       
    result2 = 0;
    for idx=1:div-1
        result2 = result + 0.5*(fx(idx)+fx(idx+1))*dx;
    end
    error2 = result2 - pi;
    errorlist2(end+1) = error2;

    result3 = 1/6*(fx(1)+4*fx(div/2)+fx(end));
    
    error3 = result3 - pi;
    errorlist3(end+1) = error3;
end

loglog(divlist,errorlist,'r')
hold on,
legend({'Riemann'});
ax2=axes('xaxislocation','top','yaxislocation','right','color','none', 'xscale', 'log', 'yscale', 'log')
hold on
% loglog(divlist,errorlist2,'g')
% set(gca,'YDir','reverse'); 
% legend({'Trapezoid'});
% hold on
loglog(divlist,errorlist3,'b')
set(gca,'YDir','reverse'); 
legend({'Simpson'});

title('Trapezoid formula integration');
xlabel('number of division');
ylabel('error');


end

ガウス＝ルジャンドル積分

モンテカルロ積分