Методы аппроксимаций функций

/

/

Введение

Из курса математики известны 3 способа задания функциональных зависимостей: 1) аналитический 2) графический 3) табличный. Табличный способ обычно возникает в результате эксперимента. Недостаток табличного задания функции заключается в том, что найдутся значения переменных которые неопределенны таблицей. Для отыскания таких значений определяют приближающуюся к заданной функцию, называемой аппроксимирующей, а действие замены аппроксимацией. Аппроксимация заключается в том что, используя имеющуюся информацию по f(x) можно рассмотреть другую функцию близкую в некотором смысле к f(x), позволяющую выполнить над ней соответствующие операции и получить оценку погрешность такой замены.

1. Постановка задачи

В ходе физического эксперимента получены следующие данные с интервалом 1с (первое наблюдение выполнено в момент t=1,0):

Подберите аппроксимирующий полином методом наименьших квадратов. Аппроксимируйте данные моделью:

Проанализируйте результаты.

2. Теоретические предпосылки

2.1 Аппроксимация функций методом наименьших квадратов

Пусть требуется аппроксимировать функцию заданную таблично, в точках x_i, ее значения в этих точках обозначим y_i.Основная идея метода наименьших квадратов - построить функцию , - вектор параметров, проходящую наиболее близко к заданной системе точек, не проходящую в общем случае через узлы таблицы, а именно, так чтобы сумма квадратов отклонений таблично заданной функции от искомой функции в узлах была минимальной. Для этого вводят критерий метода, следующим образом:

наименьший квадрат алгоритм mathcad

т.е. значение среднеквадратического отклонения должно быть минимальным.

Параметры подбирают так, чтобы выполнялся критерий.

Чаще всего используются 2-3 параметрические функции, следующих семейств:

, , , ,

дробно-линейная,

логарифмическая,

обратно-пропорциональная,

дробно-рациональная,

где a, b, c - параметры

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

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

Чтобы получить окончательный вид аппроксимирующего многочлена (2.33), нужно найти коэффициенты , выбрать вид базисных функций и определить степень обобщенного многочлена m так, чтобы среднеквадратическая погрешность (2.32) была наименьшей. При выборе степени следует руководствоваться тем, что значение среднеквадратичного отклонения (обозначим как ) с ростом степени m сначала убывает, а затем начинает возрастать. Отсюда правило выбора: за оптимальное значение степени многочлена следует принять то значение m, начиная с которого величина стабилизируется или начинает возрастать.

Применение степенных базисных функций

Очень часто для приближения по методу наименьших квадратов используются алгебраические многочлены степени mn, т.е. k=0..m. В этом случае обобщенный многочлен (2.33) имеет вид

, а коэффициенты , исходя из критерия (2.32), ищутся из системы следующего вида:

или в развернутом виде:

В случае многочлена первой степени P₁(x)=c₀+c₁x, эта система принимает вид

Для многочлена второй степени P₂(x)=c₀+c₁x+c₂x²:

3. Описание программного средства

3.1 Спецификация переменных, процедур и функций

Void XY_print() - Функция очищает экран и рисует координатные оси.

double P (double *vek, int m, double xx) - функция для вычисления значения многочлена.

double Q (double *a, int m) - функция, возвращающая среднеквадратическое отклонение со входным параметром - вектором значений *a.

void draw (double *X, int col, int count, double q) - функция, выводящая график функции аппроксимированной по методу наименьших квадратов полиномом; входные данные - вектор значений *X, цвет col, целое число count и значение среднеквадратического отклонения q.

double lsolve (double **array, double *vek, int nn, double *X) - функция, решающая СЛАУ методом Гаусса.

double mnk (double m, int color) - функция, аппроксимирующая функцию y=f(x) методом наименьших квадратов

_mnk() - Функция аппроксимирует данные моделью: y(t)=x1+x2*t+x3*sin(t) по методу наименьших квадратов.

void menu() - функция вывода меню на экран.

int main() - главная функция.

3.2 Схемы алгоритмов

double mnk (double m, int color)

/

/

/

/

double _mnk()

/

/

/

/

4. Расчеты в системы Mathcad

4.1 Контрольный пример работы

Рисунок 1 - Метод наименьших квадратов

Рисунок 2 - Аппроксимация моделью: , по методу наименьших квадратов.

4.2 График в среде в программирования

Рисунок 3 - Метод наименьших квадратов

Рисунок 4 - Метод наименьших квадратов моделью

Заключение

Проанализировав решение системы в MathCAD, можно сделать вывод, что решение, полученное с помощью программы правильное, так как графики практически совпадают по форме и по значениям.

Список использованных источников

1. Амосов А.А., Дубинский Ю.А. Вычислительные методы для инженеров: Учеб. Пособие. - М.:Высш.шк., 1994. - 544 с.:ил.

2. Лапчик М.П. Численные методы: Учебное пособие для вузов / Рагулина М.И., Хеннер Е.К.; под редакцией Лапчика М.П.: Издательский центр «Академия», 2004. - 384 с.

3. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. - Численные методы (2 изд.).

4. Тарасов В.Н. Численные методы, теория, алгоритмы, программы: Учебное пособие для вузов / Бахарева Н.Ф; под редакцией Тарасова В.П.: Издательский центр Оренбург: ИПК ОГУ, 2003. - 100 с.

Приложение А

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

#include <iostream.h>

#include <graphics.h>

#include <conio.h>

#include <stdio.h>

#include <string.h>

#include <math.h>

#include <stdlib.h>

void init_graph() {

int gdriver = DETECT, gmode, errorcode;

initgraph (&gdriver, &gmode, «»);

errorcode = graphresult();

if (errorcode!= grOk) {

printf («Graphics error:%sn», grapherrormsg(errorcode));

printf («Press any key to halt:»);

getch();

exit(1);

}

}

const n=25;

double x[n]={1. 0,2. 0,3. 0,4. 0,5. 0,6. 0,7. 0,8. 0,9. 0,10. 0,11. 0,12. 0,13. 0,14.0,

15. 0,16. 0,17. 0,18. 0,19. 0,20. 0,21. 0,22. 0,23. 0,24. 0,25.0};

double y[n]={5. 0291,6. 5099,5. 3666,4. 1272,4. 2948,6. 1261,12. 5140,10. 0502,9. 1614,7. 5677,7. 2920,10. 0357,11.0708,

13. 4045,12. 8415,11. 9666,11. 0765,11. 7774,14. 5701,17. 0440,17. 0398,15. 9069,15. 4850,15. 5112,17.6572};

int menu_output()

{

cout<<» г============================================================================¬nr»;

cout<<» ¦ 0. Approx. model (y(t)=x1+x2*t+x3*sin(t)) ¦nr»;

cout<<» ¦ 1. Real grafik ¦nr»;

cout<<» ¦ 2. Grafik po MNK ¦nr»;

cout<<» ¦ 3. Ochistit' ¦nr»;

cout<<» ¦============================================================================¦nr»;

cout<<» ¦ ESC - Exit ¦nr»;

cout<<» L============================================================================-nr»;

return 0;

}

void XY_print() { // Вырисовка осей

setlinestyle (0,0,1);

setbkcolor(0);

setfillstyle (1,0);

setcolor(WHITE);

bar (0,0,680,459);

line (20,380,320,380);

line (160,400,160,180);

line (320,380,310,375);

line (320,380,310,385);

line (160,180,155,190);

line (160,180,165,190);

outtextxy (145,180, «Y»);

outtextxy (320,390, «X»);

outtextxy (145,390, «0»);

}

void p_pixel (double x, double y, int col) {

putpixel (x+160, - y+380, col);

}

double P (double *vek, int m, double xx) { // Значение многочлена

double tem=0;

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

if(! i) tem+=vek[i];

else tem+=vek[i]*pow (xx, i);

return tem;

}

double Q (double *a, int m) { // значение среднеквадратичного отклонение многочлена P

double temp=0;

for (int i=0; i<n; i++) temp+=pow((P (a, m, x[i]) - y[i]), 2);

temp=temp/n;

temp=sqrt(temp);

temp=ddd [m-1];

return temp;

}

void draw (double *X, int col, int count, double q) { // рисует график по MNK

setfillstyle (1,0);

cout<<«Otklonenie: «<<q<<endl;

double y, temp;

for (double xx=1; xx<=25; xx+=0.001) {

y=0;

for (int t=0; t<count; t++)

if(! t) y+=X[t];

else y+=X[t]*pow (xx, t);

p_pixel (xx*6, y*6, col);

}

}

 // -Решение СЛАУ методом Гауса-

double lsolve (double **array, double *vek, int nn, double *X) {

int *IOR=new int[nn];

int i, l, k, M, j, p;

double AKK, AMAIN, delitel;

for (k=0; k<nn; k++) IOR[k]=k;

 // -Нахождение индекса-

for (k=0; k<nn; k++) {

AKK=0;

for (i=k; i<nn; i++) {

l=IOR[i];

if (fabs (array[l] [k])<AKK) continue;

M=l; p=i; AKK=fabs (array[l] [k]);

}

 // -Меняем места-

AMAIN=array[M] [k];

if (k!=p) {IOR[p]=IOR[k]; IOR[k]=M;};

 // -Прямой ход-

for (j=k; j<nn; j++) array[M] [j]=array[M] [j]/AMAIN;

vek[M]=vek[M]/AMAIN;

for (i=k+1; i<nn; i++) {

l=IOR[i];

delitel=array[l] [k];

for (j=k; j<nn; j++) array[l] [j]=array[l] [j] - delitel*array[M] [j];

vek[l]=vek[l] - delitel*vek[M];

}

}

 // -Обратный ход-

double sum;

for (k=nn-1; k>=0; k-) {

l=IOR[k];

sum=0;

for (j=k+1; j<nn; j++) sum=sum+array[l] [j]*X[j];

X[k]=vek[l] - sum;

}

delete[] IOR;

return *X;

}

double mnk (double m, int color) { // Наименьших квадратов

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

for (int i=0; i<n; i++) A[i]=new double [m];

double *b=new double [m];

double s=0;

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

s=0;

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

if(! j) s+=y[i];

else s+=y[i]*pow (x[i], j);

}

b[j]=s;

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

s=0;

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

if(! (k+j)) s+=1;

else s+=pow (x[i], k+j);

A[j] [k]=s;

}

}

double *X=new double [m];

*X=lsolve (A, b, m, X); //if (m==25) {X=a25;}

draw (X, color, m, Q (X, 25));

delete[] A;

delete[] X;

delete[] b;

return 0;

}

double _mnk() { // Аппроксимация данных моделью (y(t)=x1+x2*t+x3*sin(t) по mnk

int m=13;

double *X=new double[m];

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

for (int i=0; i<n; i++) A[i]=new double [m];

double *b=new double [m];

double s=0;

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

s=0;

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

if(! j) s+=y[i];

else s+=y[i]*pow (x[i], j);

}

b[j]=s;

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

s=0;

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

if(! (k+j)) s+=1;

else s+=pow (x[i], k+j);

A[j] [k]=s;

}

}

*X=lsolve (A, b, m, X);

clrscr();

XY_print();

menu_output();

double **matr=new double*[3];

for (i=0; i<3; i++) matr[i]=new double[3];

double *bb=new double[3];

double *cc=new double[3];

matr[0] [0]=1;

for (i=0; i<n; i++) matr[0] [1]+=x[i];

for (i=0; i<n; i++) matr[0] [2]+=sin (x[i]);

for (i=0; i<n; i++) b[0]+=p (A, 15, i);

for (i=0; i<n; i++) matr[1] [0]+=x[i];

for (i=0; i<n; i++) matr[1] [1]+=x[i]*x[i];

for (i=0; i<n; i++) matr[1] [2]+=x[i]*sin (x[i]);

for (i=0; i<n; i++) bb[1]+=x[i]*p (A, 15, i);

for (i=0; i<n; i++) matr[2] [0]+=sin (x[i]);

for (i=0; i<n; i++) matr[2] [1]+=x[i]*sin (x[i]);

for (i=0; i<n; i++) matr[2] [2]+=sin (x[i])*sin (x[i]);

for (i=0; i<n; i++) bb[2]+=p (A, 15, i)*sin (x[i]);

*cc=lsolve (matr, bb, 3, cc);

cout<< «x1=»<<cc[0]<<endl;

cout<< «x2=»<<cc[1]<<endl;

cout<< «x3=»<<cc[2]<<endl;

double yy, yy1;

yy1=0;

for (double xx=1; xx<=25; xx+=1) {

yy=cc[0]+cc[1]*xx+cc[2]*sin(xx);

line (xx*6+160, - yy*6+380, (xx-1)*6+160, - yy1*6+380);

yy1=yy;

p_pixel (xx*6, yy*6,2);

}

delete[] A;

delete[] X;

delete[] b;

return 0;

}

void real_graph() {

setcolor(4);

setfillstyle (1,0);

setlinestyle (0,0,3);

setcolor(RED);

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

circle (x[i]*6+160, - y[i]*6+380,1);

}

int main() {

char ch;

clrscr();

init_graph();

XY_print();

menu_output();

do

{

ch=getch();

switch(ch)

{

case '0':

_mnk();

break;

case '1':

real_graph();

break;

case '2':

mnk (25,2);

break;

case '3':

clrscr();

XY_print();

menu_output();

break;

default:

break;

}

if (ch==27)

{

exit(0);

}

} while (ch!=27);

return 0;

}