Построить гладкую кривую, проходящую приблизительно через заданные точки. Точки расставляются на экране с помощью мыши в произвольном порядке. Программа также позволяет с помощью мыши изменять форму кривой.
Коэффициенты аппроксимируещего многочлена вычисляются при помощи метода наименьших квадратов. Эта задача решается следующим образом.
Пусть на некотором отрезке в точках 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,®); //скрыть курсор мыши перед прорисовкой
circle(x[i],y[i],2);
reg.r_ax=1; intr(0x33,®); //снова показать курсор мыши
}
setcolor(3);
reg.r_ax=2; intr(0x33,®);
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,®);
}
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.r_ax=1; intr(0x33,®); //инициализируем мышь и отображаем курсор
while(!kbhit()) //пока не нажата клавиша на клавиатуре, выполняем цикл
{
reg.r_ax=3; intr(0x33,®);//прерывание, необходимое для определения состояния мыши
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,®);
x[correct]=reg.r_cx;
y[correct]=reg.r_dx; //перезаписываем в массиве координаты
//перемещеной точки
reg.r_ax=2; intr(0x33,®);
Calculate();
if((it++)%2)DrawCurve(); //перерисовываем кривую в новом положении
reg.r_ax=1; intr(0x33,®);
}
}
//если значение переменной correct осталось равным -1, то ставим новую узловую точку
else if(pts<14)
{
x[pts]=reg.r_cx;
y[pts]=reg.r_dx; //добавляем в массив координат кординаты новой точки
reg.r_ax=2; intr(0x33,®);
circle(x[pts],y[pts],3); //кружочек в новой точке
reg.r_ax=1; intr(0x33,®);
delay(100); pts++; //инкрементируем счетчик узловых точек
if(pts>1)
{
K=pts-1; // количество коэффициентов полинома должно быть меньше,
// чем количество узловых точек
Calculate();
DrawCurve(); //вычисляем коэффициенты и прорисовываем кривую
}
}
correct=-1; //необходимо для следующей итерации цикла
}
}
}
Ключевые слова:
аппроксимация методом наименьших квадратов, интерполяция, гладкая кривая
