Реферат: Выбор и построение интерполирующей функции
Министерство науки и образования Украины
Сумской государственный университет
кафедра информатики
Численные методы
Курсовая работа
на тему:
“ Выбор интерполирующей функции к заданной и ее построение ”
Сумы 2006
Содержание
Постановка задачи.
1. Введение.
2. Теоретическая часть.
3. Практическая реализация:
3.1 Программа на языке Pascal.
3.2 Решение в Excel.
4. Выводы.
Список использованной литературы.
Приложение.
Постановка задачи
Найти значение функции у в точке х=0.47, используя интерполяционную схему Эйткина, проверить правильность решения с помощью кубического сплайна. Значения функции у приведены в таблице:
i
1
2
3
4
5
xi
0,4
0,5
0,6
0,7
0,8
0,9
yi
0,38942
0,47943
0,56464
0,64422
0,71736
0,78333
x=
0,47
Введение
Пусть на отрезке />задано Nточек/>/>, которые называются узлами интерполирования, и значения некоторой функции />в этих точках:/>. Нужно построить функцию />( функцию, которая интерполирует), которая совпадала бы с />в узлах интерполяции и приближала ее между ними, то есть такую, что />. Геометрическая интерпретация задачи интерполяции состоит в том,что нужно найти такую кривую/>некоторого вида, что проходит через заданную систему точек />С помощью этой кривой можно найти приближенное значение />, де />Задача интерполяции становится однозначной, если вместо произвольной функции />искать многочлен />степени не выше />, который удовлетворяет условия:
/>.
Интерполяционный многочлен />всегда однозначный, поскольку существует только один многочлен степени />, который в данных точках принимает заданные значения. Существует несколько способов построения интерполяционного многочлена. Дальше мы рассмотрим основные способы подробнее.
Теоретическая часть
Интерполяционный многочлен Лагранжа
Интерполяционный многочлен Логранжа, что принимает в узлах интерполяции />соответственно значений />имеет вид:
/>(*)
С формулы видно, что степень многочлена />равна />, и многочлен Логранжа удовлетворяет все условия задачи интерполяции.
Если расстояние между всеми соседними узлами интерполирования одинаково, то есть />, формула (*) значительно упрощается. Введем новую переменную />, тогда />/> Теперь интерполяционный полином Лагранжа имеет вид:
/>. (**)
Тут />.
Коэффициенты, которые стоят перед величинами /> в формуле (**), не зависят от функции /> и от шага />, а зависят только от величин /> Поэтому таблицами составленными для различных значений />, можно воспользоватся при решении различных задач интерполирования для равноотстоящих узлов.
Возникает вопрос, на сколько близко многочлен Логранжа приближается к функции />в других точках (не узловых), то есть на сколько большой остаток. На функцию />накладывают дополнительные ограничения. А именно: предполагают, что в рассмотренной области />изменения />, которые содержат узлы интерполяции, функция />имеет все производные />до />-го порядка включительно. Тогда оценка абсолютной погрешности интерполяционной формулы Логранжа имеет вид:
--PAGE_BREAK--/>, (***)
где />.
Интерполяционный многочлен Ньютона
Разделенными разностями называются соотношения вида:
/>— первого порядка:
/>/>
— второго порядка:
/>/>(5.15)
…………………………………………………;
— n — го порядка:
/>
С помощью разделенных різностей можно построить многочлен:
/>/>(5.16)
Он называется интерполяционным многочленНьютона для заданной функции. Эта форма записи более удобна для использования, поскольку при добавлении к узлам x, x1, …, xnнового xn+1все вычесленные раньше члены остаются без изменений, а в формулу добавляется только одно слогаемое. При использовани формулы Логранжа нужно вычислять все заново.
Если значения функции заданы для равноотстоящих значений аргумента />(постоянную величину />, i=0,1,…,nназывают шагом интерполяции), то интерполяционный многочленпринимает вид:
/>(5.17)
Здесь />— конечные разности к-го порядка. Они определяются по формуле/>где />-биномиальные коэффициенты.
Сравнивая эту формулу с предыдущей, легко установить, что при />/>конечные и разделенные разности связаны соотношением вида:
/>(5.18)
Для практического использования формулу (5.17) записывают в преобразованном виде. Для этого введем новую переменную />, положив /> где /> — количество шагов />, необходимое для достижения точки /> из точки />. Таким образом получим первую интерполяционную формулу Ньютона для интерполирования вперед, то есть в начале таблицы значений:
/>(5.19)
Предположим, что точка интерполяции расположена вблизи конечной точки /> таблицы. В этом случае узлы интерполяции следует брать таким образом /> Формула Ньютона для интерполяциии назад имеет вид:
/>(5.20)
Разделенные разности можно выразить через конечные разности, если воспользоваться возможностью переставлять в них аргументы, и соотношением (5.18), откуда следует:
/>; />
/>
Введем переменную />, учитывая что /> получим для /> вторую интерполяционную формулу Ньютона для интерполяции в конце таблицы :
/>.
Как первая, так и вторая интерполяционные формулы Ньютона могут быть использованы для екстраполяции функции, то есть для вычисления значений функции />, значения аргументов /> которой лежат вне таблицы. Если /> и значение /> близко к />, то выгодно использовать первый интерполяционный многочлен Ньютона, тогда /> и /> Таким образом, первая интерполяционная формула Ньютона применяется для интерполяции вперед и екстраполяции назад, а вторая — наоборот, для интерполяции назад и екстраполяции вперед.
Отметим, что операция екстраполирования, вообще говоря, менее точная чем операция интерполяции.
Интерполяционные формулы Ньютона выгодны, поскольку при добавлении /> новых узлов интерполяции необходимые дополнительные вычисления только для /> новых членов, без изменения старых.
Схема Эйткина
Пусть дана f задана таблично в точках хi она принимает значения уi= f(хi) (i=0,1,…,n). Требуется вычислить значение функции f в некоторой точке х/>, не совпадающей с точками хi. В таком случае нет необходимости строить общее выражение многочленна Лагранжа явно, а требуется только вичислить его значение в точке х. Эти вычисления удобно выполнить по интерполяционной схеме Эйткина. Характерной чертой этой схемы является единообразие вичислений.
Если функция f задана в двух точках х0и х1 значениями у0и у1, то для вычисления ее значения в точке х/> можно воспользоваться формулой:
/>(*) линейного интерполирования.
Обозначив значение функции в точке x через />, формулу (*) можно представить в таком виде:
/>,
Где в правой части стоит определитель 2-го порядка. Эта формула эквивалентна формуле (*). Кроме того, />, />.
Пусть функция f задана в трех точках х0, х1 и х2 своими значениями у0, у1 и у2 и требуется вычислить ее значение в точке х/>. В этом случае по схеме Эйткина в точке х вычисляют сначала значения двух линейных многочленов
/>и />,
а затем значение квадратичного многочлена вида:
/>.
Непосредственной подстановкой убеждаемся, что />,
/>; />, />, />.
продолжение--PAGE_BREAK--
Покажем еще, что /> совпадает с формулой Лагранжа для трех узлов интерполирования. Поскольку
/>
/>,
то, раскрывая определитель, получаем:
/>
Эта схема обобщается на более высокие степени. Если функция f задана в четырех точках, то кубическое интерполирование выполняется по формуле
/>,
Где /> и /> — значения квадратичных многочленов в точке х/>. Непосредственной проверкой убеждаемся, что /> и />. Кроме того /> совпадает с кубическим интерполяционным многочленом Лагранжа:
/>.
Вообще, если в (n+1)-й точке хi (i=0,1,…,n) функция f принимает значения yi (i=0,1,…,n), то значение интерполяционного многочлена Лагранжа степени n в точке х/> можно вычислить по формуле
/>,
где /> и /> — значения интерполяционных многочленов, вычисленных в точке х на предшествующем шаге. Ясно, что для вычисления значения многочлена степени n в точке х необходимо по схеме Эйткина вычислить в этой точке значения n линейных, n-1 квадратичных, n-2 кубических многочленов и т. д., два многочлена степени n-1 и, наконец, один многочлен степени n. Все эти многочлены выражаются через определитель 2-го порядка, что делает вычисления единообразными.
Отметим то, что схема Эйткина применима и в случае неравноотстоящих узлов интерполирования.
Сплайн – интерполяция
В инженерной практике график функции y(xi) (i=0,N) строят в основном с помощью лекал. Если точки размещены редко, то пользуются гибкой линейкой (spline), ставят ее на ребро и изгибают так, чтобы она одновременно проходила через все точки.
Поскольку приближенное уравнение изгиба пружинистого бруса имеет вид />, то можно допустить, что ее форма между узлами есть алгебраический полином 3-й степени.
Вероятно, интерполирующую функцию между каждыми двумя узлами можно взять, например, в таком виде:
/>(*)
/>.
Неизвестные коэффициенты ai, bi, ci, di найдем с условий в узлах интерполяции.
Поскольку полиномы совпадают с табличными значениями функции y(xi) (i=1,N) в узлах интерполяции, то:
/>(А)
/>(В)
Поскольку этих уравнений в два раза меньше, чем неизвестных коэффициентов, то надо еще какие-нибудь дополнительные условия (например, условия непрерывности 1-й и 2-й производных во всех точках, в том числе и в узлах интерполирования, то есть условия гладкости угла поворота пересечения и кривизны линейки).
С условий непрерывности производных у внутренних узлах /> имеем:
/>(С)
/>(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;
продолжение--PAGE_BREAK--
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-й. И в этом случае лучше применять кубические сплайны.
Список использованной литературы
Б. П. Демидович и И. А. Марон. “Основы вычислительной математики”, Москва, 1963г.
Н.С.Бахвалов, Н.П.Жидков, Г.М.Кобельков. “Численные методы”, Москва, 1987г.
Козин А. С., Лященко Н. Я. Вычислительная математика: Пособие для факультативных занятий в 10 классе.- К.: Рад. школа, 1983. – 191 с.
Мусіяка В. Г. Основи чисельних методів механіки: підручник. – К.: Вища освіта, 2004. – 240 с.: іл.
Л. Д. Назаренко Чисельні методи. Дистанційний курс.
Приложение
Результаты работы программы EITKIN:
/>
Результаты работы программы konechn_razn:
/>