Интерполирование таблично заданных функций полиномами и сплайнами

Цель работы

Освоить интерполирование таблично заданных функций полиномами и сплайнами.

Задание

Дан массив A[0.30].

1. Рабочие точки интерполировать в пакете Mathcad с помощью канонического полинома, полинома Лагранжа, Ньютона. Оценить близость интерполяционной функции к контрольным точкам с помощью СКО:

где — разность между ординатами контрольной точки и интерполяционного полинома, n=15.

2. Рабочие точки интерполировать в пакете Mathcad с помощью линейного, квадратичного, кубического сплайнов. Оценить близость интерполяционного сплайна к контрольным точкам, рассчитав значение среднеквадратического отклонения (СКО). Выбрать наилучший с этой точки зрения сплайн.

3. Реализовать на языке высокого уровня интерполяцию линейным сплайном массива рабочих точек.

Теория

Интерполяция — операция приближения функции, заданной в отдельных точках внутри некоторого заданного промежутка.

Простейшая задача интерполяции заключается в следующем: на отрезке [a;b] заданы (n+1) точек которые называются узлами интерполяции, а также заданы значения некоторой функции f (x) в этих точках:

Требуется построить интерполирующую функцию F (x), принимающую в узлах интерполяции те же значения, что и f (x) т. е. :

Геометрически это означает, что нужно найти некоторую кривую y=F (x) определённого типа, проходящую через заданную систему точек, где .

В такой постановке задача интерполяции, может иметь бесчисленное множество решений, однако задача становится однозначной, если вместо произвольной функции F (x) искать полином P_n(x) степени n, удовлетворяющий условиям:

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

Наиболее часто используемые интерполяционные полиномы:

1. Канонический полином;

2. Полином Лагранжа;

3. Полином Ньютона.

Канонический полином

Будем искать интерполирующую функцию F (x) в виде канонического полинома степени n:

Выбор многочлена степени _n основан на том факте, что через (n+1) точку проходит единственная кривая степени n. Подставив в, получим СЛАУ:

Решая эту систему относительно переменных, найдём коэффициенты интерполяционного полинома.

Полином Лагранжа

Будем искать его в виде:

где) — функция, удовлетворяющая в узлах следующему свойству:

Таким образом, полином Лагранжа выражается следующей формулой:

Полином Ньютона

Ньютон предложил записать интерполирующую функцию в виде следующего полинома n-1 степени:

Нам известно, что. Подставим это значение в полином, вычислив тем самым значение коэффициента

Подставим в полином, для вычисления воспользуемся следующим соотношением:

отсюда

где — разделенная разность первого порядка, которая стремится к первой производной функции при. Аналогичным образом можно найти остальные разделённые разности первого порядка.

Все полученные результаты более наглядно можно представить в виде таблицы:

Основную роль в данной таблице играют диагональные элементы, которыми являются коэффициенты A₀, A₁,…, A_n. Остальные являются вспомогательными.

Все три рассмотренных полинома — это три различных способа построения одной и той же кривой:

Полином Лагранжа считается целесообразным использовать в тех случаях, когда задано небольшое количество точек. В остальных случаях рекомендуется использовать полином Ньютона.

Интерполяция сплайнами

Сплайн-функцияспециальным образом построенная гладкая кусочно-полиномиальная функция, средство интерполяции сложных кривых.

Пусть дан отрезок [a, b], который разбит точками на n частичных отрезков. Тогда сплайном степени m называется функция, обладающая следующими свойствами:

1) Функция, непрерывна на отрезке [a, b] вместе со всеми своими производными до некоторого порядка p:

2) На каждом частичном отрезке функция совпадает с некоторым полиномом

Разность (m-p) между степенью сплайна и наивысшим порядком производной называется дефектом сплайна.

Линейная интерполяция

В случае линейной интерполяции сплайном является линейная функция для отрезков .

Коэффициенты a и b в этом случае рассчитываем по формулам:

Найдя коэффициенты, можно рассчитать значения функции в любой точке интервала .

Решение СЛАУ в математическом пакете Mathcad

Задаем массивы:

Точки массива:

Выделим рабочие точки:

Сформируем матрицу, А путём создания СЛАУ и решением её методом LU разложения. Получившиеся ответы и будут значениями матрицы.

Получим интерполяционный полином:

Оценим близость интерполяционной функции к контрольным точкам с помощью СКО:

Посчитаем коэффициент r:

и подставим его:

Вычисляем коэффициенты с помощью линейного сплайна:

Вычислим значения в точках, выполнив интерполирование:

Снова оценим близость интерполяционной функции к контрольным точкам с помощью СКО:

Вычисляем коэффициенты с помощью квадратичного сплайна:

Оцениваем близость интерполяционной функции к контрольным точкам с помощью СКО:

Вычисляем коэффициенты с помощью кубического сплайна:

Оцениваем близость интерполяционной функции к контрольным точкам с помощью СКО:

— ближе всех интерполяция к контрольным точкам получилась при использовании метода кубических сплайнов.

Листинг и результат работы программы

интерполирование полином сплайн уравнение

Для решения системы линейных алгебраических уравнений я использовал алгоритмический язык Паскаль.

program Spline;

uses Crt, Graph;

{Структура, описывающая сплайн на каждом сегменте сетки}

type SplineTuple=record

a, b, c, d, x: real;

end;

const

n_max=100; {Максимальный размер массива}

var

Gd, Gm: integer; {Графические драйвер и мод}

PathToDriver: string;

splines: array[0.n_max-1] of SplineTuple; {Сплайн}

x, y: array[0.n_max-1] of real;

cur_x: real;

n:integer;

{Построение сплайна}

{x — узлы сетки, упорядоченные по возрастанию}

{y — значения функции в узлах сетки}

{n — количество узлов сетки}

Procedure BuildSpline;

var i: integer;

alpha:array[0.n_max-1] of real;

beta:array[0.n_max-1] of real;

h_i:real;

h_i1:real;

A:real;

C:real;

B:real;

F:real;

z :real;

begin

{Инициализация массива сплайнов}

for i:=0 to n-1 do

begin

splines[i]. x := x[i];

splines[i].a := y[i];

end;

splines[0].c:=0;

splines[n — 1]. c:=0;

{Решение системы линейных уравнений относительно коэффициентов сплайнов}

{Вычисление прогоночных коэффициентов — прямой ход метода}

alpha[0]: =0;

beta[0] :=0;

for i:=1 to n-2 do

begin

h_i := x[i] - x[i — 1];

h_i1:= x[i + 1] - x[i];

A:= h_i;

C:= 2.0 * (h_i + h_i1);

B:= h_i1;

F:= 6.0 * ((y[i + 1] - y[i]) / h_i1 — (y[i] - y[i — 1]) / h_i);

z:= (A * alpha[i — 1] + C);

alpha[i]: = -B / z;

beta[i]: = (F — A * beta[i — 1]) / z;

end;

{Нахождение решения — обратный ход метода прогонки}

for i:=n-2 downto 1 do

splines[i]. c := alpha[i] * splines[i + 1]. c + beta[i];

{По известным коэффициентам c[i] находим значения b[i] i d[i]}

for i:=n-1 downto 1 do

begin

h_i:= x[i] - x[i — 1];

splines[i].d:= (splines[i]. c — splines[i — 1]. c) / h_i;

splines[i]. b:= h_i * (2.0 * splines[i]. c + splines[i — 1]. c) / 6.0 + (y[i] - y[i — 1]) / h_i;

end;

end;

{Вычисление значения интерполяционной функции в произвольной точке}

function Func (x:real):real;

var

s:SplineTuple;

i, j, k:Integer;

dx:real;

begin

if (x <= splines[0]. x) then

{Если x меньше точки сетки x[0]}

s := splines[1]

else if (x >= splines[n — 1]. x) then

{Если x больше точки сетки x[n — 1] - пользуемся последним элементом массива}

s := splines[n — 1]

else {Иначе x лежит между граничными точками сетки, то производим бинарный поиск нужного элемента массива}

begin

i:= 0;

j:= n — 1;

while (i + 1 < j) do

begin

k:= i + (j — i) div 2;

if (x <= splines[k]. x) then

j:= k

else

i:= k;

end;

s := splines[j];

end;

dx:= (x — s. x);

{Вычисление значения сплайна в заданной точке по схеме Горнера}

Func:=s.a + (s.b + (s.c / 2.0 + s. d * dx / 6.0) * dx) * dx;

end;

begin

ClrScr;

Gd:=Detect;

Gm:=0;

PathToDriver:='';

InitGraph (gd, gm, PathToDriver); {Инициализируем графику}

if GraphResult<>grok then halt;

n:=16;

x[0]: =0;y[0]:=260;

x[1]:=20;y[1]:=260;

x[2]:=40;y[2]:=300;

x[3]:=60;y[3]:=270;

x[4]:=80;y[4]:=330;

x[5]:=100;y[5]:=230;

x[6]:=120;y[6]:=250;

x[7]:=140;y[7]:=220;

x[8]:=160;y[8]:=210;

x[9]:=180;y[9]:=170;

x[10]:=200;y[10]:=150;

x[11]:=220;y[11]:=230;

x[12]:=240;y[12]:=180;

x[13]:=260;y[13]:=190;

x[14]:=280;y[14]:=220;

x[15]:=300;y[15]:=150;

BuildSpline;

cur_x:=100;

MoveTo (Round (cur_x), Round (Func (cur_x)));

while (cur_x<300) do

begin

cur_x:=cur_x+0.1;

LineTo (Round (cur_x), Round (Func (cur_x)));

MoveTo (Round (cur_x), Round (Func (cur_x)));

end;

ReadKey;

CloseGraph;

END.

Результат работы программы

Вывод

В результате выполнения данной лабораторной работы были исследованы различные методы интерполяции функции, которые были заданы табличным путем. Согласно алгоритму, изложенному в книге Мудрова А. Е. «Численные методы для ПВЭМ на языках высокого уровня», написана программа на языке высокого уровня Pascal. Реализовано решение в математическом пакете Mathcad. Наиболее точным был результат, который получен при выполнении интерполяции линейным сплайном, по причине наименьшей величины СКО.

1. Мудров А. Е., «Численные методы для ПЭВМ на языках Бейсик, Фортран и Паскаль». — Томск: МП «РАСКО», 1991.-272 с.: ил.