Composite trapezoid rule not running in Octave

Question

I have the following code in Octave for implementing the composite trapezoid rule and for some reason the function only stalls whenever I execute it in Octave on f = @(x) x^2, a = 0, b = 4, TOL = 10^-6. Whenever I call trapezoid(f, a, b, TOL), nothing happens and I have to exit the Terminal in order to do anything else in Octave. Here is the code:

% 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.']);

Answer 1

I have tried your function in Matlab. Your code is not stalling. It is rather that the size of trap_1_midpoints increases exponentionaly. With that the computation time of trap_1 increases also exponentionaly. This is what you experience as stalling.

I also found a possible bug in your code. I guess the line after the if clause should be initial = trap_1. Check the missing 'i'.

With that, your code still takes forever, but if you increase the tolerance (e.g. to a value of 1) your code return.

You could try to vectorize the for loop for speed up.

Edit: I think inside your for loop, a ) is missing after f(trap_1_midpoints(i+1).