Метод наименьших квадратов

Размещено на http://www.allbest.ru/

Министерство транспорта Российской Федерации

Федеральное государственное автономное

Образовательное учреждение высшего образования

"Российский университет транспорта (МИИТ)"

РУТ (МИИТ)

Кафедра ПМ-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. Выводы:

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

Размещено на Allbest.ru