龙贝格算法（Matlab实现）

目录

​​1、原理​​

​​2、案例 ​​

​​3、代码 ​​

​​ 4、结果​​

1、原理

2、案例 

3、代码 

clear;clc;close all;
format long
b=2*pi;a=0;h=b-a;eps=10^(-6);%误差eps
kmax=10;%最大递推次数
T1=h*(a*sin(a)+b*sin(b))/2%代换公式
S1=0;C1=0;C2=0;R1=0;R2=0;
for k=1:kmax
    h=(b-a)/2^k;
    i=1:2^(k-1);
    x=a+(2*i-1)*h;
    fx=sum(x.*sin(x));%代换
    T2=T1/2+fx*h 
    S2=T2+(T2-T1)/3
    if(k<3)
        if k==2
            C2=S2+(S2-S1)/15
        end
    else
        C2=S2+(S2-S1)/15
        R2=C2+(C2-C1)/63
        if abs(R2-R1)<eps
            break;
        end
        R1=R2;
    end
    T1=T2;S1=S2;C1=C2;
end
fprintf('I=%9.8f\n',R2);
第二个积分：
clear;clc;close all;
format long
b=1;a=0;h=b-a;eps=10^(-6);%误差eps
kmax=10;%最大递推次数
T1=h*(exp(-a^2)+exp(-b^2))/2%代换公式
S1=0;C1=0;C2=0;R1=0;R2=0;
for k=1:kmax
    h=(b-a)/2^k;
    i=1:2^(k-1);
    x=a+(2*i-1)*h;
    fx=sum(exp(-x.^2));%代换
    T2=T1/2+fx*h 
    S2=T2+(T2-T1)/3
    if(k<3)
        if k==2
            C2=S2+(S2-S1)/15
        end
    else
        C2=S2+(S2-S1)/15
        R2=C2+(C2-C1)/63
        if abs(R2-R1)<eps
            break;
        end
        R1=R2;
    end
    T1=T2;S1=S2;C1=C2;
end
fprintf('I=%9.8f\n',2/pi^0.5*R2);
定积分计算：
syms x
format long
a=int(x*sin(x),0,2*pi);%定积分
i=vpa(a)%计算数值
 
b=int(exp(-x.^2),0,1);
c=2*b/pi^0.5;
j=vpa(c)

 4、结果