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

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

Процедура вычисления коэффициентов сплайна SPLINE_COEFFICIENT(x, y, N, A, B, C, D), в которой вычисляются коэффициенты трехдиагональной системы уравнений (6). В ней осуществляется привидение системы к верхней треугольной форме (прямая подстановка), нахождение (обратная подстановка) и расчет коэффициентов по
формулам (7).

Процедура SPLINE_VALUE(x, y, N, A, B, C, D, xx, P, P1), будет служить для вычисления в точке хх значения кубического сплайна и его производной . В процедуре находится рабочий интервал . Значения сплайна и его производная в процедуре вычисляются по схеме Гориера

,

, .

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

1.3. Метод наименьших квадратов.

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

Пусть задан набор экспериментальных точек с координатами , . Согласно методу наименьших квадратов аппроксимирующую функцию выбирают таким образом, чтобы была минимальна сумма

.

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

.

При этом коэффициенты полинома находятся из критерия наименьших квадратов как

.

Вычисляя соответствующие производные , , …, и приравнивая их к нулю, получаем так называемые нормальные уравнения:

,

,

………….

.

Это система линейных алгебраических уравнений относительно неизвестных . Для полинома степени получается система из уравнений. Количество различных сумм, которые надо вычислять для определения матрицы коэффициентов при неизвестных равно , для определения вектора-столбца свободных членов - . Решение такой системы нормальных уравнений может проводиться методом исключения Гаусса.

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

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

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

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

Текст программы

#include <math.h> //в файле math.h определена

//функция возведения в степень pow(double, double)

//и функция модуля числа fabs(double)

#include <iostream> //для ввода/вывода через стандартные потоки

//ввода - cin и вывода - cout

#include <stdio.h> //для использования функции sprintf()

//макросы – константы

//максимальное число точек

#define MAX_DATA_DOTS 100

//максимальная степень аппроксимирующего полинома

#define MAX_POLINOME_POWER 10

//массивы для хранения координат точек

double DotX[MAX_DATA_DOTS], DotY[MAX_DATA_DOTS];

//число точек

int DotsCount;

//степень аппроксимирующего полинома

int power;

//матрица коэффициентов для системы линейных

//уравнений

double st[MAX_POLINOME_POWER][MAX_POLINOME_POWER];

//столбец свободных членов и столбец решения

//соответственно

double fr[MAX_POLINOME_POWER], sol[MAX_POLINOME_POWER];

/* Используемые функции

DataInput() функция для ввода исходных данных

FormMatrices() функция, формирующая матрицы для решения

системы линейных уравнений методом Гаусса

GAUSS() функция, решающая систему линейных уравнений

ShowResult() функция вывода результатов

ShowTable() функция вывода таблицы результатов

main() функция, с которой начинается выполнение

программы

*/

//функция для ввода исходных данных

bool DataInput()

{

bool Err=0; //признак ошибки

//ввод числа точек

std::cout << "Enter the number of dots" << std::endl;

std::cin >> DotsCount;

//если пользователь ввел число точек больше MAX_DATA_DOTS,

//то программа завершится c сообщением об ошибке

if (DotsCount > MAX_DATA_DOTS)

{

Err=1;

goto Exit;

}

//Ввод координат точек

std::cout << "Enter coordinates of the dots" << std::endl;

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

{

std::cout << "Dot N"<<i+1<<":"<< std::endl;

std::cout << " X=";

std::cin >> DotX[i];

std::cout << " Y=";

std::cin >> DotY[i];

}

//Ввод степени аппроксимирующего полинома

std::cout << "Enter the power of a polynomial approximant" << std::endl;

std::cout<<"Power=";

std::cin >> power;

//если пользователь ввел степень больше MAX_POLINOME_POWER,

//то программа завершится c сообщением об ошибке

if (power > MAX_POLINOME_POWER)

{

Err=1;

goto Exit;

}

//Cтепень аппроксимирующего полинома не должна превышать

// число точек-1, иначе задача не имеет решения

if (power>DotsCount-1)

{

Err=1;

goto Exit;

}

Exit:

return Err;

}

//функция, формирующая матрицы для решения

//системы линейных уравнений методом Гаусса

void FormMatrices()

{

//обнуление матрицы коэффициентов

memset (st, 0, MAX_POLINOME_POWER * MAX_POLINOME_POWER);

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

{

//Задание матрицы коэффициентов

for (int j = 0; j <= power; j++)

for (int k = 0; k < DotsCount; k++)

st[i][j] += pow(DotX[k],(i+j));

//Задание столбца свободных членов

fr[i]=0;

for (int k = 0; k < DotsCount; k++)

fr[i] += (DotY[k] * pow(DotX[k],i));

}

}

//Функция, решающая систему линейных уравнений

void GAUSS()

{

int i,j,k;

double tmp;

//сканирование матрицы по столбцам

for(k = 0; k < power; k++)

{

//сканирование по элементам k-го столбца,

//обнуление всех элементов столбца ниже k-го

//путём умножения i-й строки на соответствующий коэффициент

for(i = k + 1; i < power+1; i++)

{

tmp = st[i][k] / st[k][k];

for (j = k + 1; j < power+1; j++)

st[i][j] -= tmp * st[k][j];

fr[i] -= tmp * fr[k];

}

}

sol[power] = fr[power] / st[power][power];

for (i = power - 1; i >= 0; i--)

{

tmp = 0;

for (j = i + 1; j<power+1; j++)

tmp += st[i][j] * sol[j];

sol[i] = (fr[i] - tmp) / st[i][i];

}

}

//Функция вывода результатов

void ShowResult()

{

std::cout <<std::endl<< "List of the polynomial coefficients" << std::endl;

std::cout <<"N Value" << std::endl;

int j;

for (j = 0; j < power+1; j++)

{

std::cout <<j<<" "<<sol[j]<< std::endl;

}

std::cout <<std::endl<< "Polynomial approximant:" << std::endl;

//массив символов answer содержит ответ - полином

char chtemp[20], answer[256];

sprintf(chtemp, "%6.2f", sol[0]);

strcpy (answer, "f(x) = ");

strcat (answer, chtemp);

for (j = 1; j < power+1; j++)

{

sprintf(chtemp," + (%.2f) * x^ %d", sol[j], j);

strcat (answer,chtemp);

}

std::cout << answer << std::endl;

}

//Функция, выводящая таблицу результатов:

//координаты точек и значения

//аппроксимирующего полинома в соответствующих точках

void ShowTable()

{

//Вывести заголовок таблицы

//char(9) - символ табуляции (смещения)

std::cout <<std::endl<< "X" <<char(9)<<"Y";

std::cout<<char(9)<<"Yp"<<char(9)<< std::endl;int i,j;

double x,P; //в переменной P будет сформировано

//значение полинома в точке х

//цикл по всем точкам

for (i=0; i<DotsCount;i++)

{

//Вычисление значения аппроксимирующего полинома

//в точке х=DotX[i]

P=0;

x=DotX[i];

for (j=power;j>=0;j--)

{

P=P*x+sol[j];

}

//Вывод строки таблицы

std::cout<<DotX[i]<<char(9)<<DotY[i]<<char(9)<<P<<std::endl;

}

}

//С этой функции начинается выполнение программы

void main()

{

//ввод с клавиатуры исходных данных

if (DataInput()!=0) goto Error;

//формирование матриц для решения

//системы линейных уравнений методом Гаусса

FormMatrices();

//Решение системы

GAUSS();

//Полученное решение - столбец полиномиальных коэффициентов

//Вывод результата

ShowResult();

//Вывести таблицу: координаты точек и значения

//аппроксимирующего полинома в соответствующих точках

ShowTable();

//Если ошибок не было, программа нормально завершится,

//в противном случае будет выведено сообщение об ошибке

goto OK;

Error:

std::cout <<std::endl<< "Error!" << std::endl;

OK:

//Этот код позволяет просмотреть результаты прежде, чем завершится

//исполнение программы

int dummy; //переменная, не имеющая смысловой нагрузки

std::cout <<std::endl<< "End of the program. Enter any number ";

std::cin >> dummy; //программа не завершится, пока

//пользователь не введёт произвольное число

}

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

(схема Горнера)

Такая запись позволяет минимизировать число выполняемых арифметических операций, избегая возведения в степень. Вычисление значений полинома по схеме Горнера может быть реализовано с помощью следующих операторов:

P=0;

FOR (j = power;j>=0; j--)

P=P*x+ sol[j];

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

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

2) вызов процедуры FormMatrices для формирования системы нормальных уравнений;

3) вызов процедуры GAUSS для решения системы линейных уравнений;

4) вычисление значений аппроксимирующего полинома в узловых точках;

5) вывод результатов: коэффициентов полинома, исходных значений (х_к, у_к) и значений полинома в соответствующих точках.

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

Поиск по сайту: