Лекция числово интегриране

Най-общо алгоритъмът за числово интегриране на функция може да се представи като сумиране на правоъгълни изрезки с различна ширина ω_i:

По принцип сумата в (2) ще дава точния интеграл когато N->∞, но за някои класове от функции резултата е точен и за крайно N. Различните алгоритми за интегриране се различават единствено по това как се избират точките и теглата в сумата (2).

Фиг.1 Интегралът е площта под кривата f(x) в интеграционните граници от а до b.

В таблица 1 са дадени сравнителни характеристики на трите алгоритъма, които ще разгледаме.

Таблица 1

Метод	Предимства	Недостатъци
Правило на трапеца	Простота. Оптимален за неправилни интеграли.	Нужни са много подинтервали за постигане на добра точност
Правило на Симпсон 1/3	Простота. По-висока точност от трапеца.	Работи само с четен брой подинтервали.
Квадратура на Гаус	Висока точност. Не се използват функционалните стойности в краищата на интеграционния интервал.	Извадките не са разположени еквидистантно.

1. Правило на трапеца

Нека имаме N равноотстоящи точки в интеграционния интервал [a,b], включително крайните точки т.е. имаме N-1 интервала с дължина:

(3)

Фиг.2 Правило на трапеца

Както се вижда от фиг.2 функцията f(x) се апроксимира с права линия във всеки подинтервал i, така че i-тата изрезка от площта се апроксимира с трапец, откъдето идва названието на метода. Стойноста на функцията f в подинтервала i се взима като средната височина (f_i +f_i+1)/2. Така площта на един трапец се дава от:

(4)

За да приложим правилото на трапеца за целия интеграционен интервал трябва да сумираме по всички подинтервали и получаваме:

Вижда се че понеже всяка вътрешна точка се отчита два пъти (има два съседа) нейната тежест е h. Крайните точки се отчитат веднъж и имат тежест h/2.

Според общото представяне (2) теглата в сумата се дават от вектора:

(6)

Грешката при правилото на трапеца е от порядък h².

Правилото на Симпсон 1/3 означава че апроксимацията се определя от три точки в два съседни подинтервала. Съществуват по-точни методи на Симпсон, например 3/8, където апроксимацията се определя от 4 точки в три съседни подинтервала. Ние ще се ограничим със стандартното правило 1/3. При това правило (вж. фиг.3) функцията във всеки подинтервал се апроксимира, не с линейна зависимост, както е при трапеца, а с полином от втора степен (парабола):

За да свържем параметрите α, β и γ с функцията ще разгледаме частния случай когато интеграционния интервал е от -1 до +1. Тогава имаме:

Фиг.3 Правило на Симпсон

Можем лесно да изразим параметрите α, β и γ чрез стойностите на функцията в двете крайни и в средната точки т.е.

(9)

Следователно имаме:

От (8) и (10) можем да представим интеграла като претеглена сума от стойностите на подинтегралната функция в три точки:

Понеже са необходими три точки можем да обобщим резултата (11) за два съседни подинтервала, които се определят от три точки. Така можем да запишем:

Правилото на Симпсон изисква интегрирането да бъде по двойки съседни интервали, затова общият брой на подинтервалите трябва да е четно число, съответно броя на точките да бъде нечетен. За да приложем правилото на Симпсон върху целия интеграционен интервал, ние сумираме приноса от всички съседни двойки, при което всяка нечетна точка освен първата и последната се отчитат по два пъти. Така получаваме:

Теглата в сумата на общото представяне се дават от вектора:

Сумата от теглата е добра проверка за правилноста на интегриране:

Грешката при правилото на Симпсон 1/3 е от порядък h⁴. Можем да споменем, че съществува метод на Симпсон с апроксимация с полином от трета степен, наречен правило на Симпсон 3/8. При него грешката е от същия порядък h⁴, а изчисленията се усложняват , и интервалите трябва да се сумират по тройки. Затова този метод се използва в редки случаи.

При този метод точките x_i, в оценъчната сума са N-те корени на полином на Лежандър от N-ти ред P_N(x_i)=0. Това гарантира, че между избраните точки няма особени точки. Теглата в оценъчната сума се определят от връзките:

(16)

Така избрани теглата гарантират, че апроксимационната грешка практически изчезва ако подинтегралната функция f(x) е полином от 2N-1, или по-нисък, ред.

Общо изискване е функцията f(x) да бъде гладка в интеграционния интервал и тогава точноста на метода превъзхожда предходните два метода на еквидистантни точки. Обратно, ако функцията не е гладка, примерно съдържа шум, тогава трябва да се предпочетат алгоритми от по-нисък порядък.

Точките и теглата за Гаусoва квадратура (y_i,w_i’) се изчисляват за интервала [-1,1] и за да могат да се приложат в произволен интервал [A,B] трябва да се извърши подходящо скалиране със смяна на променливите както следва:

• 
  Използвайте готова процедура за определяне корените и връзките в метода на Гаус, зададена във файла gauss.h
  
• 
  Изчислете относителната грешка за трите алгоритъма:
  

Изведете резултатите във файл, подредени в следната таблица:

N T-rule ε_T S-rule ε_S G-quad ε_G
Пробвайте различен брой на извадките N (трябва да са нечетни заради правилото на Симпсон), вариращи от 3 до 501

• 
  Използвайте MATLAB за да построите логаритмичните зависимости на грешките от броя на точките т.е. функция на log₁₀ε от log₁₀N. Забележете, че грешката на алгоритъма намалява линейно в логаритмичен мащаб докато се достигне грешката на закръгляване на компютъра.
  

Защо могат да се появят проблеми при изчисляване на интеграла:

Фиг. 4 Логаритмични графики на грешките при различните методи на интегриране на експоненциална функция в зависимост от броя на точките (извадки) при двойна точност (ляво) и единична точност (дясно).
Приложение 1. Java примерна програма за интегриране на функцията f(x)=x² с метода на трапеците:

Приложение 2. Java примерна програма за интегриране на функцията f(x)=sin(x) с метода на Симпсон: