Р: Как получить кусочные коэффициенты интерполяционного сплайна для аналитического интегрирования?

Мотивация

Я численно оцениваю глубоко вложенный множественный интеграл. На каждом уровне вложенности я получаю вектор интегралов на уровне ниже, который умножается на вектор функций плотности, чтобы получить вектор подынтегральных выражений y на этом уровне. Значения x расположены неравномерно.

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

Вопрос

Я изучал такие функции, как spline и splinefun, а также функции из пакета splines2. Но я не могу найти ничего, что подскажет мне коэффициенты ряда кубических многочленов - по одному на сегмент между узлами.

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

Спасибо.


r interpolation polynomials spline integral
person Andrew Kirk    schedule 27.07.2018    source источник

Ответы (1)


arrow_upward
4
arrow_downward

Это отличная возможность расширить свой новый ответ здесь: Как сохранить и загрузить функции сплайн-интерполяции в R? С помощью этого кусочно параметризации легко вычислить кусочный интеграл.

введите описание изображения здесь

Вот (векторизованная) функция для его вычисления:

## a function for integration on a piece
piecewise_int <- function (hi, yi, bi, ci, di) {
  yi * hi + bi * hi ^ 2 / 2 + ci * hi ^ 3 / 3 + di * hi ^ 4 / 4
  }

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

## the small example in the linked thread
set.seed(0)
xk <- c(0, 1, 2)
yk <- round(runif(3), 2)
f <- splinefun(xk, yk, "natural")  ## natural cubic spline
construction_info <- environment(f)$z

## information for integration
int_info <- with(construction_info,
                 list(h = diff(x), y = y[-n], b = b[-n], c = c[-n], d = d[-n])
                 )

## cubic spline integration on all pieces
integral <- sum(do.call(piecewise_int, int_info))
#[1] 0.81375

Мы также можем выполнить численное интегрирование, чтобы проверить этот результат.

integrate(f, 0, 2)
#0.81375 with absolute error < 9e-15
person Zheyuan Li    schedule 27.07.2018