Code for romberg integration using recursion
13 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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?
0 commentaires
Réponses (1)
Jan
le 26 Oct 2021
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
"I don't understand what's wrong with this code" - please post why you assume, that there is a problem.
4 commentaires
Jan
le 28 Oct 2021
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?
Voir également
Catégories
En savoir plus sur Symbolic Math Toolbox dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!