Метод наименьших квадратов. Простая линейная регрессия, страница 8

Построим алгоритм для интерполяции функции  на отрезке  интерполяционным линейным сплайном с равномерным шагом. Будем опираться на материал, изложенный в п.п. 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;                     // границы отрезка