MATLAB.Exponenta
–Û·Ë͇ Matlab&Toolboxes

Spline Toolbox

Интерполяционный полином в форме Ньютона

Форма Лагранжа записи интерполяционного полинома, рассмотренная в разделе Интерполяционный полином в форме Лагранжа, не является единственно возможной. Кроме этого, она имеет определенные неудобства. Так предположим, задан набора точек и по нему построен интерполяционный полином в форме Лагранжа

Тогда очевидно, что при добавлении новой точки к исходному набору точек придется заново конструировать интерполяционный полином. Этого недостатка лишен интерполяционный полином, записанный в рассмотренной ниже форме Ньютона.

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

Для вывода интерполяционного полинома в форме Ньютона понадобится понятие о разделенных разностях.

Разделенные разности и некоторые их свойства.

Пусть нам заданы различные точки и функция f(x). Первой разделенной разностью функции в точках называется величина

Разделенные разности более высоких порядков определяются рекуррентно. Например разделенная разность второго порядка в точках определяется через две разделенные разности первого порядка следующим образом:

Разделенная разность второго порядка в точках определяется так:

и т.д.

В общем случае разделенная разность s-го порядка определяется через две разделенные разности -го порядков

Разделенные разности удобно представить в виде треугольной таблицы

Для вывода в командное окно и в текстовый файл таблицы разделенных разностей некоторой функции можно воспользоваться приведенной ниже функцией DDIF. В ее входных аргументах задаются:

  1. массив точек;
  2. имя функции, для которой следует построить таблицу разделенных разностей;
  3. имя файла, в который нужно записать эту таблицу.

Поскольку в первом столбце таблицы записываются значения функции через строку, то высота таблицы равна 2*N-1, где N — число точек, т.ев вышеприведенных формулах. Ширина таблицы равна N, т.е. числу точек.

В функции DDIF использованы следующие функции MATLAB:

length — возвращает длину вектора.
feval — вычисляет функцию (по ее имени или указателю на нее) от заданного аргумента, т.е. если есть функция myfun, то вместо y = myfun(x)можно написать y = feval(‘myfun’, x) или y = feval(@myfun, x).
reshape — перестраивает массив и приводит его к заданным размерам; используется для получения вектор-столбца значений функции, поскольку заданная функция может возвращать вектор-строку.
NaN — создает массив заданных размеров, заполненный NaN (Not a Number, т.е. не число); используется для того, чтобы потом в текстовом массиве заменить NaN на пробелы для получения треугольной таблицы (т.е. для получения пустых ячеек).
num2str — преобразует числовой массив в массив строк.
strrep — заменяет в заданной строке одну подстроку на другую; используется для замены NaN на три пробела.
disp — выводит в командное окно.
fopen — открывает файл; функция fopen возвращает идентификатор файла.
fprintf — записывает данные в файл; для записи строк используется формат ‘%s’.
fclose — закрывает файл.

function DDIF(x, fun, fname)

%выводит таблицу разделенных разностей в командное окно и в файл

x - точки, в которых вычисляются разделенные разности
fun - имя или указатель на функцию
fname - имя файла

узнаем число точек
N = length(x);
вычисляем значения функции fun от x
y = feval(fun, x);
делаем из y столбец 1 на N
y = reshape(y, 1, N);
создаем массив 2*N-1 на N и заполняем его значениями NaN
D = NaN(2*N-1,  N);
заполняем первый столбец D значениями функции (через строку)
D(1:2:end, 1) = y;
проходим по столбцам массива D, в j-ый столбец записываем 
разделенную разность (j-1)-го порядка
for j = 2 : N
s = j - 1;
i = 0;
проходим по строкам массива D (через строку)
вычисляем разделенную разность через две предыдущие
for k = j : 2 : 2*N  - j 
i = i + 1;
D(k,j) = (D(k+1, j-1) - D(k-1,j-1)) / (x(i+s) - x(i));
end
end
преобразуем числовой массив D в массив строк STR
STR = num2str(D);
проходим по строкам STR и меняем NaN на три пробела
for i = 1: 2*N-1
STR(i,:) = strrep(STR(i,:), 'NaN', '   ');
end
выводим таблицу в командное окно
disp(STR)
открываем файл с именем, записанным в fname, для записи таблицы 
f = fopen(fname, 'wt+');
записываем каждую строку таблицы в файл
for i = 1 : 2*N-1
fprintf(f, '%s', STR(i,:));
end
% закрываем файл
fclose(f);

Для получения, например таблицы разделенных разностей функции sin x по шести равноотстоящим точкам на отрезке , достаточно сформировать массив точек при помощи функции linspace и вызвать функцию DIFF:

>> x = linspace(0, pi, 6);
>> DDIF(x, 'sin', 'difftable.txt')

В результате получаем в командном окне:

Та же самая таблица записана в файл difftable.txt.

Разделенные разности обладают рядом важных свойств. Одно из них — симметрия, т.е. разделенная разность

не меняет своего значения при любой перестановке . Это свойство следует, например, из представления разделенной разности в виде

которое доказывается по индукции (см., например, Н.С. Бахвалов, Н.П. Жидков, Г.М. Кобельков. Численные методы. Изд-во «Бином». 2004г.).

Получение интерполяционного полинома в форме Ньютона

Получим выражение для интерполяционного полинома в форме Ньютона для интерполяции функции f(x) по набору точек , воспользовавшись формой Лагранжа (вывод формы Лагранжа интерполяционного полинома приведен в разделе Интерполяционный полином в форме Лагранжа).

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

(очевидно, что в правой части все слагаемые, кроме сокращаются). Здесь — интерполяционный полином, построенный для f(x) по набору точек. Он записывается в форме Лагранжа следующим образом:

Ясно также, что

т.к. интерполяционный полином нулевой степени, построенный по одной точке, есть константа, значение которой равно

Наша цель — установить следующие равенства

где

Т.е., надо доказать, что

Тогда интерполяционный полином приобретает следующую форму

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

Преимущества записи интерполяционного полинома в форме Ньютона очевидны. Если в точке x вычислено значение интерполяционного полинома , построенного по набору точе , а затем к этому набору добавилась еще одна точка , то для получения значения нового интерполяционного полинома (степень которого на единицу выше) в точке достаточно добавить к значению слагаемое

Итак, докажем, что

для , где

Ясно, что разность является полиномом степени k+1. Далее, поскольку.Lk+1(x) — интерполяционный полином функции f(x), построенный по (k+2)-м точкам — интерполяционный полином функции f(x), построенный по (k+1)-ой точке , то

Тогда разность обращается в ноль в точках и так как она является полиномом степени k+1 то ее можно представить в виде

Осталось доказать, что

Для точки xk+2

и, следовательно, справедливо равенство

Запишем выражение для разности функции и интерполяционного полинома

и преобразуем его, вынося за скобки получим тогда

Воспользуемся теперь тем, что

Но поскольку разделенная разность есть симметричная функция своих аргументов (см.раздел Разделенные разности и некоторые их свойства )

.

Тогда очевидно, что

Подставим в это выражение вместо x значение xk+2, получим

Т.к. по доказанному выше

и

то

Но поскольку разделенная разность есть симметричная функция своих аргументов (см.раздел Разделенные разности и некоторые их свойства ), то

.

Таким образом мы установили, что для интерполяционного полинома функции степени n, построенного по точкам , справедливо представление в форме Ньютона:

Алгоритм вычисления интерполяционного полинома в форме Ньютона.

В предыдущем разделе Получение интерполяционного полинома в форме Ньютона была выведена формула для интерполяционного полинома, которую можно также записать в следующем более компактном виде:

.

где

Мы можем считать, что точки пронумерованы в обратном порядке, тогда

где

Тогда при вычислении значения интерполяционного полинома в форме Ньютона требуется выполнить следующие действия (см. книгу К. Де Бор. Практическое руководство по сплайнам. Изд-во «Радио и связь». Москва. 1985):

  1. найти разделенные разности, стоящие внизу треугольной таблицы (см. разд. Разделенные разности и некоторые их свойства);
  2. вычислить значение приведенного выше полинома от заданного значения x.

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

Алгоритм вычисления разделенных разностей.

Цикл по

Цикл по

После завершения цикла

При выполнении второго пункта, т.е. при вычислении интерполяционного полинома по формуле

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

сложений и

умножений.

Записав выражение для интерполяционного полинома в виде

можно найти значение интерполяционного полинома за n умножений и 2n сложений, т.е. за количество операций на порядок меньшее, чем при вычислении впрямую.

Для вычисления интерполяционного полинома в форме Ньютона в MATLAB можно воспользоваться приведенной ниже функцией newton.

function yy = newton(x, y, xx)
Вычисление интерполяционного полинома в форме Ньютона

x – массив с абсциссами точек, через которые должен проходить интерполяционный полином
y – массив ординат точек, через которые должен проходить интерполяционный полином
xx – массив значений независимой переменной, 
для которых надо вычислить интерполяционный полином
yy – вычисленные значения интерполяционного полинома

определяем число точек
N = length(x);
вычисляем разделенные разности
DIFF = y;
for  k = 1 : N-1
for i = 1: N - k
DIFF(i) = (DIFF(i+1) - DIFF(i)) / (x(i+k) - x(i));
end
end
вычисляем значения интерполяционного полинома в точках xx
с использованием операции поэлементного умножения .* 
для получения сразу всех значений полинома yy
yy = DIFF(1) * ones(size(xx));
for k = 2 : N
yy = DIFF(k) + (xx - x(k)) .* yy;
end

Например, интерполируем при помощи функции newton следующую табличную функцию, заданную массивами

x = [0.1 2.2 3.1 4.9 6.5] ;
y = [0.8 0.9 0.1 -0.7  -0.9];

при помощи интерполяционного полинома (4-ой степени) и вычислим значения интерполяционного полинома в 100 равноотстоящих точках на отрезке интерполирования [0.1, 6.5]

xx = linspace(x(1), x(end), 100);
yy = newton(x, y, xx);

Для проверки работы функции newton можно построить график, отобразив значения табличной функции маркерами, а значения интерполяционного полинома непрерывной линией

plot(x,y,'o',xx,yy)

В результате получаем

Интерполяция табличной функции полиномом

Поскольку интерполяционный полином существует и единственный (см. раздел Задача интерполяции функции, интерполяционные полиномы), то полином в форме Ньютона совпадает с интерполяционным полиномом в форме Лагранжа (с точностью до ошибок округления), построенным в разделе Интерполяционный полином в форме Лагранжа. Все результаты относительно поведения ошибки интерполяции и сходимости интерполяционного полинома, приведенные в разделах Ошибка интерполяции, чебышевские узлы и Сходимость интерполяционного полинома, остаются в силе.

Интерполяционные формулы для равноотстоящих узлов

Если узлы интерполяции отстоят друг от друга на одинаковое расстояние, т.е.

то для интерполяционных полиномов в форме Лагранжа и Ньютона можно получить специальные формулы, упрощающие вычисления (см. И.С. Березин, Н.П. Жидков. Методы вычислений. Т.1. Изд-во ФизМатЛит, 1962).

Рассмотрим, например интерполяционный полином в форме Лагранжа (см. раздел Интерполяционный полином в форме Лагранжа)

где

Введем обозначение

Тогда выражение для принимает вид

Проводя сокращения в числителе и знаменателе и предполагая, что , то соответствующий x(t) является одним из узлов интерполяции и значение интерполяционного по) получим, что

Тогда

В полученной записи интерполяционного полинома коэффициент

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

Получение интерполяционного полинома в символьном виде и его коэффициентов

Часто возникает задача не только вычислить значение интерполяционного полинома в заданных точках, но и получить его коэффициенты. Для решения этой задачи воспользуемся средствами Symbolic Math Toolbox. Приведенная ниже функция interpsymb по заданным массивам X и Y с интерполируемой табличной функцией выполняет следующие действия:

  1. строит интерполяционный полином в символьном виде;
  2. выводит в командное окно полином с коэффициентами, точно представленными рациональными числами;
  3. выводит в командное окно полином с коэффициентами, вычисленными с заданной во входном аргументе tol точностью;
  4. строит график полинома и визуализирует табличную функцию.
function L = interpsymb(X, Y, tol)
построение интерполяционного полинома в символьном виде
 
X  и Y – массивы, определяющие табличную функцию
tol – точность, с которой следует вычислить коэффициенты интерполяционного полинома. 
 
определяем символьные переменные
syms x L PHI
находим число точек, по которым надо построить интерполяционный полином
n = length(X);
вычисляем интерполяционный полином в форме Лагранжа
L = 0;
for k = 1:n
PHI = 1;
for i= [1:k-1 k+1:n]
PHI = PHI * (x-X(i)) / (X(k)-X(i));
end
L = L + Y(k)*PHI;
end
приводим подобные слагаемые
L = collect(L);
выводим в командное окно полином с рациональными коэффициентами
pretty(L)
вычисляем коэффициенты с заданной точностью
pretty(vpa(L, tol))
создаем графическое окно
figure('Color', 'w')
выводим график интерполяционного полинома
ezplot(L,[X(1) X(end)])
визуализируем табличную функцию
hold on
plot(X,Y,'r*')

При интерполяции, например следующих табличных данных

x = [1 2 4 5 6];
y = [1 -1 3 6 7];

при помощи функции interpsymb

L = interpsymb(x, y, 5);

получаем в командном окне

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

Табличная функция и интерполяционный полином


Поиск по сайту:

Система Orphus

Яндекс.Метрика