XMU《计算方法》实验三 龙贝格算法

实验三 龙贝格算法 实验报告

一、代码

clear;
fun = inline(input('请输入函数：f(x) = ', 's'));
a = input('请输入下界 a = ');
b = input('请输入上界 b = ');
e = input('请输入误差限 e = ');

h = b - a;
k = 1;
N = 1;
T(1, 1) = h / 2 * (fun(a) + fun(b));
while true
    sum = 0;
    for i = 1:N
        sum = sum + fun(a + (i - 1/2) * h);
    end
    T(k + 1, 1) = T(k, 1) / 2 + h / 2 * sum;
    for m = 1:k
        T(k - m + 1, m + 1) = (4^m * T(k + 1 - m + 1, m - 1 + 1) - T(k - m + 1, m - 1 + 1)) / (4^m - 1);
    end
    if abs(T(1, k + 1) - T(1, k - 1 + 1)) <= e
        break
    end
    h = h / 2;
    k = k + 1;
    N = N * 2;
end

format long;
disp('积分结果：')
I = T(1, k + 1)
disp('龙贝格矩阵 T:')
for i = 0:k
    for j = 0:i
        fprintf('%.15f ', T(i-j+1, j+1))
    end
    fprintf('\n');
end

二、程序运行截图

三、算法说明

龙贝格算法主要是运用理查德外推加速法来计算积分。

公式为：

\[T_m^{(k)} = \frac{4^m}{4^m-1} T_{m-1}^{(k+1)} - \frac{1}{4^m-1} T_{m-1}^{(k)} \]

计算流程如下：

1. 令 \(h = b-a, k=1\)，计算 \(T_{0, 0} = \frac h2 [f(a) + f(b)]\)。

2. 对于第 \(k\) 行，计算边界 \(T_{k, 0} = \frac 12[T_{k-1, 0} + h\sum\limits_{i=1}^{2^{k-1}} f(a + (i - \frac 12)h)]\)。

3. 枚举 \(m = 1, 2, \dots, k\)，计算

   \[T_{k-m, m} = \frac{4^mT_{k+1-m, m-1} - T_{k-m,m-1}}{4^m-1} \]

4. 若 \(|T_{0, k} - T_{0, k-1}|\leq \varepsilon\)，则说明达到误差限，停止计算，答案为 \(T_{0, k}\)。

5. 否则，令 \(k = k+1, h = \frac h2\)，从第 2 步再次计算。

四、心得体会

龙贝格求积算法是一种数值积分方法，可以有效地计算定积分的近似值。在本次实验中，我使用 Matlab 编写了龙贝格求积算法，并对其进行了测试和优化。

通过编写和测试龙贝格求积算法，我深刻认识到了数值积分方法的重要性。在实际应用中，我们经常需要对复杂函数进行积分，但往往难以找到解析解。此时，数值积分方法可以有效地解决这个问题，为实际问题的求解提供了很大的便利。

通过实验，我进一步加深了对龙贝格求积算法的理解。龙贝格求积算法是一种自适应方法，通过不断加细网格、计算更高阶的近似值，从而提高积分的精度。这种方法的优点是可以自动适应函数的特性，精度相对较高，但计算复杂度也相应较高。

通过对龙贝格求积算法进行优化，我发现在实际应用中，合理选择步长和网格数量是提高算法效率和精度的关键。同时，在处理一些特殊函数时，如存在间断点或振荡的函数，需要采用适当的技巧来处理，以确保算法的稳定性和收敛性。

通过本次实验，我不仅学会了一种新的数值积分方法，更加深了对数值计算的理解，并在实践中获得了不少经验和技巧。