close

Вход

Забыли?

вход по аккаунту

?

otchet(2)

код для вставкиСкачать
МГТУ им. Н.Э. Баумана
Лабораторная работа №2
"Метод наименьших квадратов"
Работу выполнила: Глагазина Н.
Группа: МТ4-31
Вариант №3
Москва 2013 г.
Теоретическая часть
Метод наименьших квадратов, один из методов теории ошибок для оценки неизвестных величин по результатам измерений, содержащим случайные ошибки. Он применяется также для приближённого представления заданной функции другими (более простыми) функциями. Задача заключается в нахождении коэффициентов зависимости, при которых функция: F = ∑_(i=0)^n▒( φ(a_0, a_1, ..., a_m,x_i ) - y_i )^2 = ∑_(i=0)^n▒ε_i^2 принимает наименьшее значение, причем
yi - функция значения которой известны в узлах xi , при i = 0, 1, ..., n;
ϕ(a0, a1,..., am,xi) - функция, зависящая от параметров a0, a1, ..., am. Величину εi называют невязкой.
Для начала найдём частные производные функции F по переменным ai, приравняем эти производные к нулю:
∂F/(∂a_i ) = 0, i = 1, 2, ..., m.
При ϕ(x) = amxm + am - 1xm - 1 + ... + a1x1 + a0 (1) получим следующую СЛАУ:
∂F/(∂a_0 )=2∑_(i=0)^n▒〖(a_m x_i^m + a_(m - 1) 〗 x_i^(m - 1) + ... + a_1 x_i^1 + a_0) = 0,
∂F/(∂a_1 )=2∑_(i=0)^n▒〖(a_m x_i^m + a_(m - 1) 〗 x_i^(m - 1) + ... + a_1 x_i^1 + a_0)x_i = 0,
∂F/(∂a_2 )=2∑_(i=0)^n▒〖(a_m x_i^m + a_(m - 1) 〗 x_i^(m - 1) + ... + a_1 x_i^1 + a_0)x_i^2 = 0,
.................................................................................
∂F/(∂a_m )=2∑_(i=0)^n▒〖(a_m x_i^m + a_(m - 1) 〗 x_i^(m - 1) + ... + a_1 x_i^1 + a_0)x_i^m = 0.
Преобразуем её к виду:
(n+1) a_0 + (∑_(i=0)^n▒x_i ) a_1 + (∑_(i=0)^n▒x_i^2 ) a_2 + ... + (∑_(i=0)^n▒x_i^m ) a_m= ∑_(i=0)^n▒y_i ,
(∑_(i=0)^n▒x_i ) a_0 + (∑_(i=0)^n▒x_i^2 ) a_1 + (∑_(i=0)^n▒x_i^3 ) a_2 + ... + (∑_(i=0)^n▒x_i^(m + 1) ) a_m= ∑_(i=0)^n▒〖x_i y〗_i ,
............................................................................................
(∑_(i=0)^n▒x_i^m ) a_0 + (∑_(i=0)^n▒x_i^(m + 1) ) a_1 + (∑_(i=0)^n▒x_i^(m + 2) ) a_2 + ... + (∑_(i=0)^n▒x_i^2m ) a_m= ∑_(i=0)^n▒〖〖x_i^m y〗_i.〗
Затем необходимо ввести коэффициенты:
b_pq= ∑_(i=0)^n▒〖x_i^(p+q),〗
c_p= ∑_(i=0)^n▒x_i^p y_i.
После проведенных операций получим следующую СЛАУ:
b00a0 + b01a1 + ... + b0mam = c0,
b10a0 + b01a1 + ... + b1mam = c1,
.............................................................
bm0a0 + bm1a1 + ... + bmmam = cm.
Решив полученную СЛАУ любым методом (например, методом Гаусса) и подставив полученные значения в выражение (1), получим функцию, на графике которой минимизируется сумма квадратов невязок.
Задача.
1. Написать программу, которая строит таблицу значений y_i=f(a+ih) на отрезке [a;b] с шагом h = (b - a)/n, n=10. По полученной таблице методом наименьших квадратов найти линейную функцию, параболу и кубическую функцию, на которых минимизируется сумма квадратов невязок. 2. Результаты программы оформить в виде таблицы, столбиками которой являются:
1) значения x_i=a+ih, i=1, 2,..., n;
2) значения заданной функции y_i=f(x), i=1, 2, ..., n;
3) значения получившейся линейной функции φ_1 (x)=a_0+a_1 x в точках x_i, i=1, 2,...,n;
4) значения получившейся параболы φ_2 (x)=a_0+a_1 x+a_2 x^2 в точках
x_i, i=1, 2,...,n;
5) значения получившейся кубической функции φ_3 (x)=a_0+a_1 x+a_2 x^2+a_3 x^3 в точках x_i, i=1, 2,...,n;
6) три столбца значений невязок для линейной, квадратичной и кубической функций в точках x_i, i=1, 2,...,n;
7) суммарная невязка в нижней строке для соответствующих столбцов.
Функция f(x)=sh x на отрезке [0,2].
Текст программы.
void gauss(double mas[][N+1], int nn, int otv[],double x[]);
void glavelem( int k, double mas[] [N + 1], int nn, int otv[] );
void main()
{
double a=-1;
double b=3;
double bc[R][R]; // массивы коэф. b
double bc2[R+1][R+1];
double bc3[R+2][R+2];
double ac[R]; // массивы коэф. a
double ac2[R+1];
double ac3[R+2];
double cc[R]; // массивы коэф. c
double cc2[R+1];
double cc3[R+2];
// массивы для Гаусса
double mas[N][N+1] ;
double x[N];
//
double yi[N]; // yi
double xi[N]; // xi
double linfunc[N]; //значения линейной функции в i-тых точках
double kvadfunc[N]; //значения квадратичной функции в i-тых точках
double kubfunc[N]; //значения кубической функции в i-тых точках
double elin[N]; //значения невязок линейной функции в i-тых точках
double ekvad[N]; //значения невязок линейной функции в i-тых точках
double ekub[N]; //значения невязок линейной функции в i-тых точках
int i, j, k, nn, r;
double sx=0;
for(int i=0;i<N;i++)
{
xi[i]=a + i*(b-a)/10;
yi[i]=atan(2*x[i]+3);
sx=sx+xi[i];
}
//находим коэффициенты b и c
double Sb=0, Sc=0;
for(int p=0;p<2;p++)
{
for(int q=0;q<2;q++)
{
Sc=0;
Sb=0;
for(i=0;i<N;i++)
{
Sb=pow(xi[i],(p+q)) + Sb;
Sc=pow(xi[i], p)*yi[i] + Sc;
}
bc[p][q]=Sb;
}
cc[p]=Sc;
}
//находим коэффициенты а методом Гаусса
nn=2;
for ( i = 0; i < nn; i++ )
for ( j = 0; j < nn + 1; j++ )
{
if (j<nn)
{mas[i] [j] = bc[i] [j];}
else
{mas[i] [j] = cc[i];}
}
int otv[N]; //Отвечает за порядок корней
//Сначала все корни по порядку
for ( i = 0; i < nn + 1; i++ )
otv[i] = i;
gauss(mas, nn, otv, x);//вызов метода Гаусса
r=0;
for ( i = 0; i < nn; i++ )
for ( j = 0; j < nn; j++ )
if ( i == otv[j] )
{ //Расставляем корни по порядку
ac[r]=x[j];
r=r+1;
break;
}
//находим коэффициенты b2 и c2
double Sb2=0, Sc2=0;
for(int p=0;p<3;p++)
{
for(int q=0;q<3;q++)
{
Sc2=0;
Sb2=0;
for(i=0;i<N;i++)
{
Sb2=pow(xi[i],(p+q)) + Sb2;
Sc2=pow(xi[i], p)*yi[i] + Sc2;
}
bc2[p][q]=Sb2;
}
cc2[p]=Sc2;
}
//находим коэффициенты а2 методом Гаусса
nn=3;
for ( i = 0; i < nn; i++ )
for ( j = 0; j < nn + 1; j++ )
{
if (j<nn)
{mas[i] [j] = bc2[i] [j];}
else
{mas[i] [j] = cc2[i];}
}
//Сначала все корни по порядку
for ( i = 0; i < nn + 1; i++ )
otv[i] = i;
gauss(mas, nn, otv, x);//вызов метода Гаусса
r=0;
for ( i = 0; i < nn; i++ )
for ( j = 0; j < nn; j++ )
if ( i == otv[j] )
{ //Расставляем корни по порядку
ac2[r]=x[j];
r=r+1;
break;
}
//находим коэффициенты b3 и c3
double Sb3=0, Sc3=0;
for(int p=0;p<4;p++)
{
for(int q=0;q<4;q++)
{
Sc3=0;
Sb3=0;
for(i=0;i<N;i++)
{
Sb3=pow(xi[i],(p+q)) + Sb3;
Sc3=pow(xi[i], p)*yi[i] + Sc3;
}
bc3[p][q]=Sb3;
}
cc3[p]=Sc3;
}
//находим коэффициенты а3 методом Гаусса
nn=4;
for ( i = 0; i < nn; i++ )
for ( j = 0; j < nn + 1; j++ )
{
if (j<nn)
{mas[i] [j] = bc3[i] [j];}
else
{mas[i] [j] = cc3[i];}
}
//Сначала все корни по порядку
for ( i = 0; i < nn + 1; i++ )
otv[i] = i;
gauss(mas, nn, otv, x);//вызов метода Гаусса
r=0;
for ( i = 0; i < nn; i++ )
for ( j = 0; j < nn; j++ )
if ( i == otv[j] )
{ //Расставляем корни по порядку
ac3[r]=x[j];
r=r+1;
break;
}
//
double sume1=0;
double sume2=0;
double sume3=0;
for(i=0;i<N;i++)
{
linfunc[i]=ac[0]+ac[1]*xi[i];
elin[i]=abs(linfunc[i]-yi[i]);
sume1=sume1+elin[i];
kvadfunc[i]=ac2[0]+ac2[1]*xi[i]+ac2[2]*xi[i]*xi[i];
ekvad[i]=abs(kvadfunc[i]-yi[i]);
sume2=sume2+ekvad[i];
kubfunc[i]=ac3[0]+ac3[1]*xi[i]+ac3[2]*xi[i]*xi[i]+ac3[3]*xi[i]*xi[i]*xi[i];
ekub[i]=abs(kubfunc[i]-yi[i]);
sume3=sume3+ekub[i];
}
printf(" X ");
printf(" Y ");
printf(" LIN ");
printf(" KVAD ");
printf(" KUB ");
printf(" e1 ");
printf(" e2 ");
printf(" e3\n");
for(i=0;i<N;i++) // вывод значений
{
printf("%8.2f ", xi[i]);
printf("%8.2f ", yi[i]);
printf("%8.2f ", linfunc[i]);
printf("%8.2f ", kvadfunc[i]);
printf("%8.2f ", kubfunc[i]);
printf("%8.2f ", elin[i]);
printf("%8.2f ", ekvad[i]);
printf("%8.2f ", ekub[i]);
printf("\n");
}
printf("%53.2f ",sume1);
printf("%8.2f",sume2);
printf("%9.2f\n",sume3);
}
void gauss(double mas[][N+1], int nn, int otv[],double x[])
{
//Корни системы
int i, j, k=0;
//Прямой ход метода Гаусса
for ( k = 0; k < nn; k++ )
{ //На какой позиции должен стоять главный элемент
glavelem( k, mas, nn, otv ); //Установка главного элемента
for ( j = nn; j >= k; j-- )
mas[k] [j] /= mas[k] [k];
for ( i = k + 1; i < nn; i++ )
for ( j = nn; j >= k; j-- )
mas[i] [j] -= mas[k] [j] * mas[i] [k];
}
//Обратный ход
for ( i = 0; i < nn; i++ )
x[i] = mas[i] [nn];
for ( i = nn - 2; i >= 0; i-- )
for ( j = i + 1; j < nn; j++ )
x[i] -= x[j] * mas[i] [j];
}
//----------------------------------------------
//Описание функции
//----------------------------------------------
void glavelem( int k, double mas[] [N + 1], int nn, int otv[] )
{
int i, j, i_max = k, j_max = k;
double temp;
//Ищем максимальный по модулю элемент
for ( i = k; i < nn; i++ )
for ( j = k; j < nn; j++ )
if ( fabs( mas[i_max] [j_max] ) < fabs( mas[i] [j] ) )
{
i_max = i;
j_max = j;
}
//Переставляем строки
for ( j = k; j < nn + 1; j++ )
{
temp = mas[k] [j];
mas[k] [j] = mas[i_max] [j];
mas[i_max] [j] = temp;
}
//Переставляем столбцы
for ( i = 0; i < nn; i++ )
{
temp = mas[i] [k];
mas[i] [k] = mas[i] [j_max];
mas[i] [j_max] = temp;
}
//Учитываем изменение порядка корней
i = otv[k];
otv[k] = otv[j_max];
otv[j_max] = i;
Результат работы программы.
Документ
Категория
Рефераты
Просмотров
5
Размер файла
76 Кб
Теги
otchet
1/--страниц
Пожаловаться на содержимое документа