Интерполяция методом конечных разностей

[antoha1]

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

[_Ivana]

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

[antoha1]

Рассмотрим квадратичный многочлен p(x) = 2x² − 3x + 2 и предположим, что мы знаем его табличные значения p(0), p(0,1), p(0,2), p(0,3), p(0,4) и т. д. В представленной ниже таблице первая колонка содержит табличные значения полинома, вторая — разность между двумя верхними соседними значениями из первой колонки, а третья — разность между двумя соседями из второй колонки:
p(0)=2,0
2,0−1,72=0,28
p(0,1)=1,72 0,28−0,24=0,04
1,72−1,48=0,24
p(0,2)=1,48 0,24−0,20=0,04
1,48−1,28=0,20
p(0,3)=1,28 0,20−0.16=0.04
1,28−1,12=0,16
p(0,4)=1,12

Заметим, что значения в третьем столбце — одинаковы. Это не случайность. Фактически, если мы начинаем таблицу для любого полинома степени n, колонка с номером n + 1 будет всегда содержать константу. Этот решающий факт делает работоспособным указанный метод.

Мы составляли таблицу слева-направо, но точно также её можно рассчитать и справа-налево, вычислив, таким образом, недостающие значения полинома.

Для вычисления p(0,5) мы воспользуемся значениями самой нижней диагонали. Начиная с самого нижнего значения в последней колонке 0.04. Затем продолжим вторую колонку вычтя 0,04 из 0,16 и получив значение 0,12. Таким же образом мы заполним первую колонку, вычитая из её нижнего значения 1,12 полученное нами на предыдущем шаге число 0,12 из второй колонки. p(0,5) будет равным 1,12-0,12 = 1,0. Для вычисления p(0,6) используется тот же самый алгоритм: берётся 0,04 из третьей колонки, вычитается из нижнего значения (теперь уже 0,12) во второй колонке, получившееся 0,08 прописывается в нижнюю часть второй колонки и затем вычитается из нижнего значения в первой колонке (как мы помним — 1,0). Результат — 0,92 является значением p(0,6).

Процесс может быть продолжен до бесконечности. Значения полинома получаются при этом без применения операции умножения. На этом факте, в частности, была основана работа разностной машины Чарльза Бэббиджа. Для выполнения каждого следующего цикла расчёта значений квадратичного полинома, достаточно сохранить 2 числа (последние элементы, во второй и в первой колонке); для табулирования полиномов степени n число требуемых значений больше, — а именно, требуется сохранить n значений.


Вроде тут понятно написано, но не пойму как реализовать это программно. Т.е. у пользователя будет запрашиваться значения первых трех значений полинома, а программа будет выводить все остальные? Или я все неправильно понял, извиняюсь сразу за тупизм...

[Toucan]

Тема перемещена из форума «Математика (общие вопросы)» в форум «Программирование»

[ewert]

antoha1 в сообщении #637312 писал(а):

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

Судя по приведённому тексту -- именно так.

Конечные разности определяются рекуррентно:

$\Delta^0p_i=p_i\ \ $ (разности нулевого порядка -- это сами значения функции);
$\Delta^1p_i=p_{i+1}-p_i;$
$\Delta^2p_i=\Delta^1p_{i+1}-\Delta^1p_i;$
$\Delta^3p_i=\Delta^2p_{i+1}-\Delta^2p_i$

и т.д. Если степень многочлена равна $n$, то для его определения необходимы $(n+1)$ начальных значений: $p_0,p_1,\ldots,p_n$. По ним находятся разности первого порядка: $\Delta^1p_0,\Delta^1p_1,\ldots,\Delta^1p_{n-1}$, потом второго: $\Delta^2p_0,\Delta^2p_1,\ldots,\Delta^2p_{n-2}$ и так до тех пор, пока не доберёмся до единственной разности старшего порядка $\Delta^np_0$ (она, в отличие от предыдущих, при увеличении нижнего индекса меняться не будет). Из всей этой таблички надо запомнить только последние разности каждого порядка: $\Delta^0p_n=p_n,\Delta^1p_{n-1},\Delta^2p_{n-2},\ldots,\Delta^np_{0}$. Для удобства переобозначим их как $p_n,p^{(1)}_n,p^{(2)}_n,\ldots,p^{(n)}_n$. Тогда, согласно определению конечных разностей, для каждой следующей точки этот массив пересчитывается также рекуррентно, но уже сверху вниз:

$p^{(n-1)}_{i+1}=p^{(n-1)}_{i}+p^{(n)}_{i+1},$
$p^{(n-2)}_{i+1}=p^{(n-2)}_{i}+p^{(n-1)}_{i+1},$
$\ldots,$
$p^{(1)}_{i+1}=p^{(1)}_{i}+p^{(2)}_{i+1},$
$p_{i+1}=p_{i}+p^{(1)}_{i+1}$

(напомню, что $p^{(n)}_{i+1}$ пересчитывать не надо -- она остаётся постоянной). Для программирования этой процедуры нужны несколько соотв. циклов, причём двумерного массива для таблицы не нужно -- достаточно двух одномерных: одного для разностей текущего порядка на этапе вычисления таблицы и одного для хранения последних разностей всех порядков для дальнейшей работы. Хватило бы даже и одного массива, т.к. по мере расчёта таблицы строка разностей текущего порядка укорачивается, а строка последних разностей -- удлинняется; но с двумя массивами проще и, соответственно, надёжнее.