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

Из курса математики известны 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

#include

#include

#include

#include

#include

#include

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

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

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: «<

double y, temp;

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

y=0;

for (int t=0; 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

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

for (k=0; k

AKK=0;

for (i=k; i

l=IOR[i];

if (fabs (array[l] [k])

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

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

for (i=k+1; i

l=IOR[i];

delitel=array[l] [k];

for (j=k; 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

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

double *b=new double [m];

double s=0;

for (int j=0; j

s=0;

for (int i=0; i

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

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

}

b[j]=s;

for (int k=0; k

s=0;

for (i=0; 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

double *b=new double [m];

double s=0;

for (int j=0; j

s=0;

for (int i=0; i

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

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

}

b[j]=s;

for (int k=0; k

s=0;

for (i=0; 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

for (i=0; i

for (i=0; i

for (i=0; i

for (i=0; i

for (i=0; i

for (i=0; i

for (i=0; i

for (i=0; i

for (i=0; i

for (i=0; i

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

cout<< «x1="<<

cout<< «x2="<<

cout<< «x3="<<

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

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;

}