【Matlab】Romberg求积公式求积分

我不是女神ヾ 2022-12-29 00:25 164阅读 0赞

## 1 题目 ##

用Romberg求积法计算积分 ∫ − 1 1 1 1 + 100 x 2 d x \\int\_\{-1\}^\{1\} \\frac\{1\}\{1+100 x^\{2\}\} d x ∫−111\+100x21dx的近似值，要求误差不超过 0.5 × 1 0 − 7 0.5\\times10^\{-7\} 0.5×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 运行结果 ##

![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1NlZV9TdGFy_size_16_color_FFFFFF_t_70]

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L1NlZV9TdGFy_size_16_color_FFFFFF_t_70]: /images/20221120/219e7470c95e4d298c3297356de2e57d.png