численное интегрирование

 

задачи численного интегрирования приходится решать для

для инженера это означает - всегда

 

вместо функции, которую требуется проинтегрировать, интегрируется интерполяционный полином

 

пусть требуется определить значение интеграла функции f(x) на отрезке [a,b]

этот отрезок делится точками на n равных отрезков, каждый длиной

d = (b - a) / n

обозначим через y1, y2, .. yn значение функции f(x) в точках x1, x2, .. xn, каждая из которых принадлежит одному из элементарных отрезков

далее составляем сумму:

S = y1 * d + y2 * d + .. + yn * d

эта сумма приближённо выражает интеграл для f(x) на [a,b]

 

S = d * (y1 + y2 + .. + yn)

 


 

метод прямоугольников - полином нулевого порядка

если функцию на каждом из частичных отрезков аппроксимировать прямой, параллельной оси абсцисс и проходящей через начало (или через середину отрезка, или через конец отрезка), то получим метод прямоугольников

 

площадь прямоугольника на каждом отрезке:

S = yi * d

 

полная формула в случае деления всего промежутка интегрирования на n отрезков одинаковой длины d:

S = d * (y1 + y2 + .. + yn)

очевидно, стоит рассчитывать на бо́льшую точность если брать в качестве опорной точки для нахождения высоты точку посередине промежутка

оценка погрешности


          ε   ∝   y'(ζ) * n * d2 / 4  

          ε   ∝   Ο (1/n)

 


 

метод трапеций - полином первого порядка

если функцию на каждом из частичных отрезков аппроксимировать прямой, проходящей через конечные значения, то получим метод трапеций

 

площадь трапеции на каждом отрезке:

S = d * fi = d * yi + d * (yi+1 - yi) / 2

полная формула в случае деления всего промежутка интегрирования на отрезки одинаковой длины d:

S = d * (f0 / 2 + f1 + f2 + .. + fn-1 + fn / 2)

оценка погрешности

 
          ε   ∝   y"(ζ) * n * d3 / 12  

          ε   ∝   Ο (1/n2)

 


 

метод парабол - полином второго порядка

использовав три точки отрезка интегрирования, можно заменить подынтегральную функцию квадратичной параболой, проходящей через точки. обычно в качестве таких точек используют концы отрезка и его среднюю точку

выведем формулу Симпсона:

рассматриваем функцию
  
f(x) = a * x2 + b * x + c

рассмотрим отрезок [0, 2*h]. очевидно, что 

f(0)   = c                              (1)
f(h)   = a * h2 + b * h + c             (2)    
f(2*h) = 4 * a * h2 + 2 * b * h + c     (3)

        2h
найдем  ∫  f(x) dx  
        0
 
  2h                                        2h
  ∫ f(x) dx   =  a * x3/3 + b * x2/2 + c * x |   =   (h / 3) * (8 * a * h2 + 6 * b * h + 6 * c)
  0                                          0

используем равенства (1-3), чтобы осуществить переход:
 
    2h
    ∫ f(x) dx   =  (h / 3) * (f(0) + 4 * f(h) + f(2*h))  
    0

 

итак, формула для площади на отрезке d со средней точкой h, имеет вид:

S = d * fi = d * ( yi-1 + 4 * yi + yi+1 ) / 3

 

тогда формула для [a, b] с равномерной сеткой и четным числом шагов n (при длине шага d) будет:

S = d * (ya + 4*f1 + 2*f2 + 4*f3 + 2*f4 + ... + 2*fn-2 + 4*fn-1 + yb) / 6

оценка погрешности

 ε   ∝   Ο (1/n4)  

 


 

калькулятор

независимо от выбранного метода в процессе численного интегрирования необходимо оценить погрешность. погрешность уменьшается при увеличении количества разбиений. однако при этом возрастает погрешность округления

 

функция f(x) "+ - * / ^"
нижний предел
верхний предел
шаг


результат


матричное интегрирование

function s = matinteg (a, b, n, f)
  A = tril (ones(n)) ;
                      
  t = linspace (a, b, n) ;
  y = arrayfun (f, t) ; 
  S = (((b - a)/n) * A) * y' ; 

  S(n) 
end

octave > matinteg (2,4,100, @ (x) x + 2)
 ans = 10