Выбор и построение интерполирующей функции (стр. 3 из 3)
(D)Найдем выражения для производных от сплайна S(i)(x):
(Е) 
(F)и подставим их в выражения (С) и (D). Как следствие, имеем
6 
(G) 
(H)Для получения еще двух необходимых уравнений воспользуемся условиями в конечных узлах. Например, можно считать концы линейки отпущенными, что отвечает их нулевой кривизне, то есть
(I) 
(J)Построенные при таких условиях кубические сплайны называют свободными. При наличии других известных асимптотических данных задачи, возможны и другие условия на концах отрезков.
Уравнения (A), (B), (G)-(J) составляют полную СЛАУ для определения 4N неизвестных коэффициентов. Если эту СЛАУ преобразовать, то ее решение значительно упростится.
Очевидно,
. Кроме того, из выражения (J) 
(K)а из выражения (H) –
(L)Подставив уравнение (L) в формулу (В) учитывая, что
, получим 
; (М) 
(N)Извлекая из (G) bi и bi+1 с помощью (М), а di – на основании (L), придем к такой СЛАУ относительно ci:
(**) 
Матрица этой тридиагональной, то есть нулю не равны только елементы главной и двух соседних диагоналей. Для ее решения можно воспользоваться любым методом, после чего надо найти bi и di из выражений (К) – (N).
Вообще-то можно рассмотреть задачу о нахождении сплайна n-й степени:
коэффициенты которого кусочно-постоянные и который в узлах интерполяции принимает значения заданной функции и непрерывный вместе со своими n-1 производными.
Практическая реализация
Программа на языке Pascal
В процессе выполнения работы мною была написана программа EITKIN на языке Pascal.
В данной программе есть два массива: одномерный массив X, в нем хранятся значения узлов интерполирования хi и двумерный массив Р, в нем хранятся значения многочленов степени не выше n, переменная z это, то значение для которого надо найти значение функции, n – количество узлов интерполирования. Все вычисления проводятся в одном встроенном цикле. Данные на экран выводятся в виде двухмерной матрицы.
Код программы:
program EITKIN;
uses wincrt, strings;
var x:array [1..60]of real;
p:array [1..60,1..60] of real;
z :real; i,j,n: integer;
begin
StrCopy(WindowTitle, 'Программа интерполяции функции по схеме Эйткина ');
clrscr;
write ('vvedite k-vo uzlov interpolirovanija n=');
readln (n);
write ('vvedite X dlja kotorogo nado najti znach func=');
readln (z);
writeln ('vvedite mas Xi');
for i:=1 to n do
begin
write ('vvedite elem X[',i,']=');
readln (x[i]);
end;
writeln ('vvedite mas Yi');
for i:=1 to n do
begin
write ('vvedite elem Y[',i,']=');
readln (p[1,i]);
end;
writeln ('PROCES VICHISLENIJA......');
for i:=2 to n do
begin
for j:=1 to n+1-i do
begin
p[i,j]:=1/(x[j+i-1]-x[j])*(p[i-1,j]*(x[j+i-1]-z)-p[i-1,j+1]*(x[j]-z));
end;
end;
writeln ('REZ MATRICA::::');
for i:=1 to n do
begin
write ('P^',i,'(',z:4:5,') | ');
for j:=1 to n+1-i do
begin
write (p[i,j]:4:5,' | ');
end;
writeln;
end;
writeln ('!!!!!!!!!OTVET!!!!!!!!!');
writeln ('y(',z:4:5,')=',p[n,1]:4:5);
readkey;
DoneWinCrt;
end.
Для чтобы найти значение функции у(х) в точке х с помощью этой программы нужно сначала ввести количество узлов интерполирования, значение х, для которого надо найти значение функции, а потом ввести узлы интерполирования хi и соответствующие им значения функции уi и нажать клавишу ENTER.
Также для определения степени интерполирующего многочлена я написал программу konechn_razn.
Код программы:
program konechn_razn;
uses wincrt, strings;
var y:array [1..50,1..50] of real;
i,j,n: integer;
begin
StrCopy(WindowTitle, 'Программа построения конечных разностей ');
clrscr;
write ('vvedite k-vo znachenij funcii n=');
readln (n);
writeln ('vvedite mas Yi');
for i:=1 to n do
begin
write ('vvedite elem Y[',i,']=');
readln (y[i,1]);
end;
writeln ('PROCES VICHISLENIJA......');
for j:=2 to n do
begin
for i:=1 to n+1-j do
begin
y[i,j]:=y[i+1,j-1]-y[i,j-1];
end;
end;
writeln ('REZ MATRICA::::');
writeln (' Yi | Dyi ');
for i:=1 to n do
begin
for j:=1 to n+1-i do
begin
write (y[i,j]:4:5,' | ');
end;
writeln;
end;
readkey;
DoneWinCrt;
end.
Входными данными для этой программы есть: количество узлов интерполирования и значения функции yi, для которых надо построить конечные разности.
Решение в Excel
Для проверки вычислений я решил поставленную задачу в Excel по схеме Эйткина:
Также в целях проверки вычислений я решил данную задачу с помощью кубических сплайнов:
График, отображающий значения функции, вычисленные по схеме Эйткина и с помощью кубических сплайнов:
Выводы
Все многочлены, которые надо вичислить для данного х выражаются через определитель 2-го порядка, что делает вычисления единообразными. Схему Эйткина просто программировать.
Можно отметить то, что схема Эйткина применима и в случае неравноотстоящих узлов интерполирования, то есть ее можно применять для любого шага интерполирования. Также надо отметить то, что, если в задаче требуется вычислить значение функции в одной точке, нет необходимости строить общее выражение многочленна Лагранжа или Ньютона явно, а требуется только вичислить его значение в точке х. Эти вычисления удобно выполнить по интерполяционной схеме Эйткина.
Сопоставим исходные данные, у нас имеется 6 узлов интерполирования. По этим точкам можно построить интерполяционный полином, причем 5-й степени, привлекая к исследованию интерполяцию кубическим сплайном, утверждаю, что данным методом можно построить на каждом подинтервале полином 3-й степени. Последним словом в выборе между первым и вторым методом будут конечные разности на заданном множестве узлов. Конечные разности являются аналогом производной от функции. В данном случае конечные разности использованы для определения степени полинома и для определения полином данная функция или нет, с помощью которого можно максимально приблизить данную функцию.
Данного количества узлов интерполирования не достаточно для точного определения является ли данная функция полиномом, то есть в данном случае конечные разности не являются точным критерием для выбора между двумя методами интерполирования.
| Эйткин | |
| x= | 0,47 |
| y= | 0,45289 |
| сплайн | |
| x= | 0,47 |
| y= | 0,45277 |
В результате вычисления значения функции в точке 0,47 видно что значения функции в искомой точке мало отличимые. То есть в данном случае можно применять оба метода.
Если взять точность вычисления до четвертого знака после запятой, то степень полинома по данным конечных разностей будет полином 3-й степени. Поскольку по схеме Эйткина строятся все полиномы степени не выше 6-й. И в этом случае лучше применять кубические сплайны.
Список использованной литературы
1. Б. П. Демидович и И. А. Марон. “Основы вычислительной математики”, Москва, 1963г.
2. Н.С.Бахвалов, Н.П.Жидков, Г.М.Кобельков. “Численные методы”, Москва, 1987г.
3. Козин А. С., Лященко Н. Я. Вычислительная математика: Пособие для факультативных занятий в 10 классе.- К.: Рад. школа, 1983. – 191 с.
4. Мусіяка В. Г. Основи чисельних методів механіки: підручник. – К.: Вища освіта, 2004. – 240 с.: іл.
5. Л. Д. Назаренко Чисельні методи. Дистанційний курс.
Приложение
Результаты работы программы EITKIN:
Результаты работы программы konechn_razn:
