Простейшие способы интерполяции и экстраполяции данных

Рассмотрим функцию y=f(x) непрерывную на некотором интервале [a,b] и заданную некоторыми своими значениями y_i=f(x_i), i=1,2,…,n для соответствующих значений аргумента a≤x₀<x₁<. . . <x_n ≤b. Необходимо найти значение этой функции в точке x* принадлежащей интервалу [a,b], при этом x*≠x_i и оценить погрешность полученного приближенного значения.

Нахождение неизвестных значений x* внутри отрезка [a,b] называется интерполированием.

Один из возможных путей решения поставленной задачи состоит в следующем. Для функции y=f(x) по значениям y_i в узлах x_i, i=1,2,…,n строится многочлен степени не выше n

, (1)

принимающий в точках x_i значения y_i, т.е. значения коэффициентов многочлена – a_i – находятся из условия:

. (2)

Этот многочлен называется интерполяционным, он всегда существует и единственен.

В простейшем случае можно воспользоваться многочленом 1-ой степени, т.е. интерполирование будет линейным:

. (3)

В этом случае для интерполяции необходимо построить прямую, соединяющую две точки: (x_j,y_j) и (x_j+1,y_j+1), 1≤j+1≤n. Коэффициенты a и b вычисляются следующим образом:

;

. (4)

Описание примера алгоритма линейной интерполяции:

Задаем начальные данные для интерполяции: границы интерполируемого отрезка a и b, количество узлов координатной сетки n, шаг сетки h, массив значений координат сетки z(n) и массив значений интерполируемой функции f(n).

Задаем шаг графического изображения sh интерполируемой функции f(x) и интерполяционной функции g(x).

Задаем цикл по х от a до b-sh с шагом sh*2.

Задаем x₁=x и x₂=x+sh.

Вычисляем значения интерполируемой функции в точках x₁ и x₂: f(x₁)и f(x₂).

Выводим на экран изображение интерполируемой функции в точках x₁ и x₂: f(x₁)и f(x₂).

Вычисляем значения интерполяционной функции в точках x₁ и x₂: g(x₁)и g(x₂):

Цикл по k от 1 до n-1.

Находим номера узлов интерполяционной сетки, между которыми находится значение x₁. Задаем ветвление: если x>=z(k) и x<=z(k+1), то k₁=k и k₂=k+1;если x>z(n), то k₁=n-1 и k₂=n.

По формулам (4) метода линейной интерполяции находим коэффициенты интерполяционной прямой: a=(f(k₁)-f(k₂))/(z(k₁)-z(k₂)) и b=f(k₂)-a*z(k₂).

Вычисляем значение интерполяционной функции в точке x₁: g(x₁).

Вычисляем значение интерполяционной функции в точке x₂: g(x₂).

Конец цикла по k.

Выводим на экран интерполяционную функцию в точках x₁и x₂: g(x₁)и g(x₂).

Конец цикла по x.

Конец программы.

Пример алгоритма линейной интерполяции на естественном языке:

n =количество_узлов_сетки

Определить Массив f(n), z(n)

h = шаг_сетки

Цикл по i от 1 до n

z(i) = (i - 1) * h

x = z(i)

f(i) = Sin(5 * x)

Конец Цикла

sh = 0.005

x1 = 0.05

f1 = Sin(5 * x1)

g1 = Интерполяция(x1, f, z, n)

Печать x1, f1, g1

i = 1

Цикл по x от 0.05 до 0.6 с шагом sh

i = i + 1

x2 = x + sh

f2 = Sin(5 * x2)

g2 = Интерполяция(x2, f, z, n)

Печать x2, f2, g2

x1 = x2

f1 = f2

g1 = g2

Конец Цикла

Определение Функции Интерполяция

Параметры t, f, z, n

Локальные_переменные k, a, b

Цикл по k от 1 до n - 1

Если t >= z(k) И t <= z(k + 1) Тогда

k1 = k

k2 = k + 1

Выход Из Цикла

Иначе Если t > z(n) Тогда

k1 = n - 1

k2 = n

Выход Из Цикла

Конец Если

Конец Цикла

a = (f(k1) - f(k2)) / (z(k1) - z(k2))

b = f(k2) - a * z(k2)

Возврат a * t + b

Конец Определения Функции

Пример реализации алгоритма линейной интерполяции функции f(x)=sin(5x) на отрезке от 0,05 до 5 с шагом 0,005 на языке VFP:

SET DECIMALS TO 10

n=10

DIMENSION f(n),z(n)

h=0.15

FOR i=1 TO n

z(i)=(i-1)*h

x=z(i)

f(i)= SIN(5*x)

ENDFOR

a=0.05

b=2

sh = 0.005

x1=0.05

_screen.Cls

FOR x=a TO b-sh STEP sh*2

x1=x

x2=x+sh

f1 = SIN(5*x1)

f2 = SIN(5*x2)

_screen.drawWidth=2

_screen.Line(xe(x1),ye(f1),xe(x2),ye(f2))

g1=интерполяция(x1,@f,@z,n)

g2=интерполяция(x2,@f,@z,n)

_screen.drawWidth=1

_screen.line(xe(x1),ye(g1),xe(x2),ye(g2))

WAIT WINDOW 1

ENDFOR

PROCEDURE интерполяция

PARAMETERS x,f,z,n

DIMENSION f(n),z(n)

LOCAL k,a,b

FOR k=1 TO n-1

DO CASE

CASE x>=z(k) and x<=z(k+1)

k1=k

k2=k+1

EXIT

CASE x>z(n)

k1=n-1

k2=n

EXIT

ENDCASE

ENDFOR

a=(f(k1)-f(k2))/(z(k1)-z(k2))

b=f(k2)-a*z(k2)

RETURN a*x+b

FUNCTION xe

PARAMETERS x

RETURN 50 + x*400

ENDFUNC

FUNCTION ye

PARAMETERS y

RETURN 200 - y*150

ENDFUNC

Результат работы программы (фрагмент изображение основного окна разработки VFP) приведен на рисунке 1.:

Рис. 1

На рисунке 1 хорошо видны интерполируемая функция (более жирная) и интерполяционная функции.

Реализация алгоритма линейной интерполяции функции f(x)=sin(5x) на отрезке от 0,05 до 5 с шагом 0,005 на языке VBA (для макроса MS EXCEL):

Sub Interpolation_Test()

Dim xlSheet As Excel.Worksheet

Set xlSheet = Application.ActiveSheet

n = 10

Dim f(10), z(10)

h = 0.1

For i = 1 To n

z(i) = (i - 1) * h

x = z(i)

f(i) = Sin(5 * x)

Next i

sh = 0.005

x1 = 0.05

f1 = Sin(5 * x1)

g1 = Interpolation(x1, f, z, n)

xlSheet.Cells(1, 1).Value = x1

xlSheet.Cells(1, 2).Value = f1

xlSheet.Cells(1, 3).Value = g1

i = 1

For x = 0.05 To 0.6 Step sh

i = i + 1

x2 = x + sh

f2 = Sin(5 * x2)

g2 = Interpolation(x2, f, z, n)

xlSheet.Cells(i, 1).Value = x2

xlSheet.Cells(i, 2).Value = f2

xlSheet.Cells(i, 3).Value = g2

x1 = x2

f1 = f2

g1 = g2

Next x

Call AddChart("C1", "C2", "C3")

End Sub

Function Interpolation(t, f, z, n)

Dim k, a, b

For k = 1 To n - 1

If t >= z(k) And t <= z(k + 1) Then

k1 = k

k2 = k + 1

Exit For

ElseIf t > z(n) Then

k1 = n - 1

k2 = n

Exit For

End If

Next k

a = (f(k1) - f(k2)) / (z(k1) - z(k2))

b = f(k2) - a * z(k2)

Interpolation = a * t + b

End Function

Sub AddChart(Col1, Col2, Col3)

ASheetName = Application.ActiveSheet.Name

Charts.Add

ActiveChart.ChartType = xlXYScatterSmoothNoMarkers

ActiveChart.SetSourceData Source:=Sheets(ASheetName).Range("A1")

ActiveChart.SeriesCollection.NewSeries

ActiveChart.SeriesCollection.NewSeries

ActiveChart.SeriesCollection(1).XValues = "=" & ASheetName & "!" & Col1

ActiveChart.SeriesCollection(1).Values = "=" & ASheetName & "!" & Col2

ActiveChart.SeriesCollection(1).Name = "=""f(x)"""

ActiveChart.SeriesCollection(2).XValues = "=" & ASheetName & "!" & Col1

ActiveChart.SeriesCollection(2).Values = "=" & ASheetName & "!" & Col3

ActiveChart.SeriesCollection(2).Name = "=""g(x)"""

ActiveChart.Location Where:=xlLocationAsObject, Name:=ASheetName

ActiveChart.SeriesCollection(1).Select

ActiveChart.ChartArea.Select

End Sub

В данной реализации для графического изображения интерполируемой функции f(x) и интерполяционной функции g(x)используются средства построения графиков MS EXCEL. Значения x , f(x) и g(x) заносятся в ячейки первых трех колонок рабочего листа:

xlSheet.Cells(i, 1).Value = x1

xlSheet.Cells(i, 2).Value = f1

xlSheet.Cells(i, 3).Value = g1

и в конце программы для построения графика вызывается метод AddChart(Col1, Col2, Col3).

Результаты работы программы на VBA приведен на рисунке 2.

x	f(x)	g(x)
0,05	0,247404	0,239713
0,1	0,479426	0,479426
0,15	0,681639	0,660448
0,2	0,841471	0,841471
0,25	0,948985	0,919483
0,3	0,997495	0,997495
0,35	0,983986	0,953396
0,4	0,909297	0,909297
0,45	0,778073	0,753885
0,5	0,598472	0,598472
0,55	0,381661	0,369796
0,6	0,14112	0,14112
0,65	-0,1082	-0,10483

Рис. 2