Текст программы реализующей аппроксимацию методом наименьших квадратов

 

// min_square.cpp : main project file.

 

#include "stdafx.h"

#include<iostream>

#include<conio.h>

#include<cmath>

#include<iomanip>

#define size 20

 

using namespace System;

using namespace std;

 

//объявление прототипов функций

void matrix (double[][size],int,double[]); /*вычисляет значение матрицы методом Холесского, включает в себя две функции: holess, и res*/

void holess(double[][size],int);

void res(double[][size], int,double[]);

//функция построения матрицы Грама

void gram(double[][size],double[][size],int, double[][size]);

 

//начало основной программы

int main()

{ 

  /*объявленные переменные: в массивы "х" и "у" накапливаются табличные значения аргументов и функции. В массив "с" сохраняются значения коэффициентов (результат решения матрицы Грама). Массив "а" хранит матрицу Грама, "у1"-значения функции, полученные путем вычисления аппроксимирующего многочлена. s-хранит промежуточные результаты вычисления функции апроксимирующим многочленом. р-невязка.n-количество табличных значений аргументов*/

      

double x[size][size],y[size][size],c[size],a[size][size],y1[size],s,p;

int n;

 

cout<<"Enter quantity of knots of interpolation N="; //ввод n

cin>>n;

cout<<endl;

  

for (int i=1;i<=n;i++)          //начало цикла ввода табличных значений

{

  cout<<"Inter x("<<i<<")=";

  cin>>x[i][1];      /*табличные значения аргумента сохраняются в первом           

столбце многомерного массива "х"*/

 

  for(int j=0;j<=(n-1);j++)

x[i][j]=pow(x[i][1],j); /*табличное значение "х" первого столбца i-й

строки возводится в j-ю степень и сохраняется в j-м столбце*/

      

  cout<<"Inter y("<<i<<")=";

  cin>>y[i][0];      /*сохранение табличных значений функции в нулевом

 столбце многомерного массива "у" */

 

  for(int j=0;j<=(n-3);j++)

     y[i][j]=y[i][0]*x[i][j]; /*в j-й столбец i-й строки массива "у" записывается произведение значения 0-го столбца массива "у" (табличное значение у) и х^j */

          

  cout<<endl;

}

 

gram(x,y,n,a);                         //заполнение матрицы Грама

matrix(a,n-2,c);                       //вычисление матрицы методом Холесского

  

 

 

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

{

  s=0;

  for(int j=1;j<=n-2;j++)

        s=s+c[j]*x[i][j-1]; //вычисление функции аппроксимирующим многочленом

 

  y1[i]=s; //запись вычисленного значения функции в массив выходных данных

}

s=0;

cout<<setw(2)<<"X"<<setw(12)<<"Y"<<setw(19)<<"F"<<endl;

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

{

  s=s+pow((y1[i]-y[i][0]),2);               //нахождение разницы между табличным и вычисленным значением функции

        cout<<setw(4)<<x[i][1]<<setw(14)<<y[i][0]<<setw(19)<<y1[i]<<endl; //вывод табличных и вычисленных значений

}

      

p=sqrt(s/(n+1));                       //вычисление невязки

cout<<endl;

cout<<setw(2)<<"P="<<p<<endl;          //вывод невязки на экран

   

getch();

return 0;                             //конец функции main

}

 

 

/*начало функции, строящей матрицу Грама. Входные данные: двумерный массив табличных, а также возведенных в степень n-1 значений аргументов "х"; двумерный массив табличных значений, а так же результатов произведений аргумента и функции "у"; количество узловых точек. Выходные данные двумерный массив "а" - матрица Грама*/

 

 

void gram(double x [][size],double y[][size],int n, double a[][size])

{

  int k;

            

  for(int j=0;j<=(n-1);j++)

  {

x[n+1][j]=y[n+1][j]=0;

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

x[n+1][j]=x[n+1][j]+x[i][j]; //вычисление значений членов левой части

матрицы Грама

            

        cout<<endl;

  }

      

  for(int j=0;j<=n-3;j++)                

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

y[n+1][j]=y[n+1][j]+y[i][j]; //вычисление значений членов правой части матрицы Грама

 

  k=0;

 

  for(int i=1;i<=n-2;i++)

  {

        for(int j=1;j<=(n-2);j++)

a[i][j]=x[n+1][j-1+k]; //заполнение массива "а" значениями аргументов

a[i][n-1]=y[n+1][i-1]; //заполнение столбца правых частей системы

                   нормальных уравнений

        k++;

  }

}

 

 

void matrix (double a[][size],int n,double c[])

{

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

        c[i]=0;

  holess(a,n);

res(a,n,c);

}

 

//функция реализует метод Холесского с частичной перестановкой для выбора главного элемента

//Входные данные массив коэффициентов "аrg", количество уравнений системы n

//Выходные данные массив новых коэффициентов "аrg".

 

void holess(double arg[][size],int n)

{

double max,p,s,a1[size];

int i1;

 

for(int m=1;m<=n;m++)

 {

  if(m==1)

        for(int j=2;j<=(n+1);j++)

               arg[1][j]=arg[1][j]/arg[1][1];

 

  if(m>=2)

  {

        max=fabs(arg[m][m]);

        for(int t=m;t<=n;t++)

               if(fabs(arg[t][m])>=max) //выбор главного элемента

               {

                      max=fabs(arg[t][m]);

                             i1=t;

               }

 

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

               a1[i]=arg[m][i]; //a1(n+1) - массив коэффициентов

            

        for(int i=1;i<=(n+1);i++) //замена ведущей строки на строку с выбранным главным элементом

        {

               arg[m][i]=arg[i1][i];

               arg[i1][i]=a1[i];

        }

 

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

        {

               p=0;

               for(int k=1;k<=(m-1);k++)

                      p=p+arg[i][k]*arg[k][m];

 

               arg[i][m]=arg[i][m]-p; //формирование очередного a(i,m) столбца

        }

 

        for(int j=(m+1);j<=(n+1);j++)

        {

               s=0;

               for(int k=1;k<=(m-1);k++)

                      s=s+arg[m][k]*arg[k][j];

               arg[m][j]=(arg[m][j]-s)/arg[m][m]; //формирование очередной a(m,j)строки

        }

  }

 }

}

 

 

//функция формирования решения системы уравнений

//входные данные: двумерный массив коэффициентов "arg", n-число уравнений

//выходные данные:массив решений системы х(i)

 

void res(double arg[][size], int n,double x[])

{

  double s;

  for(int i=n;i>=1;i--)

  {

        if(n==1)

               x[i]=arg[n][n+1];

            

        s=0;

 

        for(int k=i+1;k<=n;k++)

               s=s+arg[i][k]*x[k];

 

        x[i]=arg[i][n+1]-s;

  } //конец цикла

 

}       //конец функции res

 

Выводы:

 

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

 

 

Заключение:

 

В настоящей работе были рассмотрены три способа аппроксимации одной и той же функции. У всех трех способов есть своя область применения, достоинства и недостатки. Применимо к исследуемой функции, учитывая равный интервал между значениями аргументов, малый шаг, а так же то, что функция близка к линейной, наиболее предпочтительным для ее аппроксимации является метод Ньютона. Т.к. уже линейная интерполяция дала результат с высокой точностью.

 

 

Список литературы:

 

1. Гловацкая А.П. Методы и алгоритмы вычислительной математики. Уч.пос. для ВУЗОВ.-М. Радио и связь, 1999.-408с.

2. Гловацкая А.П. Сборник задач для курсовой работы по курсу Информатика. МТУСИ. 2006, 32с.

3. Дейтел Х.М. Как программировать на С++: Пятое издание. - С.Пб. ООО «Бином-Пресс», 2011.-1456с.

4. Пахомов Б.И. С/С++ и MS Visual C++ 2010 для начинающих – СПб.: БХВ - Петербург, 2011. – 736 с.