Численное интегрирование методом Симпсона

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

Один из широко известных методов — метод Симпсона, который достаточно легко реализуется. Однако, в этом деле есть кое-какая изюминка.

Сначала сам код процедуры:

procedure simpson(a,b,n,func)
local j,j1,j2,j3,i,f
j1:=(b-a)/(6.0*n)
j2:=0.0
j3:=0.0
i:=1.0
while i<((2*n)-1) do {
 f:=func(a+(((b-a)*i)/(2*n)))
 j2+:=f
 i+:=2.0
}
i:=2.0
while i<((2*n)-2) do {
  f:=func(a+(((b-a)*i)/(2*n)))
  j3+:=f
  i+:=2.0	
}
j:=j1*(func(a)+4*j2+2*j3+func(b))
return j
end

Как видно, процедура является голой реализацией известной формулы Симпсона для численного интегрирования, и ее реализация затруднений не вызывает. Но сама процедура принимает четыре аргумента: a — начальная точка интегрирования, b — конечная точка интегрирования, n — число разбиений отрезка [a;b], и далее самое интересное — аргумент func, который требует функцию, которую будем интегрировать.

И вот про последний аргумент надо сказать пару слов.

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

invocable "имя_функции"

а затем можно определять саму подынтегральную функцию.

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

invocable all

которая разрешает использование всех встреченных в файле функций f(x) в виде «f»(x), а не только f(x).

Ну а теперь на сладкое, посчитаем интеграл для функции x^2 от 1 до 5. Определяем функцию sqr(x) (от английского «square» — квадрат):

procedure sqr(x)
return x*x
end

Определяем ее как invocable (в самом начале программного файла):

invocable "sqr"

Cчитаем:

simpson(1,5,200,"sqr")

Получаем результат: 40.835996 (истинное значение — 41.333333). Что же, довольно близко при небольшом количестве итераций.

Код всего расчета:

invocable "sqr"
procedure main() 
write(simpson(1,5,200,"sqr"))
end

procedure simpson(a,b,n,func)
local j,j1,j2,j3,i,f
j1:=(b-a)/(6.0*n)
j2:=0.0
j3:=0.0
i:=1.0
while i<((2*n)-1) do {
 f:=func(a+(((b-a)*i)/(2*n)))
 j2+:=f
 i+:=2.0
}
i:=2.0
while i<((2*n)-2) do {
  f:=func(a+(((b-a)*i)/(2*n)))
  j3+:=f
  i+:=2.0	
}
j:=j1*(func(a)+4*j2+2*j3+func(b))
return j
end

procedure sqr(x)
return x*x
end

Добавить комментарий