Лабораторная работа №3 приближенное вычисление определенных интегралов. Методы прямоугольников и метод симпсона

Пример применения численного интегрирования.
Вычислим, например, интеграл по формуле Симпсона с точностью до 0,001.
Чтобы выбрать необходимое для получения заданной точности число 2n, найдем f⁽⁴⁾(x). Последовательно дифференцируя функцию f(x)= , получаем
f⁽⁴⁾(x)=4 (4х⁴-12х²+3)
Так как на отрезке [0, 1] £1, ½4х⁴-12х²+3½£5, то . Следовательно, можно взять М=20. Используя формулу оценки погрешности, имеем 20/2880n⁴<1/1000( ), откуда n⁴ >1000/144. Для того чтобы выполнялось это неравенство, достаточно взять n=2, т.е. 2n=4.
Разобьем теперь отрезок [0, 1] на четыре равные части точками х₀=0, х₁=1/4, х₂=1/2, х₃=3/4, х₄=1 и вычислим приближенно значения функции f(x)= в этих точках у₀=1,0000, у₁=0,9394, у₂=0,7788, у₃=0,5698, у₄=0,3679. Применяя формулу Симпсона, получаем
Таким образом, с точностью до 0,001. Итак, разбив отрезок [0, 1] всего на четыре равные части и заменив рассматриваемый интеграл суммой, стоящей в правой части формулы Симпсона, мы вычислили данный интеграл с необходимой точностью.


Программная реализация в MatLab.

1. 
   Вычисление определенного интеграла с помощью встроенной функции:
   

>> syms x

>> I = double(int((x^2-6.5)/(x^3+sin(7.8*x)),3.6,5.4))
>>
I =
0.26649340577279

2. 
   Программный блок, реализующий приближенное вычисление определенного интеграла:
   

//str – выбор метода приближенного вычисления определенного интеграла; a, b – пределы интегрирования; n – количество точек
function res = Integral(str, a , b, n);
h = (b-a)/n;
Sum = 0;
switch (str)


// метод Симпсона
case 'Simpson'
q1=0; q2=0;
for i = 1:n-1
if (mod(i,2)==1)
q1 = q1+f(a + i*h);
else q2 = q2+f(a + i*h);
end;
end;
Sum =(h/3)*( f(a) + 4*q1 + 2*q2 + f(b));


//метод трапеций
case 'Trapezium'
for i = 1:n-1
Sum = Sum + f(a + i*h);
end;
Sum = (h/2)*( f(a) + 2*Sum + f(b));


//метод левых прямоугольников
case 'LeftRectangles'
for i=1:n-1
Sum = Sum + h*f(a + i*h);
end;
//метод правых прямоугольников
case 'RightRectangles'
for i=1:n-1
Sum = Sum + h*f(a + i*h);
end;
//метод средних прямоугольников
case 'MediumRectangles'
for i=1:n
Sum = Sum + h*f(((a + (i-1)*h)+(a + i*h))/2);
end;

otherwise

error('This is impossible value')
end
res = Sum;
return

3. 
   Программный блок вычисления апостериорных погрешностей вычислений строится аналогично с использованием оператора выбора:
   

//вычисление погрешности для формулы трапеций
…
case 'Trapezium'
R = abs(Integral(str,a,b,n)-Integral(str,a,b,0.5*n))/3;
…

4. 
   Результаты вычислений:
   

>> format long
>> IS = Integral('Simpson',3.6,5.4,12)
IS =
0.26649801994399
>> IT = Integral('Trapezium',3.6,5.4,12)
IT =
0.26635817701576
>> ILR = Integral('LeftRectangles',3.6,5.4,12)
ILR =
0.26584089261870
>> IRR = Integral('RightRectangles',3.6,5.4,12)
IRR =
0.26687546141281
>> IMR = Integral('MediumRectangles',3.6,5.4,12)
IMR =
0.26656137337914
…
>> Priori_ErrorIT = Priori_Error('Trapezium',3.6,5.4,12)
Priori_ErrorIT =
1.040388407134533e-005
…