«Интерполяция функций с помощью сплайна» - Лабораторная работа

Содержание

1 Цель работы:

Целью работы является ознакомление студентов с задачей интерполяции функций, с методом прогонки для решения систем линейных алгебраических уравнений с ленточной матрицей, с понятием сплайна, получение навыков решения задач вычислительной математики на ЭВМ.

Введение

2 Описание метода:

Пусть отрезок [a,b] разбит на n равных частей [xi, xi+1], где xi=a+ih, i=0,.,n, xn=b, h=(b-a)/n.

Сплайном называется функция, которая вместе с несколькими производными непрерывна на всем заданном отрезке [a,b], а на каждом частичном отрезке [xi, xi+1] в отдельности является некоторым алгебраическим многочленом.

Максимальная по всем частичным отрезкам степень многочленов называется степенью сплайна, а разность между степенью сплайна и порядком наивысшей непрерывной на [a,b] производной - дефектом сплайна.

Определение. Функция Sn,(x) называется сплайном степени n дефекта  (-целое число, 0n+1) с узлами на сетке  (: a=x0< а) на каждом отрезке [xi,x i+1] функция Sn, (x) является многочленом степени n, то есть

для x[xi, xi+1] , i=0,.,n-1;

б) S n,(x)Cn-v[a,b]

(для целого k>0 через Ck =Ck[a,b] обозначается множество k раз непрерывно дифференцируемых на [a,b] функций).

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

На практике широкое распространение получили сплайны третьей степени, имеющие на [a,b] непрерывную, по крайней мере, первую производную. Эти сплайны называются кубическими и обозначаются S3(x) (без указания дефекта).

Пусть на отрезке [a,b] в узлах сетки  заданы значения некоторой функции

fi =f(xi), i=0,.,n.

Интерполяционным кубическим сплайном S3(x) называется сплайн

S3(x)=аi0 +аi1(x - xi)+аi2(x - xi)2 +аi3(x - xi)3, x[xi, xi+1], (1.1)

удовлетворяющий условиям

S3(xi)=f(xi), i=0,.,n. (1.2)

Сплайн (1.1) на каждом из отрезков [xi, xi+1], i=0,.,n-1 определяется четырьмя коэффициентами, и поэтому для его построения на всем промежутке [a,b] необходимо определить 4n коэффициентов. Для их однозначного определения необходимо задать 4n уравнений.

Условие (1.2) дает 2n уравнений, при этом функция S3(xi), удовлетворяющая этим условиям, будет непрерывна во всех внутренних узлах.

Условие непрерывности производных сплайна , r=1,2 во всех внутренних узлах xi, i=1,.,n-1 сетки  дает 2(n-1) равенств.

Вместе получается 4N-2 уравнений.

Два дополнительных условия обычно задаются в виде ограничений на значение производных сплайна на концах промежутка [a,b] и называются краевыми условиями.

Наиболее употребительны следующие типы краевых условий:

а) S\'3(а)=f\'(а), S\'(b)=f\'(b) ;

б) S\"3(а)=f\"(а), S\"(b)=f\"(b) ;

в) ;

г) S\'\'\'3(xp+0)=S\'\'\'3(xp-0), р =1, n-1.

Через краевые условия в конструкцию сплайна включаются параметры, выбирая которые можно управлять его поведением, особенно возле концов отрезка [a,b].

Условия типа в) носят названия периодических. Естественно требовать их выполнения в том случае, когда интерполируемая функция периодическая с периодом (b-a).

Если известны f\'(x) или f\"(x) в точках а и b, то естественно воспользоваться краевыми условиями типа а) или б).

Если производные неизвестны, то в большинстве случаев наилучшим решением будет применение краевых условий типа г).

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

Иногда предлагается принимать

S\"3(а)=S\"3(b)=0 .

В этом случае вблизи концов точность интерполяции функции и ее первой производной уменьшается и становится соизмеримой с точностью интерполяции сплайном первой степени, что резко ухудшает всю картину.

Выдержка из текста работы

Листинг программы:

#include

#include

#include

#include

float d[99999];

float l[99999];

float m[99999];

double f[99999];

float s[99999],mu[99999];

void main()

{

system(\"cls\");

int n;

float c,a,b,shag;

float h;

cout<<\"Vvedite shag h:=\";

cin>>h;

cout<<\"Vvedite kolichestvo znachenii n:=\";

cin>>n;

shag=h/2;

c=h/6;

a=2*h/3;

b=h/6;

cout<<\"c=\"< int i;

float x[99999],y[99999];

for(i=1;i {

x[i]=(float)i;

}

for(i=1;i {

cout<<\"Vvedite Y[\"< cin>>y[i];

}

mu[0]=m[0]=m[n]=l[0]=0;

for(i=2;i {

d[i]=((y[1+i]-y[i])/h)-(y[i]-y[i-1])/h;

}

for(i=1;i {

l[i]=-b/(a+c*l[i-1]);

mu[i]=(d[i]-c*mu[i-1])/(a+c*l[i-1]);

}

for(i=n-1;i>1;i--)

{

m[i]=l[i]*m[i+1]+mu[i];

}

double po;

for(i=2;i {

f[i]=x[i]-shag;

Заключение

Проверим также метод интерполяции на более монотонной функции, например на функции y(x)=x2. В качестве узловых точек возьмем точки из отрезка [1;8], с шагом h=1:

Vvedite shag h:=1

Vvedite kolichestvo znachenii n:=8

c=0.166667 a=0.666667 b=0.166667

Vvedite Y[1]:=1

Vvedite Y[2]:=4

Vvedite Y[3]:=9

Vvedite Y[4]:=16

Vvedite Y[5]:=25

Vvedite Y[6]:=36

Vvedite Y[7]:=49

Vvedite Y[8]:=64

S[1.5]=2.35101

S[2]= 4

S[2.5]=6.18125

S[3]= 9

S[3.5]=12.2424

S[4]= 16

S[4.5]=20.2242

S[5]= 25

S[5.5]=30.2359

S[6]= 36

S[6.5]=42.2074

S[7]= 49

Для продолжения нажмите любую клавишу . . .

Вывод: Как мы можем заметить, при интерполировании более монотонной функции, погрешность расчетов оказалось относительно небольшой.