Octave 中的复合梯形规则不是 运行
Composite trapezoid rule not running in Octave
我在 Octave 中有以下代码用于实现复合梯形规则，出于某种原因，该函数只会在我在 Octave 中执行时停顿 f = @(x) x^2, a = 0, b = 4 , 公差 = 10^-6。每当我调用 trapezoid(f, a, b, TOL) 时，什么都没有发生，我必须退出终端才能在 Octave 中执行任何其他操作。这是代码：
% INPUTS
%
% f : a function
% a : starting point
% b : endpoint
% TOL : tolerance
function root = trapezoid(f, a, b, TOL)
disp('test');
max_iterations = 10000;
disp(max_iterations);
count = 1;
disp(count);
initial = (b-a)*(f(b) + f(a))/2;
while count < max_iterations
disp(initial);
trap_0 = initial;
trap_1 = 0;
trap_1_midpoints = a:(0.5^count):b;
for i = 1:(length(trap_1_midpoints)-1)
trap_1 = trap_1 + (trap_1_midpoints(i+1) - trap_1_midpoints(i))*(f(trap_1_midpoints(i+1) + f(trap_1_midpoints(i))))/2;
endfor
if abs(trap_0 - trap_1) < TOL
root = trap_1;
return;
endif
intial = trap_1;
count = count + 1;
disp(count);
endwhile
disp(['Process ended after ' num2str(max_iterations), ' iterations.']);
我在 Matlab 中试过你的函数。
您的代码没有停止。而是 trap_1_midpoints 的大小呈指数增长。随着 trap_1 的计算时间也呈指数增长。这就是你所经历的拖延。
我还在您的代码中发现了一个可能的错误。我猜 if 子句后面的行应该是 initial = trap_1。检查缺失的 'i'.
这样，您的代码仍然需要很长时间，但是如果您增加容差（例如，将值设置为 1），您的代码 return.
您可以尝试矢量化 for 循环以加快速度。
编辑：我认为在您的 for 循环中，f(trap_1_midpoints(i+1) 之后缺少 )。
在 count=52 左右之后，算术序列 trap_1_midpoints 不再以任何有意义的方式用浮点数表示。在 count=1075 或类似的之后，步长不再表示为正浮点数 double 数。这就是说，绑定 max_iterations = 10000 是荒谬的。如下所述，count=20 之后的所有计算都没有意义。
步长 h 的理论误差为 O(T·h^2)。在 O(T/h) 个数字的总和中存在数值误差累积，即 O(mu/h) 与 mu=1ulp=2^(-52)。总的来说，对于双数 h=1e-5 或算法 count=17，数值积分的最低误差预计在 h=mu^(1/3) 左右。这可能会随着区间长度以及函数的平滑度或波浪度而变化。
人们可以预期错误除以四的行为，同时仅对于高于此边界的步长将步长减半 1e-5。这也意味着 abs(trap_0 - trap_1) 是仅在此步长范围内的 trap_0（和 abs(trap_0 - trap_1)/3 对于 trap_1）误差的可靠度量。
大约 h=1e-3 应满足误差界限 TOL=1e-6，这对应于 count=10。如果递归没有停止 count = 14（应该给出小于 1e-8 的错误），则该方法未准确实现。
我在 Octave 中有以下代码用于实现复合梯形规则，出于某种原因，该函数只会在我在 Octave 中执行时停顿 f = @(x) x^2, a = 0, b = 4 , 公差 = 10^-6。每当我调用 trapezoid(f, a, b, TOL) 时，什么都没有发生，我必须退出终端才能在 Octave 中执行任何其他操作。这是代码：
% INPUTS
%
% f : a function
% a : starting point
% b : endpoint
% TOL : tolerance
function root = trapezoid(f, a, b, TOL)
disp('test');
max_iterations = 10000;
disp(max_iterations);
count = 1;
disp(count);
initial = (b-a)*(f(b) + f(a))/2;
while count < max_iterations
disp(initial);
trap_0 = initial;
trap_1 = 0;
trap_1_midpoints = a:(0.5^count):b;
for i = 1:(length(trap_1_midpoints)-1)
trap_1 = trap_1 + (trap_1_midpoints(i+1) - trap_1_midpoints(i))*(f(trap_1_midpoints(i+1) + f(trap_1_midpoints(i))))/2;
endfor
if abs(trap_0 - trap_1) < TOL
root = trap_1;
return;
endif
intial = trap_1;
count = count + 1;
disp(count);
endwhile
disp(['Process ended after ' num2str(max_iterations), ' iterations.']);
我在 Matlab 中试过你的函数。
您的代码没有停止。而是 trap_1_midpoints 的大小呈指数增长。随着 trap_1 的计算时间也呈指数增长。这就是你所经历的拖延。
我还在您的代码中发现了一个可能的错误。我猜 if 子句后面的行应该是 initial = trap_1。检查缺失的 'i'.
这样，您的代码仍然需要很长时间，但是如果您增加容差（例如，将值设置为 1），您的代码 return.
您可以尝试矢量化 for 循环以加快速度。
编辑：我认为在您的 for 循环中，f(trap_1_midpoints(i+1) 之后缺少 )。
在 count=52 左右之后，算术序列 trap_1_midpoints 不再以任何有意义的方式用浮点数表示。在 count=1075 或类似的之后，步长不再表示为正浮点数 double 数。这就是说，绑定 max_iterations = 10000 是荒谬的。如下所述，count=20 之后的所有计算都没有意义。
步长 h 的理论误差为 O(T·h^2)。在 O(T/h) 个数字的总和中存在数值误差累积，即 O(mu/h) 与 mu=1ulp=2^(-52)。总的来说，对于双数 h=1e-5 或算法 count=17，数值积分的最低误差预计在 h=mu^(1/3) 左右。这可能会随着区间长度以及函数的平滑度或波浪度而变化。
人们可以预期错误除以四的行为，同时仅对于高于此边界的步长将步长减半 1e-5。这也意味着 abs(trap_0 - trap_1) 是仅在此步长范围内的 trap_0（和 abs(trap_0 - trap_1)/3 对于 trap_1）误差的可靠度量。
大约 h=1e-3 应满足误差界限 TOL=1e-6，这对应于 count=10。如果递归没有停止 count = 14（应该给出小于 1e-8 的错误），则该方法未准确实现。