【Matlab】Romberg求积公式求积分

1 题目

用Romberg求积法计算积分${\int }_{-1}^1\frac{1}{1+100x^2}dx$的近似值，要求误差不超过$0.5\times 10^{-7}$。

2 程序

程序1：Romberg通用程序

function[result] = romberg(a,b,F,epsilon)
    %% 初始赋值
    e = inf;
    xk = [a,b];% 
    T = [];  % 复化梯形公式
    S = [];  % 复化Simpson公式
    C = [];  % 复化Cotes公式
    R = [];  % Romberg公式
    k = 1;   % 当前次数

    %% 计算
    digits(10)
    h = b - a;
    T1 = vpa( h / 2 * ( F(a) + F(b) ) );
    T = [T,T1];
    fprintf('k=%d, 2^k=%d, T=%0.8f\n',k,2^k,T1) 

    while(e>epsilon)
        xk2= [];
        for i=1:length(xk)-1
            xk2 = [xk2, (xk(i) + xk(i+1)) / 2];
        end

        % 计算Tn
        sum = 0;
        for i=1:length(xk2)
            sum = vpa(sum + h * F(xk2(i)));
        end
        T_temp = 0.5 * ( T(end) + sum);
        T = [T,T_temp];
        k = k + 1;
        h = h / 2;
        xk = [xk, xk2];
        xk = sort(xk);
        fprintf('k=%d, 2^k=%d, T=%0.8f',k,2^k,T_temp) 

        % 计算Sn
        if k>=2
            S_temp = 4 / 3 * T(end) - 1 / 3 * T(end-1);
            S = [S, S_temp];
            fprintf(', S=%0.8f',S_temp) 
        end

        % 计算Cn
        if k>=3
            C_temp = 16 / 15 * S(end) - 1 / 15 * S(end-1);
            C = [C, C_temp];
            fprintf(', C=%0.8f',C_temp) 
        end

        % 计算Rn
        if k>=4
            R_temp = 64 / 63 * C(end) - 1 / 63 * C(end-1);
            R = [R, R_temp];
            fprintf(', R=%0.8f',R_temp) 
        end

        % 计算误差e
        if k>=5
            e = abs(R(end) -  R(end-1));
            fprintf(', e=%0.8f',e) 
        end
        fprintf('\n') 
    end
    fprintf('Romberg求积公式的结果为：%0.8f\n',R(end))
    result = R(end);
end

程序2：主函数

clc;clear all

%% 初始数据
syms x
F(x) = 1 / (1 + 100 * x^2);
a = -1;
b = 1;
epsilon = 0.5e-7;

result = romberg(a,b,F,epsilon);

3 运行结果