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

5.1.Цель работы

Познакомиться с методами решения задачи интерполяции дискретных данных. Научиться пользоваться операторами интерполяции программы Mathcad.

5.2.Задание

В промежутке [-1, 1] по оси x задано(см.Таблицу 5.1) три значения аргумента и соответствующие им значения функции y_i=f(x), где i - номер варианта задания.

По исходным данным в таблице 5.1 нужно определить y_i при двух промежуточных значениях x =-0.5 и x =0.5 c точностью до двух знаков после запятой.

Более сложные задания для 6 исходных точек представлены в таблице 5.2. В промежутке [0, 10] задано 6 дискретных значений аргумента с соответствующими значениями функции y_i=f(x), где i- номер варианта задания.

По исходным данным в промежутке [0, 10] из таблицы 5.2 необходимо определить 5 промежуточных значений функции при указанных значениях аргумента с точностью до одного знака после запятой.

Для выполнения расчётных работ предлагаются следующие варианты заданий.

1.Разработать программу на языке программирования Mathcad для кусочно-линейной интерполяции и решить одну из тестовых задач.

2.Разработать программу на языке программирования Mathcad для кусочно-квадратичной интерполяции и решить одну из тестовых задач.

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

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

5.3.Теоретические сведения

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

Решение задачи интерполирования [1-3] не является однозначным. Решение может быть получено при помощи различного вида функций. Можно соединить соседние точки прямыми линиями и использовать их для расчета промежуточных значений. Такая методика называется кусочно-линейной интерполяцией. Если весь промежуток интерполяции разделить на трёхточечные интервалы, то можно выполнить кусочно-квадратичную интерполяцию.

Разрывы непрерывности производных в узловых точках можно устранить при помощи сплайн-интерполяции. Кубическая интерполяция позволяет построить кривую, проходящую через заданные точки, таким образом, чтобы первые и вторые производные были бы непрерывными во всех точках. Такая кривая может быть построена разделением промежутка интерполяции на трехточечные интервалы, внутри которых записываются полиномы третьей степени, проходящие через эти точки. Эти полиномы затем соединяются вместе и образуют общую кривую.

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

в узлах интерполяции x₀, x₁,...x_n должен точно совпадать со значениями функции y₀, y₁,...y_n. Как известно [1], существует единственный полином степени не выше n, принимающий в точках x₀, x₁...x_n заданные значения. Поэтому m=n, при этом число узлов интерполирования на единицу больше, чем степень полинома. Степень полинома равна числу интервалов, на которые разделяется весь промежуток интерполирования. Коэффициенты a_i полинома Q_n(x) можно определить различными способами.

1.Из решения системы уравнений

(5.1)

Порядок системы (5.1) равен числу узлов интерполирования, т.е. n+1. Решение системы (5.1) может быть получено способом Крамера с помощью определителей или способом Гаусса путем исключения переменных.

2.По формуле Лагранжа[1], позволяющей определить полином Q_n(x) при произвольных узлах интерполяции.

Q_n(x)= (5.2)

3.По формуле Ньютона [1], позволяющей определить полином Q_n(x) при равноотстоящих с шагом h узлах интерполирования.

(5.3)

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

В качестве примера рассмотрим решение первого задания из таблицы 5.1. Построим квадратичный полином на промежутке по x от -1 до 1 в трех узлах интерполирования: x0=-1, y0=1/3; x1=0, y1=1; x2=1, y2=3.

Составим систему уравнений (5.1)

Последовательным исключением переменных получим

При этом интерполирующий полином примет вид

Используя формулу Лагранжа(5.2), получим тот же полином

По формуле Ньютона(5.3) при h=1 получим

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

Рассмотрим в качестве примера один из возможных алгоритмов интерполяции с применением формулы Лагранжа. Представим формулу (5.2) в виде

Алгоритм интерполяции по формуле Лагранжа можно представить в следующем виде:

1.Зададим степень полинома n, значения функции y_i и независимой переменной x_i в узлах интерполирования i=0, 1,..n.

2.Зададим x, для которого необходимо найти значение полинома y=f(x).

3.Присвоим начальные значения коэффициентам a_i=b_i=1.

4.Рассчитаем коэффициенты a_i и b_i по формулам (5.5) и (5.6).

5.Определим значение полинома по формуле (5.4).

6.Выведем полученное значение полинома в качестве значения функции y=f(x) в заданной точке.

Схема алгоритма показана на рис.5.1.

Рис.5.1. Схема вычисления интерполяционного полинома

по формуле Лагранжа

При небольшом числе узлов алгоритм (рис.5.1) может быть реализован вручную. Но для большого количества узлов лучше записать алгоритм в виде программы на одном из языков программирования.

На рис.5.2 на языке Паскаль записана программа, реализующая алгоритм (рис.5.1). Программа рассчитывает значение полинома n- порядка в произвольных точках по n+1 значениям x и y.

В программе, использующей метку, рассчитывается функция Лагранжа L. Исходные данные: n - порядок полинома, x_i и y_i - исходные точки при i=0, 1, 2...n, промежуточное значение x. Программа выводит полученное значение y вместе заданным значением x. Её можно использовать в качестве подпрограммы в сложных программах математического моделирования, когда требуется интерполяция дискретных данных.

Подпрограмма на языке программирования Mathcad (рис.5.3), реализующая тот же алгоритм при трехточечных исходных данных оформлена в виде подпрограммы-скаляра, возвращающего значение интерполяционного полинома при произвольном значении аргумента z.

Для использования подпрограммы в документе Mathcad в качестве глобальных переменных нужно задать трёхточечные векторы исходных данных. Два вложенных цикла for и опрератор ветвления if позволяют рассчитать коэффициенты a_i и b_i в элементах суммы функции Лагранжа, а также рассчитать её значение.

Вывод полученного значения функции осуществляется в основном документе Маthcad.

Program Lagrang;

{Программа рассчитывает значение полинома y=a0+a1.x+a2.x^2+...+an.x^n

в произвольной точке x по n+1 дискретным значениям y и x }

type m=array[1..(n+1)] of real;

var y, x: real; yd, xd: m; i, n: integer;

function L(x: real): real;

label 1;

var i, j: integer; c: real; a, b: m;

begin for i: =1 to (n+1) do

begin a[i]: =1; b[i]: =1

end; c: =0; for i: =1 to (n+1) do

begin for j: =1 to (n+1) do

begin if i=j then goto 1; a[i]: =a[i]*(x-xd[i]);

b[j]: =b[j]*(xd[i]-xd[j]);

1: end; c: =c+yd[i]*a[i]/b[i]

end; L: =c

end;

begin writeln(' Введите порядок полинома n');

writeln(' Введите " n+1" пар дискретных значений x и y');

for i: =1 to n do

begin read(xd[i], yd[i])

end; writeln(' Введите промежуточное значение x');

read(x); y: =L(x); writeln(x, y)

end.

Рис.5.2. Программа на языке Паскаль для

интерполяции по формуле Лагранжа.

Рис.5.3. Подпрограмма на языке программирования Маthcad для интерполяции по формуле Лагранжа.

5.4. Порядок выполнения работы

Исходные данные для расчета представить в виде двух векторов x и y и графика с кусочно-линейной интерполяцией.

Решить задачу интерполяции, используя формулу Лагранжа. Воспользоваться программой для расчета значений интерполяционного полинома и рассчитать промежуточные значения в таблицах.

Определить промежуточные значения при помощи сплайн-интерполяции.

Построить графики интерполяции: кусочно-линейной, сплайном и обычным полиномом.

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

5.5.Пример выполнения работы

. Рассмотрим выполнение лабораторной работы на примере первого задания из таблицы 5.1.Выполнение расчетов будем производить при помощи программы Маthcad. Исходные данные для расчета.

Воспользуемся графическими возможностями Маthcad для визуализации исходных данных.По оси х отложим дискретные значения x_i, а по оси y соответственно y_i.

Рис.5.4. Кусочно-линейная графическая интерполяция исходных данных средствами Mathcad.

При построении графика для дискретных значений x_i и y_i Mathcad автоматически осуществляет их кусочно-линейную интерполяцию, результаты которой показаны на рис.5.4.

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

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

Запишем имеющиеся данные при помощи двух пар векторов. Векторы z₁ и f₁ соответствуют исходным данным, а векторы z и f полному набору точек с учетом полученных промежуточных значений.

Представим записанные векторы в графическом виде с использованием автоматической кусочно-линейной интерполяции Маthcad.

Рис.5.5. Кусочно-линейная графическая интерполяция данных:
f1k - исходные, fm - дополненные расчетом по формуле Лагранжа.

На рис.5.5 заметно различие значений функции в промежуточных точках, полученных интерполяцией обычными полиномами и интерполяцией кусочно-линейными функциями.

Интерполяционный полином Лагранжа в трехточечной схеме деления промежутка интерполирования можно записать в виде (5.4)

После подстановки исходного вектора данных: x₀=-1, x₁=0, x₂=1 и y₀=0.333, y₁=1, y₂=3 интерполяционный полином примет следующий вид:

Используя символьные преобразования Mathcad, получим

(5.7)

Значения этого полинома при zx=-0.5 и zx=0.5

cоответствуют полученным ранее при помощи подпрограммы.

Интеполяционная кривая, полученная при помощи обычного полинома, показана на рис.5.6.

В Mathcad имеются операторы автоматической интерполяции дискретных данных: linterp и interp.

Оператор linterp(v_x, v_y, x) осуществляет кусочно-линейную интерполяцию векторов исходных данных v_x и v_y
для точки x.

Применяя его для нашей задачи, получим

Непрерывная кривая

для кусочно-линейной интерполяции показана на рис.5.6.

Оператор inter(v_s, v_x, v_y, x) возвращает интерполированное кубическим сплайном значение функции в точке x.

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

Вектор v_s представляет собой вторые производные в заданных точках.Он возвращается оператором cspline(v_x, v_y). Применим оператор interp и cspline для нашей задачи в точках
x=-0.5 и x=0.5.

Применение сплайн-интерполяции во всем промежутке дает кривую,

,

показанную на рис.5.6.

Рис.5.6. Интерполяция данных обычным полиномом Lagr(zx),

сплайн-интерполяцией Intrp(zx),

кусочно-линейной функцией Linterp(zx).

Выводы

Пользуясь различными методами интерполяции, мы получили три разных набора данных, представленных в таблице 5.3.

Таблица 5.3