Octave实现Romberg方法的T₀ₖ计算远慢于Python的原因排查

我在Octave里实现Romberg积分时，发现矩阵首行（T₀ₖ方法）的计算速度比自己写的Python版本慢很多。我用了数组而非字典，也没用到复杂数据结构，怀疑问题出在这段代码里：

for i = 0:2^k
    if i == 0 || i == 2^k
        sum = sum + f(x_i(a, i, hk)) / 2;
    else
        sum = sum + f(x_i(a, i, hk));
    end
end

Python实现代码 #

import math

def h_k(a, b, k):
    return (b - a)/2**k

def x_i(a, i, h_k):
    return a + i * h_k

def T_0k(k, a, b, f):
    h = h_k(a, b, k)
    sum = 0
    for i in range(2**k + 1):
        if i == 0 or i == 2**k:
            sum += f(x_i(a, i, h))/2
        else:
            sum += f(x_i(a, i, h))

    return h * sum

def T(m, k, d):
    d[(m, k)] = (4**m * d[(m-1, k+1)] - d[(m-1, k)])/(4**m - 1)

def Romberg_tab(a, b, f, m):
    d = {}
    for i in range(m+1):
        d[(0, i)] = T_0k(i, a, b, f)
    for i in range(1, m + 1):
        for j in range(m - i + 1):
            T(i, j, d)

    for m in range(21):
        print(f"T({m},{20 - m}): ", d[(m, 20 - m)])

a = -3
b = 2
def f(x):
    return 2024*x**8 - 1977*x**4 - 1981

Romberg_tab(a, b, f, 20)

Octave实现代码 #

disp("result is: ")

function hk = h_k(a, b, k)
    hk = (b - a) / 2^k;
end

function xi = x_i(a, i, hk)
    xi = a + i * hk;
end

function T_0k_result = T_0k(k, a, b, f)
    disp("T_0k_result")

    hk = h_k(a, b, k);
    sum = 0;
    for i = 0:2^k
        if i == 0 || i == 2^k
            sum = sum + f(x_i(a, i, hk)) / 2;
        else
            sum = sum + f(x_i(a, i, hk));
        end
    end
    T_0k_result = hk * sum;
end

function T_result = T(m, k, d)
    disp("T_result")

    T_result = (4^m * d{m - 1, k + 1} - d{m - 1, k}) / (4^m - 1);
end

function d = Romberg_table(a, b, f, m)
    disp("Romberg_table")
    d = cell(m+1, m+1);
    for i = 1:m+1
        d{1, i} = T_0k(i, a, b, f);
    end

    for i = 2:m+1
        for j = 1:m+1-i
            d{i, j} = T(i, j, d);
        end
    end
end

function result = f(x)
    result = 2024 * x^8 - 1977 * x^4 - 1981;
end

a = -3;
b = 2;
m = 20;

d = Romberg_table(a, b, @f, m);

核心原因分析 #

1. 函数调用开销过大：Octave的函数调用（尤其是x_i和f的频繁调用）开销远高于Python。当k=20时，循环次数达到100多万次，每次迭代都调用两个函数，累计的开销会被无限放大。而Python的函数调用经过多年优化，在这类简单计算场景下反而更高效。
2. 解释型循环效率劣势：Octave的循环是解释执行的，天生比Python的循环慢（Python的循环底层有CPython的优化加持）。百万级循环次数下，这个差异会非常明显。
3. 重复计算与分支判断：循环条件和判断里重复计算2^k，虽然Octave可能会做优化，但还是会额外消耗资源；同时每次循环都做分支判断，也会增加不必要的开销。
4. 冗余IO输出：代码里频繁的disp输出会占用大量时间，IO操作的速度远慢于计算操作。

优化方案 #

1. 内联函数调用，减少开销 #

把x_i的计算直接写到循环里，同时把首尾的判断移出循环，减少分支次数：

function T_0k_result = T_0k(k, a, b, f)
    hk = (b - a) / 2^k;
    n = 2^k;
    sum = f(a) / 2;  # 直接处理首项
    for i = 1:n-1
        sum = sum + f(a + i * hk);
    end
    sum = sum + f(b) / 2;  # 直接处理末项
    T_0k_result = hk * sum;
end

2. 用向量化操作替代循环（Octave的核心优势） #

Octave擅长向量化计算，底层是C实现，速度比解释型循环快几个数量级：

function T_0k_result = T_0k(k, a, b, f)
    hk = (b - a) / 2^k;
    n = 2^k;
    x = linspace(a, b, n+1);  # 一次性生成所有采样点
    weights = ones(1, n+1);
    weights([1, end]) = 0.5;  # 设置首尾权重
    sum_val = sum(weights .* f(x));  # 向量化加权求和
    T_0k_result = hk * sum_val;
end

3. 预计算常量 #

把2^k这类重复计算的值提前算好，避免循环里重复计算：

n = 2^k;  # 只计算一次
for i = 0:n
    ...
end

4. 移除调试用的disp输出 #

所有disp("T_0k_result")这类调试输出在正式运行时都要去掉，减少IO开销。

内容的提问来源于stack exchange，提问作者Bialy_kaloryfer