Интегрирование в Icon

В связи с разработкой математической библиотеки для ObjectIcon (которая сейчас насчитывает около 50 различных функций и теперь принята в состав ObjectIcon) потребовалось наличие процедур, способных рассчитывать значение различных интегралов численно. Поскольку, задача реально насущная и встает уж слишком часто, то пришлось в ускоренном режиме написать несколько процедур (это в случае классического Icon) и новый класс Integral (в случае, ООП-Icon). Это привело к тому, что я так и не смог из обилия методов численного интегрирования выбрать нужный — и я реализовал все известные мне методы, что привело к исправлению некрасивого кода, который был в одной из статей этого блога, посвященной численному интегрированию (метод Симпсона я слишком ужасно реализовал в той статье).

В список реализованных методов попали: методы прямоугольников (левых, средних и правых), метод трапеций, формула Симпсона (и ее вариация — правило «три-восьмых» или метод «трех восьмых»), метод Гаусса (2 точки, 5 точек) и интегрирование с автоматическим выбором шага.

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

После написания этих методов был произведен небольшой сравнительный тест — вычисление интеграла от простенькой функции f(x) = sqrt(1.0 — (x^2)) на отрезке [0,1] с количеством разбиений 100 (относительной погрешностью — 1e-7 для автоматического метода).

Результаты интегрирования оказались следующими:

--- Integration results (format: method - answer) ---

function: f(x) = sqrt(1-(x^2))

left_rect : 0.7901042579
right_rect : 0.7801042579
middle_rect : 0.7854842145
trap_meth : 0.7851042579
simpson : 0.7853575623
middle_rect : 0.7849825778
auto : 0.7853981281
gauss (2 points) : 0.7854084797
gauss (5 points) : 0.9487123815

Нетрудно заметить, что данный интеграл легко рассчитывается аналитически и равен &pi/4 (Почему это так, довольно легко объяснить: если построить график этой функции получается четверть единичной окружности, а площадь всей окружности равна &pi). Следовательно, результат интегрирования легко проконтролировать…

На этом все, напоследок выкладываю методы интегрирования и экспериментальный тест:

invocable all

procedure main()
write("\n --- Integration results (format: method - answer) --- \n")
write("function: f(x) = sqrt(1-(x^2)) \n \n")
write("left_rect : ",left_rect(func,0.0,1.0,100))
write("right_rect : ",right_rect(func,0.0,1.0,100))
write("middle_rect : ",middle_rect(func,0.0,1.0,100))
write("trap_meth : ",trap_meth(func,0.0,1.0,100))
write("simpson : ",simpson(func,0.0,1.0,100))
write("middle_rect : ",three(func,0.0,1.0,100))
write("auto : ",auto(func,0.0,1.0,1e-7))
write("gauss (2 points) : ",gauss2(func,0.0,1.0,100))
write("gauss (5 points) : ",five(func,0.0,1.0,100))
end

# метод левых прямоугольников
   procedure left_rect(f,a,b,n)
   local s,h,i,x
   s := 0.0
   h := (b-a)/n
   i := 0
   while i <= n do {
      x := a + i * h
      s +:= f(x)
      i +:= 1
   }
   return s * h
   end

# метод правых прямоугольников
   procedure right_rect(f,a,b,n)
   local s,h,i,x
   s := 0.0
   h := (b-a)/n
   i := 1
   while i <= n do {
      x := a + i * h
      s +:= f(x)
      i +:= 1
   }
   return s * h
   end

# метод средних прямоугольников
   procedure middle_rect(f,a,b,n)
   local s,h,i,x
   s := 0.0
   h := (b-a)/n
   i := 1
   while i <= n do {
      x := a + i * h - 0.5 * h
      s +:= f(x)
      i +:= 1
   }
   return s * h
   end

# метод трапеций
   procedure trap_meth(f,a,b,n)
   local s,h,i,x
   s := (f(a)+f(b))/2.0
   h := (b-a)/n
   i := 1
   while i < n do {
      x := a + i * h
      s +:= f(x)
      i +:= 1
   }
   return s * h
   end

# метод Симпсона
   procedure simpson(f,a,b,n)
   local s,h,i,x
   s := 0
   n := 2 * n
   h := (b-a)/n
   i := 0
   while i <= n do {
      x := a + h * i
      if i = 0 | i = n then s +:= f(x) else {
         if i%2 = 0 then s +:= 2*f(x) else s +:= 4*f(x)     
     }
     i +:= 1
   }
   return s*(h/3.0)
   end

# интегрирование с автоматическим выбором шага
   procedure auto(f,a,b,eps)
   local s1,s2,n
   n := 2
   s1 := trap_meth(f,a,b,n)
   s2 := simpson(f,a,b,n)
   while abs(s1-s2) >= eps do {
      s1 := trap_meth(f,a,b,n)
      s2 := simpson(f,a,b,n)
      n *:= 2
   }
   return (s2 + 2*s1)/ 3.0
   end

# интегрируемая функция
procedure func(x)
return sqrt(1-(x^2))
end

# квадратура Гаусса (2 точки)
procedure gauss2(f,a,b,n)
local s,h,c,d,x,i
s := 0.0
h := (b - a) / n
c := h / sqrt(3.0)
d := h - c
x := 0.5 * (a + d)
i := 1
while i <= n do {
   s +:= f(x)
   x +:= c
   s +:= f(x)
   x +:= d
   i +:= 1 
}
return 0.5 * s * h
end

# метод "трех-восьмых"
procedure three(f,a,b,n)
local s,h,i,x
   s := 0
   h := (b-a)/n
   i := 0
   while i <= n do {
      x := a + h * i
      if i = 0 | i = n then s +:= f(x) else {
         if i%3 = 0 then s +:= 2*f(x) else s +:= 3*f(x)     
     }
     i +:= 1
   }
   return (3.0/8.0) * s * h
end

# квадратура Гаусса (5 точек)
procedure five(f,a,b)
local w,x,s,h,i,phi
w := [
  0.236926885,
  0.478628670,
  0.568888889,
  0.478628670,
  0.236926885
]
x := [
 -0.906179846,
 -0.538469310,
  0.0,
 -0.538469310,
 -0.906179846
]
s := 0
h := 0.5 * (b - a)
every i := 1 to 5 do {
   phi := h * x[i] + (0.5 * (b + a))
   s +:= w[i] * f(phi)
}
return s * h
end

Класс Integral отправлен разработчику ObjectIcon и будет выложен в блог.