Аппроксимация функций полиномом методом наименьших квадратов

square.jpg

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

Коэффициенты аппроксимируещего многочлена вычисляются при помощи метода наименьших квадратов. Эта задача решается следующим образом.
Пусть на некотором отрезке в точках x0, x1, x2… xn нам известны значения y0, y1, y2… yn некоторой функции f(x).

Требуется определить параметры ai многочлена вида F(x)=a0+a1x+a2x2+...+akxk, где k<N
такого, что сумма квадратов отклонений значений y от значений функции f(y) в заданных точках x была минимальной, т.е.
S=Σ(yi-f(xi))2->min
Геометрический смысл заключается в том, что график найденного многочлена y = f(x) будет проходить как можно ближе к каждой из заданных точек.

Такая задача будет решена, если решить систему уравнений вида:
a0n +a1Σxi+a2Σxi2+... +akΣxik=Σyi
a0Σxi +a1Σxi2+a2Σxi3+... +akΣxik+1=Σxiyi
...
a0Σxik+a1Σxik+1+a2Σxik+2+...+akΣxi2k=Σxikyi

где везде под символом Σ подразумевается суммирование по i от 0 до n.

Для решения системы воспользуемся методом Гаусса (система уравнений приводится к треугольному виду, после чего последовательно вычисляются аn ,… а0). Имея аналитический вид функции f, нетрудно найти её значение в любой точке и построить кривую.

#include<conio.h>
#include<math.h>
#include<graphics.h>
#include<dos.h>
 
#define  N 14 		//максимальное количество точек;					               
 
#define GRAPHICS_DIR "D:\\BORLAND\\BGI"
 
int K,pts=0;		//К-степень аппроксимирующего полинома, pts-текущее количество узловых точек
//a- массив коэффициентов полинома; b- массив свободных членов системы уравнений; x,y- координаты узловых точек;
long double a[N],b[N],x[N],y[N],sums[N][N];
struct REGPACK reg;
 
long double Power(long double v,int p);	//возводит v в степень p;
void Calculate();				//вычисляет коэффициенты полинома 
void Refresh();				//обнуляет все массивы перед повторным вычислением коэффициентов
void DrawCurve();			//прорисовывает график полинома
 
long double Power(long double v,int p)
{
        if(p==0) return 1;
        if(p>1) v*=Power(v,p-1);
        return v;
}
 
void Calculate()
{
        int i,j,k;
        long double s,t,M;
        Refresh();
	//упорядочиваем узловые точки по возрастанию абсцисс
	for(i=0; i<pts; i++)		
   	{
       		for(int j=i;j>=1;j--)
       			if(x[j]<x[j-1])
       			{
                		t=x[j-1]; x[j-1]=x[j]; x[j]=t;
                		t=y[j-1]; y[j-1]=y[j]; y[j]=t;
       			}
   	}
	//заполняем коэффициенты системы уравнений
   	for(i=0; i<K+1; i++) 
	{
       		for(j=0; j<K+1; j++)
		{
	   		sums[i][j] = 0;
	   		for(k=0; k<pts; k++)
	       			sums[i][j]+=Power(x[k],i+j);
	       	}
   	}
	//заполняем столбец свободных членов
    	for(i=0; i<K+1; i++)		
       	{
                b[i]=0;
	 	for(k=0; k<pts; k++)
			b[i] +=Power(x[k],i)*y[k];
        }
	//применяем метод Гаусса для приведения матрицы системы к треугольному виду
        for(k=0; k<K+1; k++)
        {
                for(i=k+1; i<K+1; i++)
                {
	                M=sums[i][k]/sums[k][k];
	                for(j=k; j<K+1; j++)
                                 sums[i][j] -= M * sums[k][j];
           	        b[i] -= M*b[k];
                }
        }
	//вычисляем коэффициенты аппроксимирующего полинома
        for(i=K;i>=0;i--)
        {
                s=0;
                for(j = i; j<K+1; j++)
	              s+=sums[i][j]*a[j];
                a[i] = (b[i]-s)/sums[i][i];
        }
}
 
void DrawCurve()
{
 
        long double X,Y,Yp;  //Yp- у точки, из которой проводим линию в следующую, только что вычисленную
        cleardevice();	     //очишаем экран
        setcolor(WHITE);
        //в узловых точках рисуем кружочки
         for(int i=0;i<pts;i++)
        {
                reg.r_ax=2; intr(0x33,&reg);	//скрыть курсор мыши перед прорисовкой
                circle(x[i],y[i],2);
                reg.r_ax=1; intr(0x33,&reg);	//снова показать курсор мыши
        }
 
        setcolor(3);
        reg.r_ax=2; intr(0x33,&reg);
        for(X=0;X<640;X+=0.1)			//Х пробегает весь экран
        {
                Y=0;
                for(int j=0;j<=K;j++)
                        Y+=a[j]*Power(X,j); //высчитываем значение полинома, заданного своими коэффициентами, в точке Х
                if(X==0) Yp=Y;		    //при первой итерации ордината предыдущей точки равна ординате текущей
                line(X-0.1,Yp,X,Y);	    //соединяем линией предыдущую и текущую точки
                Yp=Y;			    //после отрисовки текущая точка становится предыдущей
        }
        reg.r_ax=1; intr(0x33,&reg);
}
 
void Refresh()
{
        for(int i=0;i<K;i++)
        {
                a[i]=b[i]=0;
                for(int j=0;j<K;j++)
                        sums[i][j]=0;
        }
}
 
void main()
{
        int GD,GM;
	int correct=-1;				//показывает, какая точка в данный момент переиещается мышью
        GD=DETECT;
        initgraph(&GD,&GM,GRAPHICS_DIR);	//инициализируем графику
        reg.r_ax=0; intr(0x33,&reg);
	reg.r_ax=1; intr(0x33,&reg);		//инициализируем мышь и отображаем курсор
        while(!kbhit())				//пока не нажата клавиша на клавиатуре, выполняем цикл
        {
                reg.r_ax=3; intr(0x33,&reg);//прерывание, необходимое для определения состояния  мыши
                 if(reg.r_bx==1)			//если левая кнопка нажата
		{
         //если левая кнопка по прежнему нажата, определяем, какая узловая точка перетаскивается
			if (reg.r_bx==1)
                                for(int i=0;i<pts;i++)
                                        if((abs(x[i]-reg.r_cx)<3)&&(abs(y[i]-reg.r_dx)<3))
		                        	correct=i;
			//если положение какаой-то точки откорректировано, то...
                        if (correct!=-1)
                        {
                                long it=0;
                                while(reg.r_bx==1)  //...пока левая кнопка нажата
                                {
                                        reg.r_ax=3; intr(0x33,&reg);
                                        x[correct]=reg.r_cx;		
                                        y[correct]=reg.r_dx; //перезаписываем в массиве координаты 
                                                             //перемещеной точки
 
                                        reg.r_ax=2; intr(0x33,&reg);
                                        Calculate();
                                        if((it++)%2)DrawCurve(); //перерисовываем кривую в новом положении
                                        reg.r_ax=1; intr(0x33,&reg); 
 
                                }
                        }
//если значение переменной correct осталось равным -1, то ставим новую узловую точку
                        else if(pts<14)
                        {
                                x[pts]=reg.r_cx;
                                y[pts]=reg.r_dx;   //добавляем в массив координат кординаты новой точки
                                reg.r_ax=2; intr(0x33,&reg);
		 	        circle(x[pts],y[pts],3);  //кружочек в новой точке
                                reg.r_ax=1; intr(0x33,&reg);
 
                                delay(100);    pts++;	//инкрементируем счетчик узловых точек
 
                                if(pts>1)
   			        {
				        K=pts-1; // количество коэффициентов полинома должно быть меньше, 
                                                 // чем количество узловых точек
                                        Calculate();			
                                        DrawCurve(); //вычисляем коэффициенты и прорисовываем кривую
                                }
                        }
                        correct=-1;  //необходимо для следующей итерации цикла
                }
        }
}

Ключевые слова: 
аппроксимация методом наименьших квадратов, интерполяция, гладкая кривая
ВложениеРазмер
Аппроксимация1.27 кб