【数值分析实验MATLAB】数值积分：梯形公式、辛普森公式、复化梯形公式、复化辛普森公式、龙贝格算法、自适应求积方法

awesome.st

39493人浏览 · 2020-01-21 17:06:50

awesome.st · 2020-01-21 17:06:50 发布

（1）Matlab代码:
f1:

if x==0
    y=1;
else
    %y=exp(-x^2);
    y=(sin(x))/x;
end

梯形公式：

a=0;
b=1;
tx=(b-a)/2*(f1(a)+f1(b));
tx=vpa(tx,5);
tx

simpson公式：

a=0;
b=1;
xps=(b-a)/6*(f1(a)+4*f1((a+b)/2)+f1(b));
xps=vpa(xps,5);
xps

（2）Matlab代码：
复化梯形公式：

a=0;
b=1;
e=1e-6;
m=1;
h=b-a;
T(1)=h/2*(f1(a)+f1(b));
while 1
    h=h/2;
    sum=0;
    for j=1:2^(m-1)
        sum=sum+f1(a+(2*j-1)*h);
    end
    S(1)=T(1)/2+h*sum;
    for j=1:m
        S(j+1)=S(j)+(S(j)-T(j))/(4^j-1);
    end
    if(abs(S(m+1)-T(m))<=e)
        fprintf('求得最终解为：%f\n',S(m+1)); 
        fprintf('迭代次数为:%d\n',m);
        return
    end
    T=S;
    m=m+1;
end

复化辛普森公式：

a=0;
b=1;
e=1e-6;
N=20;
n=0;
h0=(b-a)/2;
S0=h0/3*(f1(a)+f1(b)+4*f1((a+b)/2));
while n<N
    n=n+1;
    h1=h0/2;
    sum1=0;
    sum2=0;
    for k=1:2^(n-1)
        sum1=sum1+f1(a+(4*k-2)*h1);
    end
    for k=1:2^n
        sum2=sum2+f1(a+(2*k-1)*h1);
    end
    S1=S0/2+h1/3*(-2*sum1+4*sum2);
    if(abs(S1-S0)<e)
        fprintf('求得最终解为：%f\n',S1); 
        fprintf('迭代次数为:%d\n',n);
        return
    end
      S0=S1;
      h0=h1;
end

龙贝格算法：

a=0;
b=1;
e=1e-6;
h=b-a;
T(1)=h.*(f1(a)+f1(b))./2;
m=1;
while(1)
    h=h/2;
    sum=0;
    for j=1:2^(m-1)
        sum=sum+f1(a+h*(2*j-1));
    end
    S(1)=T(1)/2+h*sum;
    for j=1:m
        S(j+1)=S(j)+(S(j)-T(j))/(4^j-1);
    end
    if(abs(S(m+1)-T(m))<=e)
        break;
    end
    T=S;
    m=m+1;
end
fprintf('求得解为：%f\n',S(m+1));
fprintf('迭代次数：%d\n',m);

自适应求积方法：

format short
a=0;
b=1;
e=1e-6;
p=[a,b];
p0=p;
ep=[e];
m=0;
q=0;
I=0;
while(1)
    n1=length(ep);
    n=length(p0);
    if n==1
        break;
    end
    h=p0(2)-p0(1);
    s0=h/6*(f1(p0(1))+4*f1(p0(1)+h/2)+f1(p0(1)+h));
    s1=h/12*(f1(p0(1))+4*f1(p0(1)+h/4)+f1(p0(1)+h/2));
    s2=h/12*(f1(p0(1)+h/2)+4*f1(p0(1)+h*3/4)+f1(p0(1)+h));
    if abs(s0-s1-s2)<=15*ep(1)
        I=I+s1+s2;
        p0=p0(2:n);
        if n1>=2;
            ep=ep(2:n1);
        end
        q=q+1;
    else
        m=m+1;
        p0=[p0(1),p0(1)+h/2,p0(2:n)];
        if n1==1
            ep=[ep(1)/2,ep(1)/2];
        else
            ep=[ep(1)/2,ep(1)/2,ep(2:n1)];
        end
        if q==0
            p=p0;
        else
            p=[p(1:q),p0];
        end
    end
end
fprintf('求得的解为：%f\n',I);
fprintf('分层数为%d\n',m);
p=p';
disp('分点为');p