自选一个可积函数在某个区间上求数值积分，分别利用复化梯形公式、复化Simpson公式、复化Cotes公式进行MATLAB编程运算求得数值积分值及误差，并分别绘制几何图形表达求解过程。最后，利用Romberg求积公式进行数值积分，并对比结果。

在MATLAB中，我们可以选择一个简单的连续可积函数，比如f(x) = sin(x)，在特定区间[a, b]上计算其数值积分。这里我将简要展示如何使用复化梯形公式（Trapezoidal Rule）、复化辛普森法则（Simpson's Rule）以及复化柯特斯法则（Cotes's Formula），并利用Romberg求积公式进行比较。

首先，让我们定义函数和区间：

a = 0;
b = pi;
f = @(x) sin(x);

然后，我们编写各个积分方法的函数：

1. **复化梯形公式**

function result_Trapezoid = trapz_int(f, a, b)
    n = 100; % 网格点数
    h = (b - a) / n;
    x = a + h*(0:n); % 样本点
    result_Trapezoid = h * (f(x(1)) + f(x(end))) / 2 + h * sum(f(x(2:end-1)));
end

2. **复化辛普森法则**

function result_Simpson = simpson_int(f, a, b)
    n = 100; % 网格点数
    h = (b - a) / (n + 1);
    x = a + h*(0:2:n); % 样本点（奇数位置）
    result_Simpson = h / 6 * (4 * f(x(1:2:end)) + f(x(2:2:end-1)) + 2 * f(x(3:2:end)));
end

3. **复化柯特斯法则**

function result_Cotes = cotes_int(f, a, b)
    n = 100; % 网格点数
    h = (b - a) / n;
    % 实现具体的Cotes公式需要查阅相关资料，因为这个公式的具体形式取决于阶数
end

对于Romberg求积法，它是一个更高级的迭代积分算法，可以在给定精度下自动调整步长：

function result_Romberg = romberg_int(f, a, b, tol)
    R = @(f, a, b, N) (b - a) / N * sum((f(a + (0:N-1)*h) + f(b - (0:N-1)*h)) ./ (2*N + 1));
    h = (b - a) / 2^(0:8); % 初始步长序列
    error = inf; % 初始化误差
    for i = 1:length(h)
        if abs(R(f, a, b, i+1) - R(f, a, b, i)) < tol
            break;
        end
        error = R(f, a, b, i+1);
    end
    result_Romberg = R(f, a, b, i);
end

现在你可以按照以下步骤运行这些函数，获取积分值和误差，绘制曲线，并对比结果：

% 计算每个方法的结果
integral_Trapezoid = trapz_int(f, a, b);
integral_Simpson = simpson_int(f, a, b);
integral_Cotes = cotes_int(f, a, b);
integral_Romberg = romberg_int(f, a, b, 1e-8);

% 绘制积分误差图
errors = [integral_Trapezoid - quadl(@(x) f(x), a, b), integral_Simpson - quadl(@(x) f(x), a, b), integral_Cotes - quadl(@(x) f(x), a, b)];
plot(h, errors, 'o-', 'LineWidth', 2);
xlabel('Mesh Size');
ylabel('Error from True Integral');

% 输出结果和误差
disp(['Trapezoidal: ' num2str(integral_Trapezoid) ', Error: ' num2str(abs(integral_Trapezoid - quadl(@(x) f(x), a, b)))]);
disp(['Simpson: ' num2str(integral_Simpson) ', Error: ' num2str(abs(integral_Simpson - quadl(@(x) f(x), a, b)))]);
disp(['Cotes: ' num2str(integral_Cotes) ', Error: ' num2str(abs(integral_Cotes - quadl(@(x) f(x), a, b)))]);
disp(['Romberg: ' num2str(integral_Romberg) ', Error: ' num2str(abs(integral_Romberg - quadl(@(x) f(x), a, b)))]);