Алгоритмы и программы с интерполяцией
При обработке экспериментальных данных часто возникает необходимость в определении значения функции внутри значений исследуемых аргументов, но при несовпадении нужного аргумента с имеющимися. В математике данная процедура называется интерполяцией.
Существует много различных алгоритмов интерполяции. В настоящей статье нами рассматриваются два наиболее старых и известных метода: Интерполяция по Лагранжу и Ньютону.
Интерполяция по Лагранжу.
Сущность метода состоит в следующем. Имеется N пар значений аргумента и функции, при этом значения аргументов могут быть как равноотстоящие, так и не равноотстоящие. По имеющимся значениям строится алгебраический полином, степень которого на 1 меньше количества пар исходных заданных значений аргументов и функции. Т.е. при, например, 10 пар, строится полином девятой степени. Дальнейшие рассуждения следующие. Если построенный полином совпадает с исходной анализируемой функцией во всех имеющихся точках, то можно считать, что и в точках, между имеющимися, тоже будет значительная близость. Поэтому затем вычисляется значение найденного полинома в искомой точке.
Ниже приведены программы, реализующие данный алгоритм на
Далее приводится вариант программы на языке С (компилятор BORLAND C++ 3.1)
// LAGR.CPP
//*****************************************************************
// Программа интерполирования табличной функции по методу Лагранжа.
// Цикл ввода значения аргумента, для которого ищется значение
// функции - бесконечный. Выход по вводу значения вне введенных X.
//*****************************************************************
// Вводятся последовательно:
// - кол-во точек интерполяции
// - затем пары значений (абсцисса и ордината табличного значения)
// Внимание!!! Абсциссы не обязательно равноотстоящие.
// Контрольная точка: при вводе
// - 4
// - 0 -1 1 -3 2 3 6 1187
// - 4
// Получаю 255.0
#include <stdio.h>
#include <conio.h>
int n; // количество исходных точек интерполяции
float mx[20]; // массив исходных значений абсцисс
float my[20]; // аналогично для ординат исходной табличной функции
float x; // значение абсциссы, для которой ищется значение табличной функции
int i,j; // переменные цикла
float s; // сюда буду накапливать искомое интерполяционное значение
float q; // сюда буду накапливать очередной коэффициент Лагранжа
void main(void)
{ // 0
clrscr(); //начальная очистка экрана
puts("Введи количество исходных интерполируемых точек:");
scanf("%d",&n);
puts("Введи попарно сами исходные абсциссы (от меньших к большим) и ординаты:");
// Ввод самих исходных значений
for(i=0; i<n; scanf("%f%f",&mx[i++],&my[i]));
// бесконечный внешний цикл
while (1)
{
puts("Введи значение абсциссы, для которой надо найти значение ординаты:");
scanf("%f",&x);
// Выход из цикла при введе значения вне исходного интервала
if (x<mx[0] || x>mx[n-1]) break;
// Собственно получение искомого значения
for(s=0, j=0; j<n; s+=q*my[j++])
for(i=0, q=1; i<n;
q*=((j!=i)?(x-mx[i])/(mx[j]-mx[i]):(1)), i++);
printf("Полученное значение равно %10.3fn",s);
}
} // -0
Далее приводится текст программы на языке PASCAL. Проверено на компиляторах
PASCALABC 3.1 и PASCALABC.NET 2.2.
// LAGR.PAS
//*****************************************************************
// Программа интерполирования табличной функции по методу Лагранжа.
// Цикл ввода значения аргумента, для которого ищется значение
// функции - бесконечный. Выход по вводу значения вне введенных X.
//*****************************************************************
// Вводятся последовательно:
// - кол-во точек интерполяции
// - затем пары значений (абсцисса и ордината табличного значения)
// Внимание!!! Абсциссы не обязательно равноотстоящие.
// Контрольная точка: при вводе
// - 4
// - 0 -1 1 -3 2 3 6 1187
// - 4
// Получаю 255.0
var
n:integer; // количество исходных точек интерполяции
mx:array [0..19] of real; // массив исходных значений абсцисс
my:array [0..19] of real; // аналогично для ординат исходной табличной функции
x:real; // значение абсциссы, для которой ищется значение табличной функции
i,j:integer; // переменные цикла
s:real; // сюда буду накапливать искомое интерполяционное значение
q:real; // сюда буду накапливать очередной коэффициент Лагранжа
begin // 0
// clrscr(); //начальная очистка экрана
writeln('Введи количество исходных интерполируемых точек:');
readln(n);
writeln('Введи попарно сами исходные абсциссы (от меньших к большим) и ординаты:');
// Ввод самих исходных значений
for i:=0 to n-1 do
read(mx[i],my[i]);
// бесконечный внешний цикл
while (1>0) do
begin
write('Введи значение абсциссы, для которой надо найти значение ординаты:');
readln(x);
// Выход из цикла при введении значения вне исходного интервала
if ((x<mx[0]) or (x>mx[n-1])) then break;
// Собственно получение искомого значения
s:=0;
for j:=0 to n-1 do
begin
q:=1;
for i:=0 to n-1 do
if (i<>j) then q:=q*(x-mx[i])/(mx[j]-mx[i]);
s:=s+q*my[j];
end;
writeln('Полученное значение равно ',s);
end;
end. // -0
Далее приводится текс программы на языке QBASIC в компиляторе
MS-DOS QBASIC 1.0.
' LAGR.BAS '***************************************************************** ' Программа интерполирования табличной функции по методу Лагранжа. ' Цикл ввода значения аргумента, для которого ищется значение ' функции - бесконечный. Выход по вводу значения вне введенных X. '***************************************************************** ' Вводятся последовательно: ' - кол-во точек интерполяции ' - затем пары значений (абсцисса и ордината табличного значения) ' Внимание!!! Абсциссы не обязательно равноотстоящие. ' При вводе: ' 4 ' 0, -1 ' 1, -3 ' 2, 3 ' 6, 1187 ' 4 ' Получаю 255.0 DIM n AS INTEGER ' количество исходных точек интерполяции INPUT "Введите количество исходных пар точек ", n DIM mx(0 TO n - 1) AS SINGLE' массив исходных значений абсцисс DIM my(0 TO n - 1) AS SINGLE' аналогично для ординат исходной табличной функции DIM x AS SINGLE ' значение абсциссы, для которой ищется значение табличной функции DIM i AS INTEGER ' переменные цикла DIM j AS INTEGER ' переменные цикла DIM s AS SINGLE ' сюда буду накапливать искомое интерполяционное значение DIM q AS SINGLE ' сюда буду накапливать очередной коэффициент Лагранжа CLS 'начальная очистка экрана PRINT "Введи попарно сами исходные абсциссы (от меньших к большим) и ординаты: "; ' Ввод самих исходных значений FOR i = 0 TO n - 1 INPUT mx(i), my(i) NEXT i ' бесконечный внешний цикл DO WHILE (1 > 0) INPUT "Введи значение абсциссы, для которой надо найти значение ординаты: ", x ' Выход из цикла при введении значения вне исходного интервала IF ((x < mx(0)) OR (x > mx(n - 1))) THEN EXIT DO ' Собственно получение искомого значения s = 0 FOR j = 0 TO n - 1 q = 1 FOR i = 0 TO n - 1 IF (i <> j) THEN q = q * (x - mx(i)) / (mx(j) - mx(i)) NEXT i s = s + q * my(j) NEXT j PRINT "Полученное значение равно ", s LOOP END
Ниже приводится текст этой же программы составленный для интерпретатора Python 3.4.3.
# LAGR.PY
#*****************************************************************
# Программа интерполирования табличной функции по методу Лагранжа.
# Цикл ввода значения аргумента, для которого ищется значение
# функции - бесконечный. Выход по вводу значения вне введенных X.
#*****************************************************************
# Вводятся последовательно:
# - кол-во точек интерполяции
# - затем пары значений (абсцисса и ордината табличного значения)
# Внимание!!! Абсциссы не обязательно равноотстоящие.
# При вводе:
# 4
# 0, -1
# 1, -3
# 2, 3
# 6, 1187
# 4
# Получаю 255.0
#import pdb; pdb.set_trace()
n= int(input('Введи количество исходных точек '))
mx = [] # объявляю пустой список (в данном случае одномерный массив)
my = [] # объявляю пустой список (в данном случае одномерный массив)
print('Введи попарно сами исходные абсциссы (от меньших к большим) и ординаты: ')
# Заполнение исходных массивов
for i in range(n):
(xm,ym)= list(map(float, input().split()))
mx.append(xm)
my.append(ym)
# бесконечный внешний цикл
while(1>0):
x=float(input('Введи значение абсциссы, для которой надо найти значение ординаты: '))
# Выход из цикла при введении значения вне исходного интервала
if ((x < mx[0]) or (x > mx[n - 1])): break
# Собственно получение искомого значения
s = 0
for j in range(n):
q = 1
for i in range(n):
if (i!=j): q = q * (x - mx[i]) / (mx[j] - mx[i])
s = s + q * my[j]
print('Полученное значение равно ', s)
input() # команда ДЕРЖИТ экран
Интерполяция по Ньютону
Другим, часто используемым алгоритмом интерполяции, является интерполяция по Ньютону. Данный метод может использоваться как для равноотстоящих узлов, так и для неравномерно расположенных. Нами реализован алгоритм для равноотстоящих узлов. Суть алгоритма состоит в том, что значение функции в любой точке внутри заданного диапазона определяется не только самими исходными значениями функции, но и так называемыми конечными разностями различных порядков. Конечные разности первого порядка - это разность двух рядом стоящих значений исходной функции. Конечная разность второго порядка - это разность двух рядом стоящих значений конечных разностей первого порядка и так далее.
Ниже приведены программы, реализующие данный алгоритм на трех языках.
Далее приводится вариант программы на языке С (компилятор BORLAND C++ 3.1)
// NEWT.CPP
//*******************************************************
// Программа интерполяции по НЬЮТОНУ.
// Реализован алгоритм равноотстоящих аргументов.
// При работе с уже имеющимися в программе данными
// надо ввести 7 - кол-во исходных пар точек
// 2.05 - абсцисса, для которой ищу
// Программа выдает 0.04877852.
//*******************************************************
#include <conio.h>
#include <stdio.h>
int main(void)
{
float x,y,s,q;
int i,j,n;
float mx[20];
float my[20];
float dy[20][20];
printf("Введите количество исходных пар точек ");
scanf("%d",&n);
//PRINT "Введи исходные пары X и Y через запятую в строку "
//FOR i = 0 TO n-1
// INPUT mx[i], my[i]
//NEXT
mx[0]=2; mx[1]=2.1; mx[2]=2.2; mx[3]=2.3;
mx[4]=2.4; mx[5]=2.5; mx[6]=2.6;
my[0]=0.054; my[1]=0.044; my[2]=0.0355; my[3]=0.0283;
my[4]=0.0224; my[5]=0.0175; my[6]=0.0136;
// Перепись массива my в первый столбец массива dy.
for( i=0;i<n;i++)
dy[i][ 0]=my[i];
// Заполнение массива конечных разностей
for(j=0;j<= n - 2;j++)
for(i=0;i<=n - j - 2; i++)
dy[i][ j + 1]=dy[i + 1][ j] - dy[i][ j];
printf("Введи значение абсциссы, для которой надо найти значение ординаты: ");
scanf("%f",&x);
// Собственно получение искомого значения
s=dy[0][ 0];
y=q=(x - mx[0]) / (mx[1] - mx[0]);
for( j=1 ;j<=n - 2;j++)
{
y/=j;
s+=dy[0][ j] * y;
y*=(q - j);
}
printf("Полученное значение равно %10.8f", s);
}
Далее приводится текст программы на языке PASCAL. Проверено на компиляторах
PASCALABC 3.1 и PASCALABC.NET 2.2.
// NEWT.PAS
//*******************************************************
// Программа интерполяции по НЬЮТОНУ.
// Реализован алгоритм равноотстоящих аргументов.
// При работе с уже имеющимися в программе данными
// надо ввести 7 - кол-во исходных пар точек
// 2.05 - абсцисса, для которой ищу
// Программа выдает 0.04877852.
//*******************************************************
var x,y,s,q:real;
i,j,n:integer;
mx: array[0..20] of real;
my: array[0..20] of real;
dy: array[0..20,0..20] of real;
begin
write('Введите количество исходных пар точек ');
readln(n);
//PRINT "Введи исходные пары X и Y через запятую в строку "
//FOR i := 0 TO n-1
// INPUT mx[i], my[i]
//NEXT
mx[0]:=2; mx[1]:=2.1; mx[2]:=2.2; mx[3]:=2.3;
mx[4]:=2.4; mx[5]:=2.5; mx[6]:=2.6;
my[0]:=0.054; my[1]:=0.044; my[2]:=0.0355; my[3]:=0.0283;
my[4]:=0.0224; my[5]:=0.0175; my[6]:=0.0136;
// Перепись массива my в первый столбец массива dy.
FOR i:=0 TO n - 1 do
dy[i, 0]:=my[i];
// Заполнение массива конечных разностей
FOR j:=0 TO n - 2 do
FOR i:=0 TO n - j - 2 do
dy[i, j + 1]:=dy[i + 1, j] - dy[i, j];
write('Введи значение абсциссы, для которой надо найти значение ординаты: ');
readln(x);
// Собственно получение искомого значения
s:=dy[0, 0];
q:=(x - mx[0]) / (mx[1] - mx[0]);
y:=q;
FOR j:=1 TO n - 2 do
begin
y:=y / j;
s:=s + dy[0, j] * y;
y:=y * (q - j);
end;
writeln('Полученное значение равно ', s:10:8);
END.
Далее приводится текс программы на языке QBASIC в компиляторе
MS-DOS QBASIC 1.0.
' NEWT.BAS '******************************************************* ' Программа интерполяции по НЬЮТОНУ. ' Реализован алгоритм равноотстоящих аргументов. ' При работе с уже имеющимися в программе данными ' надо ввести 7 - кол-во исходных пар точек ' 2.05 - абсцисса, для которой ищу ' Программа выдает 0.04877852. '******************************************************* DIM x AS SINGLE, y AS SINGLE, s AS SINGLE, q AS SINGLE DIM i AS INTEGER, j AS INTEGER INPUT "Введите количество исходных пар точек ", n DIM mx(0 TO n - 1) AS SINGLE DIM my(0 TO n - 1) AS SINGLE DIM dy(0 TO n - 1, 0 TO n - 1) AS SINGLE 'PRINT "Введи исходные пары X и Y через запятую в строку " 'FOR i = 0 TO n-1 ' INPUT mx(i), my(i) 'NEXT mx(0) = 2!: mx(1) = 2.1: mx(2) = 2.2: mx(3) = 2.3: mx(4) = 2.4: mx(5) = 2.5: mx(6) = 2.6 my(0) = .054: my(1) = .044: my(2) = .0355: my(3) = .0283 my(4) = .0224: my(5) = .0175: my(6) = .0136 ' Перепись массива my в первый столбец массива dy. FOR i = 0 TO n - 1 dy(i, 0) = my(i) NEXT ' Заполнение массива конечных разностей FOR j = 0 TO n - 2 FOR i = 0 TO n - j - 2 dy(i, j + 1) = dy(i + 1, j) - dy(i, j) NEXT i, j INPUT "Введи значение абсциссы, для которой надо найти значение ординаты:", x ' Собственно получение искомого значения s = dy(0, 0) q = (x - mx(0)) / (mx(1) - mx(0)) y = q FOR j = 1 TO n - 2 y = y / j s = s + dy(0, j) * y y = y * (q - j) NEXT PRINT "Полученное значение равно ", PRINT USING "##.########"; s END