mysimpson.m

来自「自己编写的matlab程序 数值分析里的程序」· M 代码 · 共 25 行

clc
a=1;b=6;n=100;
if mod(n,2)~=0
    return
else
    m=n/2;
    h=(b-a)/n;
%     fplot('fun(x)',[a b],'r*')
    Sm=h/3*(fun(a)+fun(b));
    for k=1:m
        x(2*k-1)=a+(2*k-1)*h;
        x(2*k)=a+2*k*h;
        Sm=Sm+h/3*(4*(fun(x(2*k-1)))+2*(fun(x(2*k))));        
    end
    Sm=Sm-h/3*2*fun(b)
%     Sm=h/3*(fun(a)-fun(b));
%     for k=1:m
%         x(2*k-1)=a+(2*k-1)*h;
%         x(2*k)=a+2*k*h;
%         Sm=Sm+h/3*(4*(fun(x(2*k-1)))+2*(fun(x(2*k))));        
%     end
%     Sm
    quadl('fun(x)',a,b)
end