function [calka, teoretyczny_blad] = romberg(f, a, b, dokladnosc)
M(1, 1) = (b - a) * (f(a) + f(b)) / 2;
n = 1;
h = (b - a);
m = 1;
do
n = n + 1;
h = h / 2;
S = sum(f(a + (2 * (1 : 1 : m) - 1) * h));
M(n, 1) = (M(n - 1, 1) + h * S * 2) / 2;
for i = 2 : 1 : n;
M(n, i) = M(n, i - 1) + (M(n, i - 1) - M(n - 1, i - 1)) / (4 ^ (i - 1) - 1);
end
teoretyczny_blad = abs(M(n, n) - M(n - 1, n - 1));
m = m * 2;
until(teoretyczny_blad < dokladnosc)
calka = M(n, n);
fun.m
function y = fun(x)
y = x .^ 2 .* exp(-x);
skrypt.m
format long;
a = 0;
b = 1;
f = @fun;
disp("dokladna wartosc calki");

1

calka = 1/(sqrt(4*x-1));
disp(calka);
disp("obliczenia dla dopuszczalnego bledu 1e-4");
[calka1, teoretyczny1] = romberg(f, a, b, 1e-4);
disp("wartosc kwadratury Romberga");
disp(calka1);
disp("szacunkowy blad");
disp(teoretyczny1);
disp("blad rzeczywisty");
rzeczywisty1 = abs(calka - calka1);
disp(rzeczywisty1);
disp("blad rzeczywisty mniejszy od dopuszczalnego?");
disp(rzeczywisty1 <= 1e-4);
disp("blad rzeczywisty nie wiekszy od teoretycznego?");
disp(rzeczywisty1 <= teoretyczny1);
disp("obliczenia dla dopuszczalnego bledu 1e-6");
[calka2, teoretyczny2] = romberg(f, a, b, 1e-6);
disp("wartosc kwadratury Romberga");
disp(calka2);
disp("szacunkowy blad");
disp(teoretyczny2);
disp("blad rzeczywisty");
rzeczywisty2 = abs(calka - calka2);
disp(rzeczywisty2);
disp("blad rzeczywisty mniejszy od dopuszczalnego?");
disp(rzeczywisty2 <= 1e-6);
disp("blad rzeczywisty nie wiekszy od teoretycznego?");
disp(rzeczywisty2 <= teoretyczny2);