Интегрирование функции на основе численных методов Matlab

[Olivia]

Добрый день! Необходимо вычислить интеграл методами средних прямоугольников, трапеции и Симпсона. При этом встроенными функциями пользоваться нельзя, кол-во разбиений должно накапливаться, а не изначально задаваться. Программу написала и она работает, но сказали, что формулы для I_old и I_new неверные. Кто может помочь, как исправить, потому что работаю с Matlab всего неделю Очень нужна ваша помощь

%Метод средних прямоугольников
f=@(x)(3*x.^5+4*x.^3+1);
a=0;
b=2;
eps = 0.001; % точность

Sum=0; % метод средних прямоугольников
n=1;
h=(b-a);
I_old=(f(a)+f(b))*h/2;
n=2;
h=h/2;
I_new=I_old/2+h*f(a+h)*h;
eps_pr=abs(I_new-I_old)/3;
while(eps_pr>eps)
I_old=I_new;
Sum=0;
h=h/2;
n=2*n;
for j=0:n-1
Sum=Sum+f(a+j*h+h/2);
end;
I_new=Sum*h;
eps_pr=abs(I_new-I_old)/3;
eps_teor=(528*h.^2*(b-a))/24;
end;
I_trap = (4*I_new-I_old)/3;
disp('Алгоритм средних прямоугольников');
disp(['Значение интеграла: ',num2str(I_new)])
disp(['Шаг: ',num2str(h)])
disp(['Погрешность: ',num2str(eps_pr)])
disp(['Кол-во разбиений: ', num2str(n)])
disp(['Теоретическая погрешность:', num2str(eps_teor)])
disp(['Эсктраполированное значение: ', num2str(I_trap)])




%Метод трапеции
f=@(x)(3*x.^5+4*x.^3+1);
a=0;
b=2;
eps = 0.001;

Sum=0; % метод трапеции
n=1;
h=(b-a);
I_old=(f(a)+f(b))*h/2;
n=2;
h=h/2;
I_new=I_old/2+h*f(a+h)*h; eps_pr=abs(I_new-I_old)/3;
while(eps_pr>eps)
I_old=I_new;
Sum=0;
h=h/2;
n=2*n;
for j=0:n-1
Sum=Sum + (f(a+j*h) + f(a+(j+1)*h));
end;
I_new=Sum*h/2;
eps_pr=abs(I_new-I_old)/3;
eps_teor=(528*h.^2*(b-a))/12;
end;
I_trap = (4*I_new-I_old)/3;
disp('Метод трапеций');
disp(['Значение интеграла: ',num2str(I_new)])
disp(['Шаг: ',num2str(h)])
disp(['Погрешность: ',num2str(eps_pr)])
disp(['Кол-во разбиений: ', num2str(n)])
disp(['Теоретическая погрешность:', num2str(eps_teor)])
disp(['Эсктраполированное значение: ', num2str(I_trap)])



%Метод Симпсона

f=@(x)(3*x.^5+4*x.^3+1);
a=0;
b=2;
eps = 0.001;

Sum_n=0; % метод Симпсона
Sum_m=0;
n=1;
h=(b-a);
I_old=(f(a)+f(b))*h/2;
n=2;
h=h/2;
I_new=I_old/2+h*f(a+h)*h;
eps_pr=abs(I_new-I_old)/15;
while(eps_pr>eps)
I_old=I_new;
Sum=0;
Sum_s=0;
Sum_n=0;
Sum_m=0;
h=h/2;
n=2*n;
for j=2:2:n-1
Sum_n=Sum_n+f(a+j*h);
end;

for j=1:2:n
Sum_m=Sum_m+f(a+j*h);
end;

Sum_s=f(a)+f(b);
Sum=Sum_s+2*Sum_n+4*Sum_m;
I_new=Sum*(h/3);
eps_pr=abs(I_new-I_old)/15;
end;
Itrap = (16*I_new-I_old)/15;
disp('Алгоритм Симпсона');
disp(['Значение интеграла: ',num2str(I_new)])
disp(['Шаг: ',num2str(h)])
disp(['Погрешность: ',num2str(eps_pr)])
disp(['Кол-во разбиений: ', num2str(n)])
disp(['Теоретическая погрешность:', num2str(eps_teor)])
disp(['Эсктраполированное значение: ', num2str(I_trap)])