#include<conio.h>
#include "matrix.h"
typedef matrix<double> dmatrix;
typedef vector<double> dvector;
/*Определим вначале некоторые вспомогательные подпрограммы*/
//Для вычисления числа Фибоначчи с заданным номером
long Fib(long n)
{
long F[3]={0,1};
if(n<0)
throw xmsg("Некорректный номер числа");
if(n<2)
return F[n];
for(long i=2;i<=n;i++)
{
F[2]=F[0]+F[1];
F[0]=F[1];
F[1]=F[2];
}
return F[2];
}
/*Простая парабола для тестирования одномерного поиска*/
double gv(double x)
{
return (4*x*x+5*x+1);
}
/*Функциональный класс объектов для поиска экстремума; мы включили в него и некоторые функции определения значения вещественного корня*/
class MonoExtremum
{
double c, d;//Границы для одномерного поиска
double eps;//Для допустимой погрешности поиска
/*Указатель на подпрограмму, возвращающую значение исследуемой функции*/
double (*getv)(double x);
public:
//Конструктор для одномерного поиска
MonoExtremum(double C, double D, double EPS, double (*GETV)(double X)): c(C), d(D), eps(EPS), getv(GETV) { }
//Прототипы функций поиска вещественного корня
double RealRoot(long cod); /*сокращением интервала неопределенности*/
double DNewton(); //Метод секущих
double Stoh(); //Стохастическая аппроксимация
/*Прототипы функций решения задачи условной оптимизации*/
double Dihot(); //метод дихотомии
double Fibon(long cnt); //метод Фибоначчи
double Gold(); //метод золотого сечения
dvector Coord(); //метод покоординатного спуска
dvector Grad(); //метод градиента
dvector PSM(); //последовательный симплекс-метод
};
/*Обобщенная функция уточнения корня последовательным сокращением интервала неопределенности*/
double MonoExtremum::RealRoot(long cod=0)
{
double ct=c,dt=d,x,vc,vd,vx,k;
for(;;)
{
switch(cod)
{
case 0: //Больцано
k=0.5;
vc=getv(ct);
break;
case 1: //Хорды
vc=getv(ct);
vd=getv(dt);
k=vc/(vd+vc);
if(k<0)
k=-k;
break;
case 2: //Случайный поиск
randomize();
k=double(1+random(1111))/1111.0;
vc=getv(ct);
break;
}
x=ct+k*(dt-ct);
vx=getv(x);
if(fabs(vx)<eps)
return x;
if(vx*vc>0)
ct=x;
else
dt=x;
}
}
//Метод секущих - дискретный аналог метода Ньютона
double MonoExtremum::Dnewton()
{
double sx[3];
sx[0]=(c+d)/2;
sx[1]=sx[0]+1e-4;
sx[2]=sx[1]-(getv(sx[1]))*(sx[1]-sx[0])/
(getv(sx[1])-getv(sx[0]));
while(1)
{
if(getv(sx[2])<eps)
return sx[2];
sx[0]=sx[2];
sx[1]=sx[0]+1e-4;
sx[2]=sx[1]-(getv(sx[1]))*(sx[1]-sx[0])/
(getv(sx[1])-getv(sx[0]));
}
}
//Стохастическая аппроксимация
double MonoExtremum::Stoh()
{
double scnt=1;
//Для определения знака приращения
double sx[3],fa[3];
sx[0]=(c+d)/2.0;
fa[0]=fabs(getv(sx[0]));
if(fa[0]<eps)
return sx[0];
sx[1]=sx[0]+1e-2;
fa[1]=fabs(getv(sx[1]));
if(fa[1]<eps)
return sx[1];
for(scnt=1;;scnt++)
{
sx[2]=sx[1]-(20.0/scnt)*fa[1]*(sx[1]-sx[0])/ (fa[1]-fa[0]);
fa[2]=fabs(getv(sx[2]));
if(fa[2]<eps)
return sx[2];
sx[0]=sx[1];
sx[1]=sx[2];
fa[0]=fa[1];
fa[1]=fa[2];
}
}
//Реализация функций одномерного поиска экстремума
//Метод дихотомии
double MonoExtremum::Dihot()
{
double xl,xr,yl,yr,dt=d,ct=c;
for(;(dt-ct)>eps;)
{
/*задаем координаты 2-х точек эксперимента вблизи средины интервала*/
xl=(ct+dt)/2-0.01*(dt-ct);
xr=(ct+dt)/2+0.01*(dt-ct);
//Вычисляем значения функции в этих точках
yl=getv(xl);
yr=getv(xr);
if(sign(yl-yr)>0)
/*Бывшая правая - теперь левая, а правая - симметрично*/
{
ct=xl;
xl=xr;
xr=ct+dt-xl;
}
else
{
dt=xr;
xr=xl;
xl=ct+dt-xr;
}
}
return (ct+dt)/2;
}
//Метод Фибоначчи для поиска экстремума
double MonoExtremum::Fibon(long cnt)
{
double k=(double)Fib(cnt-1)/(double)Fib(cnt);
/*вычисляем координаты 2-х первых точек эксперимента*/
double ct=c, dt=d, xl=d-k*(d-c), xr=c+k*(d-c), yl, yr;
/*Далее в цикле в худшую точку переносим ближайшую границу, а добавочную точку определяем симметрично оставшейся относительно середины интервала */
for(long i=cnt;i>1;i--)
{
yl=getv(xl);
yr=getv(xr);
if(sign(yl-yr)>0)
{
ct=xl;
xl=xr;
xr=ct+dt-xl;
}
else
{
dt=xr;
xr=xl;
xl=ct+dt-xr;
}
}
return (ct+dt)/2;
}
//Метод золотого сечения
double MonoExtremum::Gold()
{
double k=0.618033989;
/*вычисляем координаты 2-х первых точек эксперимента*/
double ct=c, dt=d, xl=d-k*(d-c), xr=c+k*(d-c), yl, yr;
/*Далее в цикле в худшую точку переносим ближайшую границу, а добавочную точку определяем симметрично оставшейся относительно середины интервала*/
for(;(dt-ct)>eps;)
{
yl=getv(xl);
yr=getv(xr);
if(sign(yl-yr)>0)
{
ct=xl;
xl=xr;
xr=ct+dt-xl;
}
else
{
dt=xr;
xr=xl;
xl=ct+dt-xr;
}
}
return (ct+dt)/2;
}
//Программа тестирования методов уточнения корней и
//одномерного поиска экстремума
main()
{
/*Тестирование методов уточнения корня последовательным сужением интервала неопределенности*/
MonoExtremum root1(-0.5,0.5,1e-3,gv); //Сконструируем объект
double res0=root1.RealRoot(0); //0-метод Больцано
cout<<endl<<"Bolzano: "<<res0<<" Func : "<< gv(res0);
double res1=root1.RealRoot(1); //1-метод хорд
cout<<endl<<"Chord : "<<res1<<" Func : "<< gv(res1);
double res2=root1.RealRoot(2); /*2-метод Монте _ Карло*/
cout<<endl<<"Rand :"<<res2<<" Func : "<<gv(res2);
//Тестирование поиска корня методом секущих
double resDNewton=root1.DNewton();
cout<<endl<<"DNewton : "<<resDNewton<<" Func : "
<<gv(resDNewton);
/*Тестирование поиска корня методом стохастической аппроксимации*/
double resStoh=root1.Stoh();
cout<<endl<<"Stoh : "<<resStoh<<" Func : "
<<gv(resStoh);
//Тестирование методов поиска экстремума
MonoExtremum extr1(-3.0,+5.0,1e-3,gv);/*Сконструируем объект*/
//Вызываем методы решения
double resDihot=extr1.Dihot(), //Дихотомии
resFib=extr1.Fibon(25), //Фибоначчи
resGold=extr1.Gold(); //Золотого сечения
cprintf("\r\nDihot : %lf Fibon : %lf Gold : %lf ", resDihot,resFib,resGold);
getch();
}
