Министерствонауки и образования Украины
Сумскойгосударственный университет
кафедраинформатики
Численные методыКурсовая работа
на тему:“ Выбор интерполирующейфункции к заданной и ее построение ”Сумы 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 точек />/>,которые называются узлами интерполирования, и значения некоторой функции /> в этих точках: />. Нужно построить функцию /> ( функцию, котораяинтерполирует), которая совпадала бы с /> вузлах интерполяции и приближала ее между ними, то есть такую, что />. Геометрическаяинтерпретация задачи интерполяции состоит в том, что нужно найти такуюкривую/> некоторого вида,что проходит через заданную систему точек /> Спомощью этой кривой можно найти приближенное значение />, де /> Задача интерполяциистановится однозначной, если вместо произвольной функции /> искать многочлен /> степени не выше />, который удовлетворяетусловия:
/>.
Интерполяционныймногочлен /> всегда однозначный,поскольку существует только один многочлен степени />,который в данных точках принимает заданные значения. Существует несколькоспособов построения интерполяционного многочлена. Дальше мы рассмотрим основныеспособы подробнее.
Теоретическая часть
Интерполяционныймногочлен Лагранжа
Интерполяционныймногочлен Логранжа, что принимает в узлах интерполяции /> соответственно значений /> имеет вид:
/> (*)
С формулы видно, чтостепень многочлена /> равна />, и многочлен Логранжа удовлетворяетвсе условия задачи интерполяции.
Если расстояние между всемисоседними узлами интерполирования одинаково, то есть />, формула (*) значительно упрощается.Введем новую переменную />, тогда /> /> Теперь интерполяционныйполином Лагранжа имеет вид:
/>. (**)
Тут />.
Коэффициенты, которыестоят перед величинами /> в формуле (**),не зависят от функции /> и от шага />, а зависят только отвеличин /> Поэтому таблицамисоставленными для различных значений />, можновоспользоватся при решении различных задач интерполирования для равноотстоящихузлов.
Возникает вопрос, насколько близко многочлен Логранжа приближается к функции /> в других точках (неузловых), то есть на сколько большой остаток. На функцию /> накладывают дополнительныеограничения. А именно: предполагают, что в рассмотренной области /> изменения />, которые содержат узлыинтерполяции, функция /> имеет всепроизводные /> до />-го порядка включительно. Тогдаоценка абсолютной погрешности интерполяционной формулы Логранжа имеет вид:
/>, (***)
где /> .
Интерполяционныймногочлен Ньютона
Разделенными разностяминазываются соотношения вида:
— первого порядка:/>
/> />
— второго порядка:
/>/> (5.15)
…………………………………………………;
— n — го порядка:
/>
С помощью разделенныхрізностей можно построить многочлен:
/>/> (5.16)
Он называетсяинтерполяционным многочленНьютона для заданной функции. Эта форма записи более удобна для использования,поскольку при добавлении к узлам x0, 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 и требуется вычислить ее значение в точке х/>. В этом случае по схемеЭйткина в точке х вычисляют сначала значения двух линейных многочленов
/> и />,
а затем значениеквадратичного многочлена вида:
/>.
Непосредственнойподстановкой убеждаемся, что />,
/>; />,/>, />.
Покажем еще, что /> совпадает с формулой Лагранжадля трех узлов интерполирования. Поскольку
/>
/>,
то, раскрываяопределитель, получаем:
/>
Эта схема обобщается наболее высокие степени. Если функция 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]ofreal;
p:array [1..60,1..60] ofreal;
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 ndo
begin
write('vvedite elem X[',i,']=');
readln (x[i]);
end;
writeln('vvedite mas Yi');
for i:=1 to ndo
begin
write('vvedite elem Y[',i,']=');
readln (p[1,i]);
end;
writeln ('PROCESVICHISLENIJA......');
for i:=2 to ndo
begin
for j:=1 ton+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 ('REZMATRICA::::');
for i:=1 to ndo
begin
write('P^',i,'(',z:4:5,') | ');
for j:=1 ton+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 ndo
begin
write('vvedite elem Y[',i,']=');
readln (y[i,1]);
end;
writeln ('PROCESVICHISLENIJA......');
for j:=2 to ndo
begin
for i:=1 ton+1-j do
begin
y[i,j]:=y[i+1,j-1]-y[i,j-1];
end;
end;
writeln ('REZMATRICA::::');
writeln (' Yi |Dyi ');
for i:=1 to ndo
begin
for j:=1 ton+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:
/>