Министерство транспорта Российской Федерации
Федеральное государственное автономное
Образовательное учреждение высшего образования
"Российский университет транспорта (МИИТ)"
РУТ (МИИТ)
Кафедра ПМ-1
Лабораторная работа
Метод наименьших квадратов
Выполнила: Илюшина Е.Р.
Проверил: Посвянский В.П.
Москва 2019
Постановка задачи
1. Составить таблицу значений функции y = на отрезкес шагом . Полученную таблицу значений функции рассматривать в дальнейшем в качестве исходных данных для последующей среднеквадратичной аппроксимации.
2. Составить системы уравнений для определения коэффициентов многочленов наилучшего среднеквадратичного приближения (м.н.с.п.) и решить их с точностью не хуже 0.0001. Под точностью решения СЛАУ понимается наибольшее значение модуля разности между значением левой части уравнений после подстановки полученного решения в систему и соответствующей правой частью уравнений.
Найти эти многочлены порядка не более 1, 2, 5, 6, напечатать их коэффициенты и достигнутую точность.
3. Составить таблицу значений многочленов наилучшего среднеквадратичного приближения не более 1, 2, 5, 6 порядка на отрезке [-1.5,1.5] с шагом h1 = 3/50. Протабулировать заданную функцию на том же отрезке с шагом h1. Напечатать таблицу.
4. Построить графики аппроксимируемой функции, заданной на дискретном множестве точек п.1. и интерполяционных многочленов, используя полученные значения в таблицах п.3. Вывести на печать таблицы п.3.
Теория
Многочлен наилучшего среднеквадратичного приближения (м.н.с.п.):
Система линейных алгебраических уравнений для определния коэффициентов м.н.с.п.:
Таблица узлов:
|
xi |
-1,5 |
-1,2 |
-0,9 |
-0,6 |
-0,3 |
0 |
0,3 |
0,6 |
0,9 |
1,2 |
1,5 |
|
|
yi |
4,813 |
2,634 |
1,846 |
1,77 |
1,918 |
2 |
1,918 |
1,77 |
1,846 |
2,634 |
4,813 |
Коэффициенты м.н.с.п. ,
|
n=1 |
2,542 |
1,557е-05 |
||||||
|
n=2 |
1,417 |
1,684е-05 |
1,25 |
|||||
|
n=5 |
2 |
-5,209е-07 |
-0,9995 |
4,576е-05 |
0,9998 |
-2,075е-05 |
||
|
n=6 |
1,999 |
-5,209е-07 |
-0,9955 |
4,576е-05 |
0,9949 |
-2,075е-05 |
0,001442 |
Код программы:
#include "stdafx.h"
#include <iostream>
#include <conio.h>
#include <iomanip>
#include <math.h>
#include <fstream>
using namespace std;
ofstream fout;
ifstream fin;
int acc = 0.0001;
class coor
{
double x[51];
double y[51];
int len;
public:
coor() {
len = 0;
}
coor(double a, double b, int divider) {
input(a, b, divider);
}
void input(double a, double b, int divider) {
int i;
double h = (abs(a) + abs(b)) / divider;
len = divider + 1;
for (i = 0; i < len; i++) {
x[i] = a + i * h;
y[i] = x[i] * pow(2.71828183, (-1) * x[i]);
}
return ;
}
void output() {
int i;
fout " setw(10) " "x";
for(i = 0; i < len; i++) {
fout " setw(10) " x[i];
}
fout " "\n";
fout " setw(10) " "y";
for(i = 0; i < len; i++) {
fout " setw(10) " y[i];
}
fout " "\n";
return ;
}
void set(double u, double v) {
x[len] = u;
y[len] = v;
len++;
return ;
}
double get_x(int number) {
return x[number];
}
double get_y(int number) {
return y[number];
} уравнение коэффициент многочлен
};
class tabl
{
coor real_values[2];
coor MNK[4];
double **a, **x, *r;
public:
tabl() {
int i, j;
x = new double *[4];
for(j = 0; j < 4; j++) {
x[j] = new double [7];
}
a = new double *[7];
for (i = 0; i < 7; i++) {
a[i] = new double[7];
}
r = new double[7];
}
void input_real_values(double a, double b) {
real_values[0].input(a, b, 10);
real_values[1].input(a, b, 50);
return ;
}
void find() {
int i, j, k;
double sum_a = 0, sum_r = 0;
for (i = 0; i < 7; i++) {
for (j = 0; j < 7; j++) {
for (k = 0; k < 11; k++) {
sum_a = sum_a + pow(real_values[0].get_x(k), i + j);
}
a[i][j] = sum_a;
sum_a = 0;
}
for (k = 0; k < 11; k++) {
sum_r = sum_r + pow(real_values[0].get_x(k), i) * real_values[0].get_y(k);
}
r[i] = sum_r;
sum_r = 0;
}
a[0][0] = 11;
return ;
}
void check(double *A[], double X[], double R[], int n)
{
int i, j, k;
double sum, max = 0;
for (i = 0; i < n; i++) {
sum = 0;
for (j = 0; j < n; j++) {
sum = sum + A[i][j] * X[j];
}
if (abs(sum - R[i]) > max) {
max = abs(sum - R[i]);
}
cout " "Проверка выполнена. Невязка равна " " max " "\n";
}
return ;
}
void gauss_method_main_element(int acc, int n, int number)
{
double **A, *X, *R;
int i, j, k;
int *z;
double sum;
z = new int[n];
R = new double [n];
X = new double [n];
A = new double *[n];
for (i = 0; i < n; i++) {
A[i] = new double[n];
}
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
A[i][j] = floor(a[i][j] * pow(10.0, acc)) / pow(10.0, acc);
z[j] = j;
}
R[i] = floor(r[i] * pow(10.0, acc)) / pow(10.0, acc);
}
double c, s;
int f, l, p, q, d;
for (i = 0; i < n - 1; i++) {
s = 0, f = i, l = i;
for (p = i; p < n; p++) {
for (q = i; q < n; q++) {
if (abs(A[p][q]) > s) {
s = abs(A[p][q]);
f = p;
l = q;
}
}
}
for(q = 0; q < n; q++) {
c = A[i][q];
A[i][q] = A[f][q];
A[f][q] = c;
}
s = R[i];
R[i] = R[f];
R[f] = s;
for (p = 0; p < n; p++) {
c = A[p][i];
A[p][i] = A[p][l];
A[p][l] = c;
}
d = z[i];
z[i] = z[l];
z[l] = d;
for (j = i + 1; j < n; j++) {
A[i][j] = A[i][j] / A[i][i];
for (k = i + 1; k < n; k++) {
A[k][j] = A[k][j] - A[k][i] * A[i][j];
}
}
R[i] = R[i] / A[i][i];
for (k = i + 1; k < n; k++) {
R[k] = R[k] - R[i] * A[k][i];
}
}
if (A[i][i] != 0) {
R[i] = R[i] / A[i][i];
}
for (i = n - 1; i >= 0; i--) {
sum = 0;
for (k = i + 1; k < n; k++) {
sum = sum + A[i][k] * X[k];
}
X[i] = R[i] - sum;
}
double *help = new double[n];
for(i = 0; i < n; i++) {
help[z[i]] = X[i];
}
for(i = 0; i < n; i++) {
X[i] = help[i];
}
check(a, X, r, n);
for (i = 0; i < n; i++) {
x[number][i] = X[i];
}
for (i = 0; i < n; i++) {
delete [] A[i];
}
delete [] A;
delete [] X;
delete [] R;
delete [] z;
return ;
}
void M_N_K(int n[]) {
int i, j;
int number;
double sum = 0;
for (number = 0; number < 4; number++) {
for (j = 0; j < 51; j++) {
for (i = 0; i < n[number]; i++) {
sum = sum + x[number][i] * pow(real_values[1].get_x(j), i);
}
MNK[number].set(real_values[1].get_x(j), sum);
sum = 0;
}
}
return ;
}
void output(int n[]) {
int i, j;
fout.open("Tabl.txt");
fout " setprecision(4);
fout " "Таблица узлов интерполирования:\n";
real_values[0].output();
fout " "\n";
fout " "Коэффициенты м.н.с.п.:\n";
for (i = 0; i < 4; i++) {
for (j = 0; j < n[i]; j++) {
fout " x[i][j] " ' ';
}
fout " "\n";
}
fout " "Таблица значений функции и м.н.с.п.\n";
fout " setw(10) " "x" " setw(10) " "y" " setw(10) " "P1" " setw(10) " "P2" " setw(10) " "P5" " setw(10) " "P6\n";
for (i = 0; i < 51; i++) {
fout " setw(10) " real_values[1].get_x(i) " setw(10) " real_values[1].get_y(i);
for (j = 0; j < 4; j++) {
fout " setw(10) " MNK[j].get_y(i);
}
fout " "\n";
}
fout.close();
return ;
}
};
int _tmain(int argc, _TCHAR* argv[])
{
int i, j, k, n[] = {2, 3, 6, 7};
double sum;
double a, b;
tabl worker;
setlocale(LC_ALL, "rus");
fin.open("inf.txt");
fin " a " b;
fin.close();
worker.input_real_values(a, b);
worker.find();
for (i = 0; i < 4; i++) {
worker.gauss_method_main_element(4, n[i], i);
}
worker.M_N_K(n);
worker.output(n);
system("pause");
return 0;
}
1. Распечатка результатов:
2. Графики:
3. Выводы:
Как видно из результатов, метод наименьших квадратов позволяет аппроксмировать функцию. С увеличением порядка многочлена наилучшего среднеквадратичного приближения увеличивается точность данного метода.