Таким чином, формула Сімпсона , називають також формулою парабол, має вид
(2.3.1)Покладемо F±1=F(±h), де Fфункція (2.1.2). Оскільки F(0)=0, F(k)(x)=
f(k-1)(x), 1£k£5то згідно формули Тейлора з залишковим членом в інтегральній формі маємо
Звідси одержуємо
(2.3.2)тому що інші члени взаємно знищуються.
Оскільки
, tÎ[0,h] то застосовуючи до інтеграла (2.3.2) узагальнену теорему про середнє, знаходимо (2.3.3)де hÎ[0,h], xÎ[-h,h] - деякі точки. Приймаючи до уваги, що
з (2.3.2), (2.3.3) приходимо до формули
(2.3.4)тобто до формули Сімпсона з залишковим членом.
3. Чисельні методи знаходження визначеного
Інтеграла зі змінною верхньою межею
У деяких випадках необхідно обчислити такі інтеграли
Можна, звичайно, розглядати його для кожного значення верхньої границі х як інтеграл зі сталими границями і обчислювати за однією з квадратурних формул, що невигідно у випадку великої кількості значень x. Краще вибрати деяку сітку і скласти таблицю значень інтеграла на цій сітці Fn=F(x) за квадратурної формули високої точності. Тоді
(3.1)причому останній інтеграл можна одчислювати за простими квадратурними формулами.
Окрім того, маючи таблицю F(xn), можна знаходити F(x) інтерполяцією за цією таблицею. Природно, маючи і похідну інтегралаF¢(x)=r(x)f(x). Краще скористатись інтерполяційним поліномом Ерміта.
4.Опис обчислювального алгоритму
При реалізаціі алгоритму обчислення визначеного інтеграла зі змінними границями інтегрування використовуються процедури та функцiї, для того щоб скоротити витрати машинного часу при обчислюваннi, та для компактностi программи. Программа для знаходження написана на мовiDelphi5, стан пограмми – вiдлажена.
5. Обговорювання результатів
Таблиця 1
Формула (3.1) | Формула Сімпсона | Формула трапецій | Дійсне значення інтеграла | |
a=0; b=1; | -0.7974398040Різниця 0.0000012883 | -0.7974386790Різниця 0.0000001633 | -0.7993252434Різниця 0.00188672780 | -0.7974385156 |
a=0; b=2; | 3.9190337956 Різниця 0.0000062805 | 3.9190353338 Різниця 0.0000047422 | 3.90875628130 Різниця 0.01028379486 | 3.9190400761 |
a=0; b=3; | 10.5498688094 Різниця 0.00002744251 | 10.5498688094 Різниця 0.00002744251 | 10.5247085565 Різниця 0.02518769537 | 10.5498962519 |
a=0; b=4; | 17.8842287345 Різниця 0.0000804723 | 17.8842201707 Різниця 0.00008903613 | 17.8382724576 Різниця 0.0460367491 | 17.8843092068 |
a=0; b=5; | 25.5043003647 Різниця 0.0001835185 | 25.5042688642 Різниця 0.00021501907 | 25.4318420115 Різниця 0.0726418717 | 25.5044838833 |
a=0; b=6; | 33.2576007639 Різниця0.00035637138 | 33.2575244054 Різниця 0.00043272988 | 33.1529684530 Різниця 0.1049886822 | 33.2579571352 |
Таблиця 1 була отримана при наступних вхідних даних:
Кількість вузлів при побудові таблиці значень інтегралу (1) =20
Кількість вузлів при застосуванні формули трапецій =20
Кількість вузлів при застосуванні формули Сімпсона =20
Таким чином з таблиці 1 видно, що чим більший проміжок ми беремо тим кращу точність отримаємо, навіть краще за формулу Сімпсона, але загальна похибка (відносно дійсного значення) також збільшується. Формулу (1) доцільно використовувати, якщо потрібно обчислити інтеграл на відносно великому проміжку та якщо треба обчислити відразу декілька інтегралів.
Список посилань
1. Каліткін Н.Н.‘Чисельні методи’ – М.: Наука, 1978. – 512 с.
2. Балашова С.Д. ‘Тексты лекций по курсу “ Численные методы”’. – Днепропетровск: Из – во ДГУ, 1989. – 206 с
3. Мусiяка В.Г. Основи чисельних методiв механiки. – Днiпропетровськ: Вид – во ДДУ, 1993. – 156 с.
4. Методические рекомендации по курсу “ Методы вычислений в инженерных расчётах”/ Составитель В.Г. Мусияка. – Днепропетровск: Из – во ДГУ, 1992. – 40 с.
5. Фіхтегольц Г.М. ‘Основи математичного аналізу’– М.: Наука, 1968. – 440 с.
Д О Д А Т К И
А Опис вихiдних даних та результатiв розрахунку
Вихiднi данi
Кількість вузлів при побудові таблиці значень інтегралу (1) nGrid - integer;
Кількість вузлів при застосуванні формули трапецій nTrap - integer;
Кількість вузлів при застосуванні формули Сімпсона nSim - integer;
Границі інтегрування a і b – real;
Наслiдки виконання програми друкуються у виглядi:
Вихідні дані це функції типа real.
FullIntegral(L,R);
integralSimpsona(L,R);
integralTrapeciay(L,R);
first(L,R);
B Схемаобчислювального алгоритму
Функція y(x) тип real Функція first(x1,x2) тип real
Процедура createGrid Процедура createGridOfInt
Процедура setN(nG, nT, nS )
Функція integralSimpsona(aSim,bSim) тип real
Функція integralTrapeciay(aTrap,bTrap) тип real
так ні
так ні
так ні
так ні
С Лiстiнг програми
unit funct;
interface
uses
Windows, Messages, SysUtils, Classes, Graphics, Controls, ExtCtrls, Math;
const
a = 0;
b = 6;
type
one_array = array of real;
Parametrs = record
nSimp, nTrap :integer;
end;
function y(x:real):real;
function first(x1,x2:real):real;
procedure setN(nG, nT, nS :integer);
procedure createGrid;
procedure createGridOfInt;
function integralSimpsona(aSim,bSim:real):real;
function integralTrapeciay(aTrap,bTrap:real):real;
function FullIntegral(aFull,bFull:real):real;
var
Xgrid:one_array; GridOfInt:one_array; nGrid, nSim, nTrap :integer;
implementation
//--------------------------------------------------------------
function y(x:real):real;
begin
Result:=2*arctan(x)-11*Power(3,(-x))+5;
end;
//--------------------------------------------------------------
function first(x1,x2:real):real;
begin
result:=5*x2+2*x2*ArcTan(x2)+11*Power(3,(-x2))/Ln(3)-Ln(1+x2*x2)-(5*x1+2*x1*ArcTan(x1)+11*Power(3,(-x1))/Ln(3)-Ln(1+x1*x1));
end;
//--------------------------------------------------------------
procedure setN(nG, nT, nS :integer);
begin
nGrid:=nG; nSim:=nS; nTrap:=nT;
end;
//--------------------------------------------------------------
procedure createGrid;
var i:integer; h:real;
begin
h:=(b-a)/nGrid;
SetLength(Xgrid,nGrid+1);
for i:=0 to nGrid do
Xgrid[i]:=a+h*i;
end;
//--------------------------------------------------------------
procedure createGridOfInt;
var i, n :integer;
begin
n:=High(Xgrid);
SetLength(GridOfInt,n+1);
for i:=1 to n do
GridOfInt[i]:=integralSimpsona(a,Xgrid[i]);
end;
//--------------------------------------------------------------
function integralSimpsona(aSim,bSim:real):real;
var X:one_array; i:integer; sum1, sum2, h:real;
begin
if(aSim<>bSim)then
begin
h:=(bSim-aSim)/nSim;
SetLength(X,nSim+1);
for i:=0 to nSim do
X[i]:=aSim+h*i;
sum1:=0;
sum2:=0;
for i:=1 to (nSim div 2) do
sum1:=sum1+y(X[2*i-1]);
for i:=1 to ((nSim div 2)-1) do
sum2:=sum2+y(X[2*i]);
Result:=(bSim-aSim)*(y(X[0])+y(X[nSim])+4*sum1+2*sum2)/(3*nSim);
end
else
Result:=0;
end;
//--------------------------------------------------------------
function integralTrapeciay(aTrap,bTrap:real):real;
var i:integer; sum, h :real; X:one_array;
begin
h:=(bTrap-aTrap)/nTrap;
SetLength(X,nTrap+1);
for i:=0 to nTrap do
X[i]:=aTrap+h*i;
sum:=(y(X[0])+y(X[nTrap]))/2;
for i:=1 to (nTrap-1) do
sum:=sum+y(X[i]);
Result:=sum*h;
end;
//--------------------------------------------------------------
function FullIntegral(aFull,bFull:real):real;
var z1, z2, z3, raznost :real; i, ai, bi :integer;
begin
if(aFull<>bFull)then
begin
for i:=0 to (High(Xgrid)-1) do
begin
if((Xgrid[i]<=aFull)and(aFull<=Xgrid[i+1]))then ai:=i;
if((Xgrid[i]<=bFull)and(bFull<=Xgrid[i+1]))then bi:=i;
end;
raznost:=GridOfInt[ai]+integralTrapeciay(Xgrid[ai],aFull);
Result:=GridOfInt[bi]+integralTrapeciay(Xgrid[bi],bFull)-raznost;
end
else
Result:=0;
end;
//--------------------------------------------------------------