Code for romberg integration using recursion

format long
eps_step = 1e-5;
N = 11;
a = 0;
b = 1;
R = zeros( N + 1, N + 1 );
h = b - a;
pow=(1/3);
R(0 + 1, 0 + 1) = 0.5*(cos(a) + cos(b))*h;

for i=1:N
    h = h/2;
    R( i + 1, 1 ) = 0.5*(cos(a) + 2*sum(  cos( (a + h):h:(b - h) ) ) + cos(b))*h;
    for j = 1:i
        R(i + 1, j + 1) = (4^j*R(i + 1, j) - R(i, j))/(4^j - 1);
    end

    if abs( R(0 + 1, 0 + 1) -  R( i + 1, 1 ) ) < eps_step 
       break;
    elseif fprintf( 'The composite trapezoidal rule did not converge' );
        continue;
    end

   %R(0 + 1, 0 + 1) = R( i + 1, 1 );
   if abs( R(i + 1, i + 1) - R(i, i) ) < eps_step
        R(i + 1, i + 1)
       break;
    elseif i == N + 1
       error( 'Romberg integration did not converge' );
    end
end

I don't understand what's wrong with this code. I want to find what 'N' value is required for both romberg and trapezoidal and see which is greater. But the code does not run if the N value does not work for trapezoidal rule. Can anyone figure out the error in the code?

Answer 1

[Jan]

The posted code runs without an error in R2021b.

eps_step = 1e-5;
N = 11;
a = 0;
b = 1;
R = zeros( N + 1, N + 1 );
h = b - a;
pow=(1/3);
R(0 + 1, 0 + 1) = 0.5*(cos(a) + cos(b))*h;
for i=1:N
    h = h/2;
    R( i + 1, 1 ) = 0.5*(cos(a) + 2*sum(  cos( (a + h):h:(b - h) ) ) + cos(b))*h;
    for j = 1:i
        R(i + 1, j + 1) = (4^j*R(i + 1, j) - R(i, j))/(4^j - 1);
    end
    if abs( R(0 + 1, 0 + 1) -  R( i + 1, 1 ) ) < eps_step 
       break;
    elseif fprintf( 'The composite trapezoidal rule did not converge\n' );
        continue;
    end
   %R(0 + 1, 0 + 1) = R( i + 1, 1 );
   if abs( R(i + 1, i + 1) - R(i, i) ) < eps_step
        R(i + 1, i + 1)
       break;
    elseif i == N + 1
       error( 'Romberg integration did not converge' );
    end
end

The composite trapezoidal rule did not converge
The composite trapezoidal rule did not converge
The composite trapezoidal rule did not converge
The composite trapezoidal rule did not converge
The composite trapezoidal rule did not converge
The composite trapezoidal rule did not converge
The composite trapezoidal rule did not converge
The composite trapezoidal rule did not converge
The composite trapezoidal rule did not converge
The composite trapezoidal rule did not converge
The composite trapezoidal rule did not converge

"I don't understand what's wrong with this code" - please post why you assume, that there is a problem.

[Prajakta Borse]

@JanI want the code to print the statement "The composite trapezoidal rule did not converge" once and then compute the romberg integration. BUt it stops at composite trapezoidal. Why is it happening?

[Jan]

What does this mean: "BUt it stops at composite trapezoidal"?

In the question you mention an "error in the Romberg integration". Which error occurs where?

[Prajakta Borse]

It does not give an error, I want the code to compute value using Trapezoidal rule and romberg integration in the same code. But I'm not able to do that.

[Jan]

Split your problem into parts. Start with the trapezoidal rule and modify it until it works. Afterwards do the same with the Romberg integration. I do not see, that this is working currently. You find many examples in the net for working versions. Finally join the two methods.

The description "I'm not able to do that" is not usefeul. You do see a problem. Then explain, what you see. Which values do not match which expectations?