Алгоритм для интерполяции функции естественным интерполяционным кубическим сплайном

РАСЧЕТНО-ГРАФИЧЕСКОЕ ЗАДАНИЕ 1

Задание 1. Построить алгоритм для интерполяции функции  на отрезке  интерполяционным линейным сплайном  с равномерным шагом. Разработать программу, которая реализует этот алгоритм и выводит на печать значение сплайна  и значение функции , значение производной сплайна и производной функции при . Размерность интерполяционной таблицы равна 21.

Задание 2. Построить алгоритм для интерполяции функции  на отрезке  естественным интерполяционным кубическим сплайном  с равномерным шагом. Разработать программу, которая реализует этот алгоритм и выдает на печать значение сплайна  и функции , значение производной сплайна и производной функции при . Размерность интерполяционной таблицы равна 21.

Варианты заданий

Номер варианта
1		-1	1	0.15
2		-2	2	0.1
3		-1.5	1.5	0.2
4		-1.2	1.2	0.3
5		-1.6	1.6	0.4
6		-1	1	0.25
7		-2	2	-0.1
8		-1.5	1.5	-0.2
9		-1.2	1.2	-0.3
10		-1.6	1.6	-0.4
11		-1	1	-0.15
12		-2	2	-0.5
13		-1.5	1.5	0.1
14		-1.2	1.2	0.2
15		-1.6	1.6	0.6
16		-1	1	-0.25
17		-2	2	0.3
18		-1.5	1.5	-0.1
19		-1.2	1.2	-0.2
20		-1.6	1.6	-0.6
21		-1	1	0.35
22		-2	2	-0.3
23		-1.5	1.5	0.4
24		-1.2	1.2	0.1
25		-1.6	1.6	0.2
26		-1	1	-0.35
27		-2	2	0.5
28		-1.5	1.5	-0.4
29		-1.2	1.2	-0.1
30		-1.6	1.6	-0.2

ПРИМЕР ВЫПОЛНЕНИЯ ЗАДАНИЯ 1

Построим алгоритм для интерполяции функции  на отрезке  интерполяционным линейным сплайном с равномерным шагом. Будем опираться на материал, изложенный в п.п. 2.1.

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

Узлы сетки  находятся по формуле ,  . Значения  находятся по формуле ,  . Таким образом, по формуле функции  на отрезке  построена интерполяционная таблица  размерности 21 с равномерным шагом. Так как все узлы сетки различны, то существует единственный интерполяционный линейный сплайн , удовлетворяющий этой таблице. Для построения сплайна  сначала найдем его коэффициенты по формулам

,   ,   ,  ,   .

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

.

Программа

#include <stdio.h>

#include <conio.h>

#include <math.h>

float A, B;          // границы отрезка [а, b]

float Xi[21], Yi[21]; // интерполяционная таблица

float Ai[21], Bi[21]; // коэффициенты сплайна ai, bi

float h;             // шаг сетки

int   n;             // число отрезков

// Функция sin(x)

float fn (float x)  { return (sin (x));}

// Вычисление значений сплайна L(x)

float L (float x)

{

int i=int ((x-A)/h);   // номер интервала

return Ai[i]+Bi[i]*(x-Xi[i]);

}

void main ()

{

// Инициализация начальных значений

A=-2, B=2;          // отрезок [-2,2]

n=20;               // число отрезков n

h=float (B-A)/n;    // шаг h=(b-a)/n

// Создание интерполяционной таблицы

for (int i=0; i<=n; i++)

{

Xi[i]=A+i*h;

Yi[i]=fn (A+i*h);

};

// Расчет коэффициентов сплайна L(x)

for (i=0; i<=n-1; i++)

{

Ai[i]=Yi[i];

Bi[i]=(Yi[i+1]-Yi[i])/h;

};

Ai[n]=Yi[n],

Bi[n]=Bi[n-1];

printf ("Значение функции f в точке x=%.2g равно %.2g\n", 0.1, fn (0.1));

printf ("Значение сплайна L в точке x=%.2g равно %.2g\n", 0.1, L (0.1));

getch ();

}

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

Значение функции f в точке x=0.1 равно 0.1

Значение сплайна L в точке x=0.1 равно 0.099

ПРИМЕР ВЫПОЛНЕНИЯ ЗАДАНИЯ 2

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

Как и в задании 1, число отрезков , шаг, где , . Введем сетку на отрезке : .

Узлы сетки  находятся по формуле ,  . Значения  находятся по формуле ,  . Таким образом, по формуле функции  на отрезке  построена интерполяционная таблица  размерности 21 с равномерным шагом. Так как все узлы сетки различны, то существует единственный интерполяционный кубический сплайн , удовлетворяющий этой таблице.

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

  ,

Эта система линейных уравнений имеет единственное решение, которое находится методом прогонки. Запишем формулы метода прогонки для данного примера. Сначала в порядке возрастания находятся прогоночные коэффициенты  и :

,   ,   ,   ,   , где   .

Затем в порядке убывания находятся неизвестные :

,   ,  

После нахождения коэффициентов  и , находим по явным формулам коэффициенты  и :

,    ,    ,

,   

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

Значение производной сплайна при  находим по формуле

.

Программа

#include <stdio.h>

#include <conio.h>

#include <math.h>

float A,B;                     // границы отрезка

float Xi[21],Yi[21];           // интерполяционная таблица

float Ai[21],Bi[21],Ci[21],Di[21];   // коэффициенты ai, bi, ci, di

float Alpha[21],Betta[21];     // прогоночные коэффициенты

float h;                       // шаг сетки

float n;                       // число отрезков

// Функция sin(x)

float fn (float x) { return (sin(x)); }

// Производная функции sin(x)

float fsh (float x) { return (cos(x)); }

// Вычисление значений сплайна S(x)

float Sx (float x)

{

int i=int((x-A)/h);

return Ai[i]+Bi[i]*(x-Xi[i])+Ci[i]*(x-Xi[i])*(x-Xi[i])+Di[i]*(x-Xi[i])*

(x-Xi[i])*(x-Xi[i]);

}

// Вычисление значений производной сплайна

float Sshx (float x)

{

int i=int((x-A)/h);

return Bi[i]+2*Ci[i]*(x-Xi[i])+3*Di[i]*(x-Xi[i])*(x-Xi[i]);

}

void main()

{

clrscr();

float Ti=0;       // вспомогательная  переменная ti

// Инициализация начальных значений

A=-2,B=2;         // отрезок [-2,2]

n=20;             // число отрезков n

h=float(B-A)/n;   // шаг

// Создание интерполяционной таблицы

for(int i=0;i<=n;i++)

{

Xi[i] = A+i*h;

Yi[i] = fn(Xi[i]);

}

// Находим коэффициенты сплайна ai, bi, ci, di

// Расчет коэффициентов ai

for(i=0;i<=n;i++) Ai[i]=Yi[i];

// Решение СЛУ с трехдиагональной матрицей методом прогонки

//                           для нахождения коэффициентов ci

// Вычисление прогоночных коэффициентов (прямая прогонка)

Alpha[1] = 0;

Betta[1]  = 0;

for(i=1;i<=n-1;i++)

{

Alpha[i+1] = (-1)/(4+Alpha[i]);

Ti         = 3*(Yi[i+1]-2*Yi[i]+Yi[i-1])/(h*h) ;

Betta[i+1] = (Ti-Betta[i])/(4+Alpha[i]);

}

// Вычисление коэффициентов ci

// Обратная прогонка

Ci[n]=0;

for(i=n-1;i>=0;i--) { Ci[i] = Alpha[i+1]*Ci[i+1]+Betta[i+1]; }

// Находим коэффициенты bi,di

for(i=0;i<=n-1;i++)

{

Bi[i] = ((Yi[i+1]-Yi[i])/h)-(h*(Ci[i+1]+2*Ci[i])/3);

Di[i] = (Ci[i+1]-Ci[i])/(3*h);

}

Bi[n] = Bi[n-1]+2*Ci[n-1]*h+3*Di[n-1]*h*h;

Di[n] = Di[n-1];

printf("Значение функции f в точке x=%.2g равно %.4g\n",0.1,fn(0.1));

printf("Значение сплайна S в точке x=%.2g равно %.4g\n",0.1,Sx(0.1));

printf("Значение производной функции f' в точке x=%.2g равно

%.4g\n",0.1,fsh(0.1));

printf("Значение производной сплайна S' в точке x=%.2g равно

%.4g\n",0.1,Sshx(0.1));

getch();

}

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

Значение функции f в точке x=0.1 равно 0.09983

Значение сплайна S в точке x=0.1 равно 0.09983

Значение производной функции f' в точке x=0.1 равно 0.995

Значение производной сплайна S' в точке x=0.1 равно 0.995

РАСЧЕТНО-ГРАФИЧЕСКОЕ ЗАДАНИЕ 2

Задание 1. Построить алгоритм для аппроксимации функции  на отрезке  регрессионной функцией  методом наименьших квадратов. Разработать программу, которая реализует этот алгоритм и выводит на печать значение среднеквадратичной погрешности и таблицу значений функций  и  на отрезке  в узлах сетки. Размерность таблицы равна 21, шаг равномерный.

Замечание. Для каждого варианта функция  выбирается та же, что и в расчетно-графическом задании 1.

Варианты заданий

Номер варианта	g(x)	a	b
1		-0.2	0.2
2		0	0.5
3		-0.2	0.2
4		-0.18	0.22
5		-0.2	0.2
6		-0.2	0.2
7		0	0.6
8		-0.4	0
9		0	0.5
10		-0.5	0.5
11		-0.2	0.2
12		0.5	1.5
13		-0.5	0.5
14		-0.5	0.5
15		-0.2	0.2
16		-0.5	0.5
17		0.5	1.5
18		0.5	0.7
19		-1	0.2
20		0.73	0.89
21		-0.2	0.2
22		-0.2	0.2
23		-0.2	0.2
24		0	0.2
25		-0.5	0.5
26		0.4	0.9
27		0.6	0.9
28		-0.5	0.5
29		0.5	1
30		0.4	0.8

Задание 2. Задано  уровней некоторого временного ряда , , . Построить алгоритм для нахождения прогноза с  по  уровень временного ряда методом наименьших квадратов для параболической зависимости. Разработать программу, которая реализует этот алгоритм и выводит на печать значения прогноза  и среднеквадратичной погрешности.

Варианты заданий

Номер варианта
1	20	-1.35; -1.66; -1.96; -2.26; -2.56; -2.86; -3.16; -3.45; -3.75; -4.04; -4.33; -4.61;-4.90; -5.18; -5.45; -5.73; -6.00; -6.26; -6.52; -6.78
2	21	0.75; 0.44; 0.14; -0.16; -0.46; -0.76; -1.06; -1.35; -1.65; -1.94; -2.23; -2.51; -2.80; -3.08; -3.35; -3.63; -3.90; -4.16; -4.42; -4.68; -4.93
3	17	1.81; 2.56; 3.31; 4.06; 4.80; 5.53; 6.25; 6.95; 7.64; 8.32; 8.98; 9.62; 10.24; 10.84; 11.41; 11.96; 12.48
4	19	1.07; 1.14; 1.21; 1.28; 1.35; 1.42; 1.48; 1.55; 1.61; 1.67; 1.72; 1.77; 1.82; 1.87; 1.91; 1.95; 1.98; 2.01; 2.04
5	18	0.69; 0.32; -0.04; -0.40; -0.76; -1.11; -1.46; -1.80; -2.13; -2.46; -2.77; -3.08;-3.38; -3.67; -3.95; -4.21; -4.46; -4.70
6	20	1.36; 1.42; 1.48; 1.53; 1.59; 1.65; 1.71; 1.77; 1.82; 1.88; 1.94; 2.00; 2.05; 2.11; 2.17; 2.22; 2.28; 2.33; 2.38; 2.44
7	21	0.93; 0.65; 0.38; 0.11; -0.16; -0.43; -0.69; -0.94; -1.19; -1.44; -1.67; -1.91; -2.13; -2.34; -2.55; -2.75; -2.93; -3.11; -3.27; -3.43; -3.57
8	19	1.56; 2.07; 2.58; 3.09; 3.58; 4.08; 4.56; 5.04; 5.51; 5.97; 6.41; 6.85; 7.27; 7.67; 8.06; 8.43; 8.78; 9.12; 9.43
9	18	1.78; 2.05; 2.33; 2.59; 2.85; 3.09; 3.33; 3.54; 3.75; 3.93; 4.09; 4.23; 4.35; 4.45; 4.52; 4.57; 4.60; 4.60
10	17	-2.07; -2.63; -3.19; -3.74; -4.28; -4.80; -5.31; -5.80; -6.27; -6.72; -7.14; -7.53; -7.89; -8.23; -8.53; -8.79; -9.02
11	20	2.21; 2.91; 3.61; 4.29; 4.96; 5.62; 6.25; 6.86; 7.45; 8.01; 8.53; 9.02; 9.47; 9.89; 10.26; 10.59; 10.88; 11.12; 11.31; 11.45
12	21	1.42; 1.44; 1.46; 1.48; 1.50; 1.51; 1.53; 1.55; 1.57; 1.59; 1.61; 1.63; 1.64; 1.66; 1.68; 1.70; 1.72; 1.73; 1.75; 1.77; 1.78
13	19	-1.63; -1.66; -1.69; -1.72; -1.75; -1.78; -1.81; -1.83; -1.86; -1.89; -1.92; -1.95; -1.98; -2.01; -2.03; -2.06; -2.09; -2.12; -2.14
14	18	0.11; -0.52; -1.16; -1.78; -2.40; -3.02; -3.62; -4.22; -4.80; -5.36; -5.91; -6.45; -6.97; -7.46; -7.94; -8.40; -8.83; -9.24
15	17	2.39; 3.02; 3.66; 4.28; 4.90; 5.52; 6.12; 6.72; 7.30; 7.86; 8.41; 8.95; 9.47; 9.96; 10.44; 10.90; 11.33
16	20	1.27; 1.22; 1.17; 1.13; 1.08; 1.03; 0.98; 0.93; 0.88; 0.84; 0.79; 0.74; 0.69; 0.65; 0.60; 0.55; 0.51; 0.46; 0.42; 0.37
17	21	0.68; 0.11; -0.46; -1.03; -1.60; -2.16; -2.72; -3.27; -3.82; -4.36; -4.89; -5.42; -5.94; -6.44; -6.94; -7.43; -7.91; -8.38; -8.84; -9.28; -9.71
18	19	-1.82; -2.39; -2.96; -3.53; -4.10; -4.66; -5.22; -5.77; -6.32; -6.86; -7.39; -7.92; -8.44; -8.94; -9.44; -9.93; -10.41; -10.88; -11.34
19	18	2.76; 4.14; 5.51; 6.87; 8.21; 9.55; 10.86; 12.15; 13.42; 14.66; 15.86; 17.04; 18.18; 19.28; 20.34; 21.35; 22.32; 23.24
20	17	-1.57; -1.90; -2.23; -2.55; -2.86; -3.16; -3.44; -3.70; -3.95; -4.17; -4.37; -4.54; -4.68; -4.80; -4.89; -4.95; -4.98
21	20	-2.11; -2.88; -3.65; -4.41; -5.15; -5.88; -6.59; -7.28; -7.95; -8.59; -9.19; -9.77; -10.31;-10.82; -11.28; -11.71; -12.09; -12.43; -12.72; -12.97
22	21	-1.82; -2.61; -3.40; -4.17; -4.94; -5.70; -6.44; -7.16; -7.86; -8.55; -9.21; -9.84; -10.45; -11.02; -11.57; -12.08; -12.56; -13.00; -13.40; -13.77; -14.09
23	19	0.74; 0.66; 0.58; 0.50; 0.42; 0.34; 0.26; 0.18; 0.10; 0.02; -0.06; -0.13; -0.21;-0.29; -0.36; -0.44; -0.51; -0.58; -0.65
24	18	0.43; 0.11; -0.21; -0.52; -0.83; -1.12; -1.40; -1.67; -1.93; -2.17; -2.39; -2.59; -2.77; -2.93; -3.07; -3.19; -3.28; -3.34
25	17	0.84; 0.78; 0.71; 0.65; 0.58; 0.52; 0.45; 0.39; 0.33; 0.26; 0.20; 0.14; 0.07; 0.01; -0.05; -0.11; -0.17
26	20	-0.59; -0.53; -0.48; -0.42; -0.36; -0.30; -0.25; -0.19; -0.13; -0.08; -0.02; 0.03; 0.09; 0.14; 0.20; 0.25; 0.31; 0.36; 0.41; 0.46
27	21	-0.72; -0.78; -0.85; -0.91; -0.98; -1.05; -1.11; -1.18; -1.24; -1.31; -1.37; -1.43; -1.50; -1.56; -1.62; -1.69; -1.75; -1.81; -1.87; -1.93; -1.99
28	19	-0.58; -0.51; -0.44; -0.37; -0.31; -0.24; -0.17; -0.10; -0.03; 0.03; 0.10; 0.17; 0.23; 0.30; 0.36; 0.42; 0.49; 0.55; 0.61
29	18	-0.40; 0.05; 0.49; 0.94; 1.38; 1.82; 2.26; 2.70; 3.14; 3.57; 4.00; 4.42; 4.84; 5.25; 5.66; 6.06; 6.46; 6.85
30	17	-1.30; -1.75; -2.19; -2.64; -3.08; -3.52; -3.96; -4.40; -4.84; -5.27; -5.70; -6.12; -6.54; -6.95; -7.36; -7.76; -8.16

ПРИМЕР ВЫПОЛНЕНИЯ ЗАДАНИЯ 1

Задание. Построить алгоритм для аппроксимации функции  на отрезке  регрессионной функцией  с использованием метода наименьших квадратов. Разработать программу, реализующую этот алгоритм. Размерность равномерной сетки равна 21.

Решение. Построение алгоритма в данном случае заключается в построении системы линейных уравнений и обоснования существования и единственности решения этой системы линейных уравнений. В дальнейшем будем опираться на материал, изложенный в п.п. 4.3.

Сначала построим сетку и таблицу данных. Так как сетка равномерная и число узлов сетки равно 21, то число отрезков  равно 20 (). Обозначим через  шаг, который вычисляется по формуле . В нашем случае , .

Введем сетку на отрезке : . Узлы сетки находятся по формуле ,  . Значение находится по формуле ,  . Таким образом, по формуле функции  на отрезке  построена таблица данных