quadvav

В воскресенье, с настроением поработать над jnum, занимался поиском одномерного адаптивного интегратора, способного интегрировать вектор-функции. Удивительно, но такого не нашел... Зато нашел алгоритм quadva (статья называется "Vectorized Adaptive Quadrature in Matlab") от проф. L.F. Shampine, который "имеет" quadl из матлаба по всем параметрам. Подобную стратегию векторизации адаптивных квадратур я применял и раньше, при вычислении конформных отображений, а потому алгоритм мне сразу понравился и я решил сделать на его основе свой интегратор на чистом J с поддержкой вектор-функций.

В результате получился интегратор quadvav (с ДОСовскими концами строк). Он, используя линейки Гаусса-Кронрода (кстати, Александра Семеновича Кронрода), быстро (особенно за счет векторизации) считает интегралы с конечными и бесконечными пределами, допуская слабые интегрируемые сингулярности на концах интервала. Работает так (предварительно нужно записать quadvav.ijs в подкаталог user установочного каталога интерпретатора J):

   load <'~user/quadvav.ijs'
([: % *:) quadvav 1 _
1 1.47143e_8
   (^ @ - @ *:) quadvav __ _
1.77245 1.2648e_5
   ([: %@%: 1 - *:) quadvav _1 1
3.14159 9.3627e_9
   ((i. 5)&(^~"0 1)) quadvav 0 1
          1         0.5    0.333333        0.25        0.2
4.96565e_17 3.52687e_17 3.80486e_17 2.59796e_17 5.43357e_7

По ходу спешу сообщить, что (того не предполагая) улучшил Матлаб. ;-) Дело в том, что я решил явно одно элементарное алгебраическое уравнение, возникающее при ослаблении потенциальных сингулярностей в случае конечного интревала интегрирования, которое в quadva решалось итерациями Ньютона. На свое письмо по этому поводу получил следующий ответ:

...The program is being modified for addition to Matlab itself at this time. ... I did consider analytical solution of the cubic, but Newton's method was so easy and effective that I implemented it. I had not appreciated how simple the solution you provide below is, so I'll have to reconsider analytical solution.

Thanks very much, Larry Shampine

Так что ENJOY ! ;-) В следующем релизе я добавлю quadvav к jnum. Может быть сделать на J cuhre ?

update (01.02.2007): Добавил интегрирование по ломаной в комплексной плоскости (в этом случае на контуре не должно быть сингулярностей) так же как и в версии quadva, которая будет включена в следующий Матлаб.