Аппроксимация функций

Государственное образовательное  учреждение

Высшего профессионального  образования

Уфимский государственный  нефтяной технический университет

Филиал ГОУ ВПО УГНТУ  в г. Стерлитамаке

 

 

 

 

Кафедра АТИС

 

 

 

Курсовая работа

По ПОА

Вариант 14.

 

 

 

 

 

Выполнил:        студент группы АК-09-31

Мухаметгалина А. А.

Проверил:       к.ф-м.н. , доцент

Биккулов И.М.

 

 

 

 

 

 

 

Стерлитамак

2010

 

Содержание

1. Введение          3
2. Задание          4
3. Алгоритм аппроксимации функции степенным регрессионным полиномом М-ого порядка       5
4. Блок-схема основной программы.      7
5. Листинг программы, реализующей предложенный алгоритм

на языке С++.         13

6. Блок-схема программы, вычисляющей степенной полином с

помощью схемы Горнера.       19

7. Листинг программы, вычисляющей степенной полином с

помощью схемы Горнера.       20

8. Заключение.         22
9. Список используемой литературы.      23

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Введение

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

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

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

Данная работа связана  с аналитической градуировкой датчиков, сопряженных с управляющей вычислительной машиной (УВМ).

Функциональная зависимость  между измеряемой величиной Y и выходным сигналом датчика Х в общем  случае определяется зависимостью:

Для задач контроля и  управления необходимо знать истинное значение измеряемой величины Y. При  неизвестном значениее Х, Y может  быть найден как:

X=(F)^-1(X)=f(X).

Если датчик имеет  линейную характеристику Х=ay+b,то определить Y легко:

.

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

 

 

ЗАДАНИЕ

 

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

Y – измеряемая величина,

X – выход термопары,

Еmax – максимально-допустимая ошибка аппроксимации.

Таблица 1 – Значения для варианта №14

№ п\п	Y	X
1	0	0,00
2	20	2,96
3	50	7,81
4	60	10,32
5	75	12,01
Emax	Допустимая среднеквадратичная ошибка аппроксимации	1

 

В соответствии с вариантом  задания необходимо:

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

2. Разработать программу, реализующую предложенный алгоритм на языке С++.

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

4. Используя алгоритм схемы Горнера, разработать программу вычисления полученного степенного полинома.

5. Пояснительную записку к курсовой работе оформить с помощью текстового редактора типа “Word”.

 

Алгоритм аппроксимации  табличной функции степенным  регрессионным полиномом М-ного порядка (по методу наименьших квадратов).

 

Опишем алгоритм, использую средства MathCAD 14.

 

Дискретные значения X и Y определяются значениями:

Аппроксимируем заданную табличную функцию простейшим полиномом  вида:

P(X) = a₀ + a₁X

и выделим основные этапы  расчета:

1. Формируем матрицу Х, размерности 5 х 2:

2. Получаем транспонированную  матрицу:

3. Вычисляем произведения X^TX и X^TY:

4. Вычисляем обратную  матрицу:

5. Вычисляем вектор  коэффициентов:

Таким образом, аппроксимирующий полином будет иметь вид:

P(X) = –0,13122 + 0,16466 X

Для оценки ошибки аппроксимации  воспользуемся выражением:

Если рассчитанная ошибка аппроксимации окажется больше допустимой ошибки Emax, то следует увеличить степень полинома до выполнения неравенства E ≤ Emax. При этом следует иметь ввиду, что повышать степень полинома можно только до тех пор, пока m+1 ≤ N.

 

Блок-схема основной программы

 

Блок-схема функции вычисления определителя матрицы

 

Блок-схема функции вычеркивания строки и столбца

 
Программа, реализующая предложенный алгоритм на языке С++.

// kursov1.cpp: определяет точку входа для  консольного приложения.

#include "stdafx.h"

#include <math.h>

#include <conio.h>

#include <windows.h>

#include <stdio.h>

#include <iostream>

using namespace std;

 

double tabl_X[5] = { 0, 20, 50, 60, 75 };    //исходные данные X

double tabl_Y[5] = { 0, 2.96, 7.81, 10.32, 12.01 };  //исходные данные Y

 

const int N = 5; //число точек

const double Emax=1;//допустимая среднеквадратичная ошибка аппроксимации

int m;    //ширина матрицы Х, длина вектора А (m <= N)

int i, j;   //счетчики

double **X;   //матрица Х - рабочая

double **Xt;  //матрица Х - транспонированная

double **XtX;  //произведение XtX

double *XtY;  //произведение XtY

double **XtXb;  //обратная XtX

double *A;   //вектор коэффициентов

double e[5];  //вектор квадратичных ошибок

double E;   //ошибка аппроксимации

 

//Функция печати  матрицы

void print_matr(double **X, int n, int m, int k){

cout.precision(k);  

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

for(j = 0; j < m; j++)

cout << X[i][j] << "  ";

cout << endl;

}

}//========================================================

 

//Функция печати вектора

void print_vect(double *X, int n, int k){

cout.precision(k);  //число выводимых после точки знаков

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

cout << X[i] << "  ";

cout << endl;

}//========================================================

 

//Функция  освобождения памяти из-под матрицы

void free_mem(double **X, int n){

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

delete[] X[i];

delete[] X;

}//========================================================

 

//функция вычеркивания  строки и столбца

void get_matr(double **matr, int n, double **temp_matr, int indRow, int indCol){

int ki = 0, kj;

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

if(i != indRow){ 

kj = 0;

for (j = 0; j < n; j++){

if (j != indCol){ 

temp_matr[ki][kj] = matr[i][j];  

kj++;

} }

ki++;          

} }

}//========================================================

 

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

double determ(double **X, int m){

double sum = 0;  

int n;

if (m == 1)

sum = X[0][0]; //расчет определителя для матрцицы единичной размерноси

else if (m == 2)

sum = X[0][0] * X[1][1] - X[1][0] * X[0][1]; //расчет определителя для матрцицы размерноси 2

else{

for(int i = 0; i < m; i++){  //расчет определителя для матрцицы размерности больше двух

n = m - 1;

double **temp_matr = new double*[n]; 

for(j = 0; j < n; j++)

temp_matr[j] = new double[n];

get_matr(X, m, temp_matr, 0, i); 

double det = determ(temp_matr, n); //Вызывается функция расчета определителя

sum = sum + pow(-1.0,2+i) * X[0][i] * det; //расчет непосредственно определителя

free_mem(temp_matr, n); //освобождаем память из под матрицы temp_matr

} }

return sum;

}//========================================================

 

int _tmain(int argc, _TCHAR* argv[])

{

m = 1; 

do

{ 

m++; //на первой итерации степень полинома равна 1, соответственно длина вектора А равна 2

 

//формируем матрицу Х

X = new double*[N];

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

X[i] = new double[m];

 

for(i = 0; i < N; i++){ 

X[i][0] = 1;   

X[i][1] = tabl_X[i]; 

if(m > 2)

for(j = 2; j < m; j++)

X[i][j] = pow(tabl_X[i], j); 

}

 

//получаем транспонированную матрицу X

Xt = new double*[m];

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

Xt[i] = new double[N];

 

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

for(j = 0; j < N; j++)

Xt[i][j] = X[j][i]; 

 

//вычисляем произведения XtX и XtY

XtX = new double*[m];

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

XtX[i] = new double[m];

XtY = new double[m];

 

double sum; 

int k;

for(i = 0; i < m; i++){

for(j = 0; j < m; j++){

sum = 0;

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

sum = sum + Xt[i][k]*X[k][j]; //перемножение строки Хt на столбец Х

XtX[i][j] = sum; 

}

sum = 0;

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

sum = sum + Xt[i][k]*tabl_Y[k]; //перемножение Xt на вектор Y

XtY[i] = sum;

}

 

//вычисляем обратную  матрицу

XtXb = new double*[m];

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

XtXb[i] = new double[m];

 

double det = determ(XtX, m); //Вызывается функция расчета определителя

if(det != 0){ 

for(int i = 0; i < m; i++){

for(int j = 0; j < m; j++){

int n = m - 1;  

double **temp_matr = new double*[n]; 

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

temp_matr[k] = new double[n];

get_matr(XtX, m, temp_matr, i, j); 

double det1 = determ(temp_matr, n);

XtXb[i][j] = pow(-1.0, i + j + 2) * det1 / det; //рассчитываем элементы обратной матрицы XtXb

free_mem(temp_matr, n); //освобождаем память из под матрицы temp_matr

}

}      

}

else

cout << "Невозможно найти обратную матрицу" << endl;

       

//вычисляем вектор  коэффициентов

A = new double[m];

for(i = 0; i < m; i++){

sum = 0;

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

sum = sum + XtXb[i][k]*XtY[k]; //вектор А рассчитывается как произведение матрицы XtXb на вектор XtY

A[i] = sum;

}

 

//Вычисляем ошибку  аппроксимации

for(i = 0; i < N; i++){

sum = 0;

for(j = 0; j < m; j++)

sum = sum + pow(tabl_X[i],j) * A[j];

e[i] = pow(tabl_Y[i] - sum, 2);  //расчет вектора квадратичных ошибок

}

sum = 0;

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

sum = sum + e[i];

E = pow(sum/(1.0*N), 0.5); //расчет среднеквадратичной ошибки аппроксимации

 

}

while((E > Emax) && (m < N)); 

 

setlocale(LC_ALL, ".ACP");  

printf("Исходные данные:\n\n");

printf("  X |   Y  \n");

printf("-----------\n");

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

printf(" %2.0lf | %5.2lf  \n",tabl_X[i],tabl_Y[i]);

printf("\nEmax = %.2lf\n\n",Emax);

printf("============================================\n\n");

 

cout << fixed;

 

printf("Результаты расчетов:\n\n");

printf("Порядок регрессионного полинома M = %d\n\n",m-1);

printf("Матрица Х:\n");

print_matr(X,N,m,0);

printf("\nМатрица Хt - транспонированная Х:\n");

print_matr(Xt,m,N,0);

printf("\nМатрица ХtX = Xt * X:\n");

print_matr(XtX,m,m,0);

printf("\nВектор ХtY = Xt * tabl_Y:\n");

print_vect(XtY,m,3);

printf("\nМатрица ХtYb - матрица, обратная XtY:\n");

print_matr(XtXb,m,m,7);

printf("\nВектор А - вектор коэффициентов регрессионного полинома:\n");

print_vect(A,m,9);

printf("\nCреднеквадратичная ошибка аппроксимации E = %.5lf\n",E);

 

_getch();

return 0;

}

 

 

 

 

 

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

 

Как видим, ошибка E = 0.32, причем E < Emax. Значит первой степени полинома для аппроксимации достаточно.

 

Используя алгоритм схемы  Горнера, разработаем программу  вычисления полученного степенного полинома.

 

Вычисление степенных  полиномов на ЭВМ обычно осуществляется по формуле Горнера, которая для полинома m-го порядка в виде:

y = P(X) = (…((a₁X + a₂)X + a₃)X + … + a_m)X + a_m+1

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

Y = a_kX + a_k+1, где k = 1, 2, …, m.

Коэффициенты полинома сводятся в массив, включающий (m+1) элемент. Начальное значение переменной Y, задаваемое перед циклом, должно быть равно коэффициенту a₁ при Х в старшей степени, а параметр цикла должен изменяться от 2 до m+1.

Следует отметить, что  при использовании данного алгоритма  степенной полином следует переписать в виде:

P(X) = a₁X^m + a₂X^m-1 + a₃X^m-2 + … + a_mX + a_m+1.

 

Блок-схема основной программы

 

 

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

// kursov2.cpp: определяет точку входа для  консольного приложения.

#include "stdafx.h"

#include <conio.h>

#include <windows.h>

#include <stdio.h>

#include <iostream>

using namespace std;

 

double tabl_X[5] = { 0, 20, 50, 60, 75 };   //исходные данные X

double tabl_Y[5] = { 0, 2.96, 7.81, 10.32, 12.01 };  //исходные данные Y

 

const int N = 5; //число точек

int m;   //длина вектора А (m <= N)

int i, j;  //счетчики

double *A;  //вектор коэффициентов

double Y[5];  //вектор расчетных Y

 

int _tmain(int argc, _TCHAR* argv[])

{

setlocale(LC_ALL, ".ACP"); //установка поддержки русского языка для вывода на экран