Решение задачи с использованием метода покоординатного спуска

Описание метода покоординатного спуска

 

Изложим этот метод на примере функции трех переменных .

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

Затем из новой точки сделаем спуск по направлению, параллельному оси , т. е. рассмотрим , найдем ее минимум и обозначим его через . Второй шаг приводит нас в точку . Из этой точки делаем третий шаг – спуск параллельно оси  и находим минимум функции . Приход в точку  завершает цикл спусков.

Будем повторять циклы. На каждом спуске функция не возрастает, и при этом значения функции ограничены снизу ее значением в минимуме . Следовательно, итерации сходятся к некоторому пределу . Будет ли здесь иметь место равенство, т. е. сойдутся ли спуски к минимуму и как быстро?

Это зависит от функции и выбора нулевого приближения.

Метод спуска по координатам несложен и легко программируется на ЭВМ. Но сходится он медленно. Поэтому его используют в качестве первой попытки при нахождении минимума.

Алгоритм решения

 

Будем перебирать с с шагом какой-либо малой длины.

Зададим нулевое приближение, например:

 

 

Найдем частные производные .

 

 

Пусть при каком-то j

Тогда k-ое приближение считаем по формулам:

 

 

Шаг t будем выбирать таким образом, чтобы

 

 и .

В результате, закончив процесс по критерию , где -заданная точность мы придем к набору , при котором функция f максимальна.

Подставим найденный набор  и соответствующее  в функцию f1= и перебрав все с, выберем те , при которых f1 минимальна.

Заключение

 

В ходе решения поставленной задачи рассмотрены случаи, когда n=4,5,6. Были найдены две основные области наборов :

1) x_i=0 или 1(количество 0 и 1 одинаково)

2) x_i=0.5, .

Причем участок 1<p<1.5 покрыт первой областью, при q>>p  –– из первой области, при q≈p  –– из второй области, а при p→∞ область делилась между ними примерно пополам.

На участке p>2 появлялись области вида:

i<k => x_i=0;

i>n-k => x_i=1;

k-1<i<n-k+1=> x_i=0.5.

На участке 2<q<3 p 2 существует область, в которой максимум достигается при  вида:

x_i=a => x_n-i=1-a, 0<a<1.

Список использованных источников

 

1. Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. М.: Высшая школа, 1994. 543с.

2. Березин И.С. и Жидков Н. П. Методы вычислений. т.1. М.: “Наука”, 1965. 633c.

3. Подбельский В.В. и Фомин С.С. Программирование на языке Си. М.: “Финансы и статистика”, 2000. 599с.

Приложение 1. Листинг программы №1

 

//вывод на экран областей максимума функции

#include "stdafx.h"

#include "KE2.h"

#include "math.h"

#include "KE2Doc.h"

#include "KE2View.h"

 

#ifdef _DEBUG

#define new DEBUG_NEW

#undef THIS_FILE

static char THIS_FILE[] = __FILE__;

#endif

 

IMPLEMENT_DYNCREATE(CKE2View, CView)

 

BEGIN_MESSAGE_MAP(CKE2View, CView)

    //{{AFX_MSG_MAP(CKE2View)

              // NOTE - the ClassWizard will add and remove mapping macros here.

              // DO NOT EDIT what you see in these blocks of generated code!

    //}}AFX_MSG_MAP

    // Standard printing commands

    ON_COMMAND(ID_FILE_PRINT, CView::OnFilePrint)

    ON_COMMAND(ID_FILE_PRINT_DIRECT, CView::OnFilePrint)

    ON_COMMAND(ID_FILE_PRINT_PREVIEW, CView::OnFilePrintPreview)

END_MESSAGE_MAP()

CKE2View::CKE2View()

{

}

 

CKE2View::~CKE2View()

{

}

 

BOOL CKE2View::PreCreateWindow(CREATESTRUCT& cs)

{

 

    return CView::PreCreateWindow(cs);

}

 

void CKE2View::OnDraw(CDC* pDC)

{

    CKE2Doc* pDoc = GetDocument();

    ASSERT_VALID(pDoc);

    drawPoint(pDC);

    // TODO: add draw code for native data here

}

 

BOOL CKE2View::OnPreparePrinting(CPrintInfo* pInfo)

{

    // default preparation

    return DoPreparePrinting(pInfo);

}

 

void CKE2View::OnBeginPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)

{

    // TODO: add extra initialization before printing

}

 

void CKE2View::OnEndPrinting(CDC* /*pDC*/, CPrintInfo* /*pInfo*/)

{

    // TODO: add cleanup after printing

}

 

#ifdef _DEBUG

void CKE2View::AssertValid() const

{

    CView::AssertValid();

}

 

void CKE2View::Dump(CDumpContext& dc) const

{

    CView::Dump(dc);

}

 

CKE2Doc* CKE2View::GetDocument() // non-debug version is inline

{

    ASSERT(m_pDocument->IsKindOf(RUNTIME_CLASS(CKE2Doc)));

    return (CKE2Doc*)m_pDocument;

}

#endif //_DEBUG

int sgn(float a)

{

int sg;

if (a>0) sg=1;

if (a==0) sg=0;

if (a<0) sg=-1;

return(sg);

}

#define n 6

void CKE2View::drawPoint(CDC *pDC)

{

 

     double **c,*f1,*f,*x,*w,*e,max,p=2,q=2,xx,yy;

 int i=0,j=0,k,m,a,b,*l,bb=0;

 c=new double*[10000];

 for(i=0;i<10000;i++)

 {

 c[i]=new double[3];

 memset(c[i],0,sizeof(double)*3);

 }

 

 f=new double[10000];

 e=new double[10000];

 w=new double[10000];

 f1=new double[10000];

 x=new double[n];

 l=new int[10000];

 

 for(xx=0.5;xx<1;xx+=0.01)

    for(yy=xx;yy<1);yy+=0.01)

    {

p=1./(1.-xx);

q=1./(1.-yy);

memset(w,0,sizeof(double)*10000);

memset(e,0,sizeof(double)*10000);

memset(f1,0,sizeof(double)*10000);

memset(x,0,sizeof(double)*n);

 x[n-1]=1;

 j=0;

 

 for(i=0;i<10000;i++)

 {j=0;

 f1[i]=1;c[i][0]=0;c[i][1]=1;c[i][2]=0.5;

 while(fabs(f1[i])>0.00000001)

 {

    f1[i]=0;

 for(k=0;k<n;k++)

 { f1[i]+=pow((fabs(x[k]-c[i][2])),(p-1))*sgn(x[k]-c[i][2]);}

 if (f1[i]<-0.00000001)

    {max=c[i][2];c[i][2]=c[i][2]-(fabs(c[i][2]-c[i][1])/2.0);c[i][1]=max;}

 if (f1[i]>0.00000001)

    {max=c[i][2];c[i][2]=c[i][2]+(fabs(c[i][2]-c[i][1])/2.0);c[i][1]=max;}

 if (fabs(f1[i])<=0.00000001)

 {c[i][0]=c[i][2];goto B;}

 

 }

B:

 c[i][0]=c[i][2];

 

 for(a=0;a<n;a++)

 {

 

 for(b=0;b<n;b++)

 w[i]+=pow((fabs(x[a]-x[b])),q);

 e[i]+=pow((fabs(x[a]-c[i][0])),p);

 }

 f[i]=pow((e[i]/n),(1./p))/pow((w[i]/(n*n)),(1./q));

 x[n-2]+=0.1;

 for(k=2;k<n;k++)

 {

 if(x[n-k]>1.04)

 {

 x[n-k-1]+=0.1;

 x[n-k]=x[n-k-1];

 for(m=2;m<k;m++)

 x[n-m]=x[n-k-1];

 }

 if (x[0]!=0) goto A;

 }

 

 }

 

A:

 

max=f[0];k=0;

for(m=0;m<i;m++)

{

if (fabs(max-f[m])<0.001) {k++;l[k]=m;}

if (max<f[m]) {k=0;max=f[m];l[k]=m;}

}

 

for(a=0;a<n-1;a++)

    x[a]=0;

for(a=0;a<l[0];a++)

{

    x[n-2]+=0.1;

 for(k=2;k<n;k++)

 if(x[n-k]>1.04)

 {

 x[n-k-1]+=0.1;

 x[n-k]=x[n-k-1];

 for(m=2;m<k;m++)

 x[n-m]=x[n-k-1];

 }

}

 

b=0;

for(k=0;k<n;k++)

    {

    if((x[k]==0)||(fabs(x[k]-1)<0.04)) b++;

    else

    {

              if(fabs(x[k]-0.5)<0.04) b+=2;

              else b=-n;

    }

    }

b-=n;

if (b<0) b=24;

if (b==0) b=58;

if(b==bb) continue;

bb=b;

c=%f\n",p,q,l[0],l[k],k+1,max,c[l[0]][0]);

COLORREF cr(RGB((b%3)*127,(b%4)*85,(b%5)*63));

CBrush r(cr);

CPen rp(PS_SOLID,0,cr);

pDC->SelectObject(&rp);

pDC->SelectObject(&r);

CPoint r1[3]={CPoint(0,360),CPoint(int(720./p),-int(720./q)+360),CPoint(int(720./p),360)};

pDC->Polygon(r1,3);

 

    }

}

 

Приложение 2. Листинг программы №2.

//Покоординатный спуск

#include<stdAfx.h>

#include<stdio.h>

#include<iostream.h>

#include<conio.h>

#include<math.h>

#define n 4

void main()

{

    //double ff(double v,double vv);

    int sgn(float a);

    double *aa,**x,*f1,f,e,w,w1,e1,q=7,p=7,*r,*z,f2,*r1,max,m=0,c,f20,f21;

    int bb,i,MAX,k,j,a,b,mm,ss;

    x=new double*[n];

    for(i=0;i<n;i++)

              x[i]=new double[2];

    z=new double[3];

    aa=new double[n-2];

    r=new double[n];

    r1=new double[n];

    //f2=new double[n];

    f1=new double[n];

        

    for(i=1;i<n-1;i++) //начальное прибл - все х от 1-го до n-2-го равны

              x[i][0]=0.1; //

    x[n-1][0]=1;x[n-1][1]=1;x[0][1]=0;x[0][0]=0;x[1][0]=0.9;//x[2][0]=0;//x[3][0]=0;x[4][0]=1;//начальное приближение

        

    for(c=0.5;c<0.7;c+=0.1) //Цикл по c        

    {       

              bb=0;

              for(k=1;;k++) //Цикл по k

              {

              //  for(i=0;i<n;i++)

              //            cerr<<"x["<<i<<"]="<<x[i][0]<<"\n";

              //  cerr<<"\n";

                       w=0;e=0;w1=0;e1=0;

                       for(a=0;a<n;a++)

                                 r[a]=0;

                       for(a=0;a<n;a++)

                       {

                                 for(b=0;b<n;b++)

                                 {

                                          w+=pow((fabs(x[a][0]-x[b][0])),q);

                                          r[a]+=pow((fabs(x[a][0]-x[b][0])),q-1)*sgn(x[a][0]-x[b][0]);

                                 }

                                 e+=pow((fabs(x[a][0]-c)),p);

                       }       

                       w=pow(w/(n*n),1./q);//знаменатель исх ф-ции

                       e=pow(e/n,1./p);//числитель исх ф-ции

 f=e/w;

                       //cerr<<"\n"<<f<<"\n";

                       f1[0]=0;f1[n-1]=0;

                       MAX=0;

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

                       {

                                 f1[j]=pow(n,(2./q-1./p))*(pow(e,(1./p-1))*pow(fabs(x[j][0]-c),(p-1))*w

                                 *sgn(x[j][0]-c)-2.*pow(w,(1./q-1))*r[j])/(w*w); //производная исх. функции по x[j][0] в точке с координатами x[i][0] i=0..n-1 на k-ом приближении

              //  cerr<<f1[j]<<"\n";

                                 for(a=0;a<bb;a++)

                                          if(aa[a]==j) break;

                                 if(a!=bb) continue;

                        if (fabs(f1[j])>fabs(f1[MAX])) MAX=j;

                       }

                       // т.к. x[0]=0 и x[n-1]=1 - cosnt

 mm=0;ss=0;

                       for(i=0;;i++)

                       {

                                 if (mm==0) z[0]=100000000./pow(1.2,i);//шаг

                                 if(mm==1) 

                                 {

                                          z[0]=-1000000000./pow(1.2,ss);

                                          ss++;

                                 }

                                 /*if(z[0]<0.000000000000001)

                                 {

                                          z[0]=-0.5/pow(1.5,mm);

                                  mm++;

                                 }*/

                                 for(a=0;a<n;a++)

                                          r1[a]=0;

                                 w1=0;e1=0;

                                 for(a=0;a<n;a++)

                                 {

                                          if(a==MAX)

                                          {

                                                   for(b=0;b<n;b++)

                                                   {

                                                             w1+=pow(fabs(x[a][0]+f1[a]*z[0]-(x[b][0]+f1[b]*z[0])),q);

                                                             r1[a]+=pow( fabs(x[a][0]+f1[a]*z[0] - (x[b][0]+f1[b]*z[0]) ),q-1)*

                                                             sgn( x[a][0]+f1[a]*z[0] - (f1[b]*z[0]+x[b][0]) );

                                                   }

                                                   e1+=pow((fabs(x[a][0]+f1[a]*z[0]-c)),p);

 }                                  

                                          else

                                          {

                                                   for(b=0;b<n;b++)

                                                   {

                                                             w1+=pow((fabs(x[a][0]-x[b][0])),q);

                                                             r1[a]+=pow((fabs(x[a][0]-x[b][0])),q-1)*sgn(x[a][0]-x[b][0]);

                                                   }

                                              e1+=pow((fabs(x[a][0]-c)),p);

                                          }       

                                 }

                                 w1=pow(w1/(n*n),1/q);e1=pow(e1/n,1/p);

                                 //printf("\n f=%f f[a]=%f",e/w,e1/w1);

                                 a=0;

                                 //for(j=1;j<n-1;j++)

                                 //if (((x[j][0]+z[0]*f1[j])<1)&&((x[j][0]+z[0]*f1[j])>0)) a++;

                                 //if(a<n-2) continue;

                                 if    (((x[MAX][0]+z[0]*f1[MAX])<1)&&((x[MAX][0]+z[0]*f1[MAX])>0)) a++;

                                 if(a!=1) continue;

                                 if (e1/w1>e/w) break;

                                 if ((e1/w1-e/w)>-0.00000001)

                                 {

                                          if(z[0]>0)

                                          {

                                                   mm=1;

               continue;

                                          }

                                          for(a=1;a<n-1;a++)

                                          x[a][1]=x[a][0];

                                          goto A;

                                 }

                       }

        

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

                                 if(j==MAX) x[j][1]=x[j][0]+z[0]*f1[j];

                                 else x[j][1]=x[j][0];

 if(fabs(x[MAX][1]-x[MAX][0])<0.00000001)

                       {

                                 aa[bb]=MAX;

                                 bb++;

                       }

                       if(bb==n-2) break;

 

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

                                 x[j][0]=x[j][1];

              } //закрытие цикла по k

 

              A:

              cerr<<"K-vo iteratsiy: "<<k-1<<"\n\n";

              for(i=0;i<n;i++)

                       cerr<<"x["<<i<<"]="<<x[i][0]<<"\n";

                       cerr<<"\nf="<<f<<"\n";

    }// закрытие цикла по с

}

 

int sgn(float a)

{

    int sg;

    if (a>0) sg=1;

    if (a==0) sg=0;

    if (a<0) sg=-1;

    return(sg);

}

Приложение 3. Листинг программы №3

 

#include <stdAfx.h>

#include <stdlib.h>

#include <stdio.h>

#include <conio.h>

#include <math.h>

#include <iostream.h>

#define n 4

int sgn(float);

main()

{

 

 double **c,*f1,*f,*x,*w,*e,max,p=2,q=2,xx,yy;

 int i=0,j=0,k,m,a,b,*l;

 

 c=new double*[10000];

 for(i=0;i<10000;i++)

 c[i]=new double[4];

 

 f=new double[10000];

 e=new double[10000];

 w=new double[10000];

 f1=new double[10000];

 x=new double[n];

 l=new int[10000];

 

 for(xx=0.5;xx<1;xx+=0.01)

    for(yy=xx;yy<1;yy+=0.01)

    {

p=1./(1.-xx);

q=1./(1.-yy);

 for(i=0;i<10000;i++)

 {

     w[i]=0;

     e[i]=0;

     f1[i]=0;

 }

 

 for(i=0;i<10000;i++)

 for(j=0;j<3;j++)

 {c[i][j]=0;

 }

 

 for(i=0;i<n;i++)

 x[i]=0;

 x[n-1]=1;

 j=0;

 

 for(i=0;i<10000;i++)

 {j=0;

 f1[i]=1;c[i][0]=0;c[i][1]=1;c[i][2]=0.5;

 while(fabs(f1[i])>0.000000001)

 {

    f1[i]=0;

 for(k=0;k<n;k++)

 { f1[i]+=pow((fabs(x[k]-c[i][2])),(p-1))*sgn(x[k]-c[i][2]);}

 if (f1[i]<-0.000000001)

    {max=c[i][2];c[i][2]=c[i][2]-(fabs(c[i][2]-c[i][1])/2.0);c[i][1]=max;}

 if (f1[i]>0.000000001)

    {max=c[i][2];c[i][2]=c[i][2]+(fabs(c[i][2]-c[i][1])/2.0);c[i][1]=max;}

 if (fabs(f1[i])<=0.000000001)

 {c[i][0]=c[i][2];goto B;}

 

 }

B:

 c[i][0]=c[i][2];

 

 for(a=0;a<n;a++)

 {

 

 for(b=0;b<n;b++)

 w[i]+=pow((fabs(x[a]-x[b])),q);

 e[i]+=pow((fabs(x[a]-c[i][0])),p);

 }

 f[i]=pow((e[i]/n),(1./p))/pow((w[i]/(n*n)),(1./q));

 x[n-2]+=0.1;

 for(k=2;k<n;k++)

 {

 if(x[n-k]>1.04)

 {

 x[n-k-1]+=0.1;

 x[n-k]=x[n-k-1];

 for(m=2;m<k;m++)

 x[n-m]=x[n-k-1];

 }

 if (x[0]!=0) goto A;

 }

 

 }

A:

max=f[0];k=0;

for(m=0;m<i;m++)

{

if (fabs(max-f[m])<0.0001) {k++;l[k]=m;}

if (max<f[m]) {k=0;max=f[m];l[k]=m;}

}

printf("p=%f q=%f max=%f\n",p,q,max);

x[n-1]=1;

for(a=0;a<=n-2;a++)

    x[a]=0;

for(a=0;a<l[0];a++)

    {

    x[n-2]+=0.1;

 for(k=2;k<n;k++)

 {

 if(x[n-k]>1.04)

 {

 x[n-k-1]+=0.1;

 x[n-k]=x[n-k-1];

 for(m=2;m<k;m++)

 x[n-m]=x[n-k-1];

 }

 }

}

for(a=0;a<n;a++)

printf("Nabor:x[%d]=%f ",a,x[a]);

    }

 getch();

 return 0;

}

int sgn(float a)

{

int sg;

if (a>0) sg=1;

if (a==0) sg=0;

if (a<0) sg=-1;

return(sg);

}

Приложение №4 Результаты работы программы №1:

 

n=6

 

 

 

 

n=5

 

n=4

 

0