Текст подпрограммы и версий ( Фортран )
mnk3r.zip 		Тексты тестовых примеров ( Фортран )
tmnk3r.zip
Текст подпрограммы и версий ( Си )
mnk3r_c.zip 		Тексты тестовых примеров ( Си )
tmnk3r_c.zip
Текст подпрограммы и версий ( Паскаль )
mnk3r_p.zip 		Тексты тестовых примеров ( Паскаль )
tmnk3r_p.zip

Подпрограмма:  MNK3R

Назначение

Решение задачи минимизации диффеpенциpуемой функции многих переменных по заданному направлению (задача одномерной минимизации) методами линейной и квадратичной аппроксимации производной по направлению.

Математическое описание

Для решения задачи 

    min  f(x0 - h p0) ,   0 ≤ h ≤ hL < +∞ ,   x0 ∈ EN ,   p0 ∈ EN ,
      h

x0 и  p0 заданы, используются последовательно методы линейной и квадратичной аппроксимации производной по направлению  ( f ' (x0 - h p0), p0 ) .

Минимум  f (x) по направлению  p0 считается найденным, если выполнено одно из следующих условий:

1.  | hK - hK + 1| < EPS1 , где  hK, hK + 1 - два последовательных шага по направлению  p0,  EPS1 - заданная абсолютная точность одномерной минимизации;
2.  | hK - hK + 1| / hK < EPS2, где EPS2 - заданная относительная точность одномерной минимизации;
3.  | fi' (x0 - hK p0)| < EPSGi для всех  i = 1, 2, ..., N, где EPSG - вектоp покомпонентной точности вычисления градиента в точке  x0 - hK p0;
4.  
| ( fi' (x0 - hK p0), p0) | < max   { max  EPSGi, 10 -1 -EM/2} ,
                                                  i = 1, 2, ..., N
   где  EM - максимальный порядок абсолютной величины
   вещественной константы;
5.  K ≥ Kmax ,  где  Kmax - максимально допустимое количество вычислений градиента.

Алгоритм является модификацией алгоритма одномерной минимизации, изложенного в работе:

B.A.Cкоков, Стандартная программа минимизации диффеpенциpуемой функции многих переменных методом сопряженных градиентов, серия "Стандартные программы решения задач математического программирования", вып.10, Изд - во МГУ, 1968, ротапринт.

Использование 

    SUBROUTINE  MNK3R (N, I0, P, HL, H, X, F, FE, G, GE, FUNC,
                                             GRAD, UP, IP, RM) 
   Параметры
N - размерность пространства переменных (тип: целый);
I0 - целый вектоp длины  N, задающий фиксированные на время минимизации компоненты вектоpа  X: если  I0 (I) = 0, то  I - ая компонента вектоpа  X остается равной своему начальному значению;
P - вещественный вектоp длины  N, задающий направление одномерной минимизации  p0;
HL - правый конец отрезка одномерной минимизации  hL (тип: вещественный);
H - вещественная переменная, содержащая на входе заданный начальный шаг одномерного поиска, а на выходе - расстояние от начальной точки поиска до вычисленной точки минимума;
X - вещественный вектоp длины  N, содержащий на входе начальную точку  x0 поиска минимума  f (x) по направлению  p0, а на выходе - точку минимального вычисленного значения  f (x);
F - вещественная переменная, содержащая на входе значение функции в начальной точке  x0, а на выходе - вычисленное минимальное значение функции вдоль направления  p0;
FE - вещественная переменная, содержащая на входе заданную точность вычисления значений функции, а на выходе - точность вычисления функции в точке минимума;
G - вещественный вектоp длины  N, содержащий на входе градиент функции в начальной точке  x0, а на выходе - градиент функции в вычисленной точке минимума;
GE - вещественный вектоp длины  N, содержащий на входе заданную покомпонентную точность вычисления градиента, а на выходе - точность вычисления градиента в точке минимума;
FUNC - имя подпрограммы вычисления значения функции  f (x) (см. замечания по использованию);
GRAD - имя подпрограммы вычисления градиента функции (см. замечания по использованию);
UP - вещественный вектоp длины 6, задающий упpавляющие параметры алгоритма:
UP(1) - заданный признак; если  UP (1) = 2, то для построения линейной аппроксимации используется алгоритм удвоения пробных шагов; если  UP (1) = 1, то используется постоянный шаг (см. замечания по использованию);
UP(2) - заданная точность одномерной минимизации (относительная при  UP (3) = 1, абсолютная при  UP (3) ≠ 1);
UP(3) - признак:  UP (3) = 1, если задана относительная точность одномерной минимизации;  UP (3) ≠ 1, если задана абсолютная точность;
UP(4) - заданное максимально допустимое количество вычислений градиента;
UP(5) - максимальный порядок абсолютной величины вещественной константы;
UP(6) - максимальное количество десятичных значащих цифр в вещественной константе;
IP - целый вектоp длины 4, содержащий на выходе:
IP(1) - количество вычислений функции при одномерной минимизации;
IP(2) - количество вычислений градиента при одномерной минимизации;
IP(3) - признак: если  IP (3) = 0, то вычисление градиента  G в точке минимума не производилось; если  IP (3) = 1, то в  G содержится вычисленное значение градиента в точке минимума;
IP(4) - признак: если  IP (4) = 0, то одномерный минимум функции найден и лежит внутpи отрезка [0, HL]; если  IP (4) = 1, то минимум достигается в точке  HL; если  IP (4) = 2, то функция  f (x) не убывает по направлению  p0 и алгоритм не может обеспечить сходимость (см. замечания по использованию);
RM - вещественный вектоp длины 2 * N, используемый в подпрограмме как рабочий. 
   Версии:  нет 
   Вызываемые подпрограммы:  нет 
   Замечания по использованию
 

Подпрограммы FUNC и GRAD составляются пользователем.

Первый оператор подпрограммы вычисления функции должен иметь вид: 

           SUBROUTINE  FUNC (X, F, FE)

        Параметры
        X   - вещественный вектор длины  N, задающий точку
                пространства, в которой вычисляется значение функции;
        F   - вещественная переменная, содержащая
                вычисленное значение функции в точке  X;
        FE - вещественная переменная, содержащая на входе
                заданную точность вычисления значения функции
                в точке  X; на выходе - достигнутую точность;

Если значение достигнутой точности вычисления функции не известно, то в теле подпрограммы FUNC параметр FE не должен переопределяться.

Первый оператор подпрограммы вычисления градиента функции  f (x) должен иметь вид: 

           SUBROUTINE  GRAD (X, G, GE, I0)

        Параметры
        X   - вещественный вектор длины  N, задающий точку
                пространства, в которой вычисляется градиент;
        G   - вещественный вектоp длины  N, содержащий
                градиент функции в точке  X;
        GE - вещественный вектоp длины  N, содержащий  на
                входе заданную покомпонентную точность
                вычисления градиента функции; на выходе -
                достигнутую точность вычисления градиента;
        I0  - целый вектоp фиксированных компонент,
                упpавляющий вычислением компонент градиента:
                если  I0(I) = 0, то полагается  G(I) = 0 .

Если значение достигнутой точности GE (I) для некоторого  I не известно, то в теле подпрограммы GRAD параметр GE не должен переопределяться.

Идентификаторы подпрограмм вычисления значения функции  f (x) и ее градиента должны быть определены в вызывающей программе оператором EXTERNAL.

Если решается задача безусловной минимизации функции по направлению  p0, то следует задать  HL = C, где  C - достаточно большое представимое в машине вещественное число.

B общем случае наибольшее влияние на эффективность программы оказывает выбор начального шага одномерной минимизации  H и значения параметра  UP (1).

При одномерной минимизации функции  f (x) на направлении  p0 вначале ищется отрезок для построения линейной аппроксимации производной по направлению с помощью постоянного шага (при  UP (1) = 1) или алгоритма удвоения пробных шагов (при  UP (1) = 2). Затем последовательно осуществляется линейная и квадратичные аппроксимации производной по направлению.

Если функция многоэкстремальна, то при некоторых значениях  H алгоритм может не обеспечить сходимость к ближайшей к  x0 точке локального одномерного минимума. В этом случае, а также в случае неубывания функции по направлению  p0 (IP (4) = 2) следует уменьшить начальный шаг  H и задать  UP (1) = 1 .

Используются служебные подпрограммы: MNKU1, MNKU3, MNKU5, MNKU6, MNKO2, MNKO3.

Пример использования 

    min  f(x) ,  x ∈ E1 ,  x0 = 0.0 ,  p0 = -1.0 ,
    f(x)  =  100 e -x + x .
    Точка минимума  f(x) вдоль  p0 pавна   x*  =  - ln(0.01) ,
    f(x*)  =  1 - ln(0.01)

       EXTERNAL  FUNC, GRAD
       DIMENSION  P(1), X(1), G(1), GE(1), UP(6), RM(2), I0(1), IP(4)
       DATA  N, I0(1), P(1), X(1), H, HL, FE, GE(1) /1, 1, -1., 0., 0.5, 
      *                                             1.E + 10, 1.E-07, 1.E-07/
       DATA  UP /2., 1.E-06, 2., 50., 18., 12./
       CALL  FUNC (X, F, FE)
       CALL  GRAD (X, G, GE, I0)
       CALL  MNK3R (N, I0, P, HL, H, X, F, FE, G, GE, FUNC, GRAD, UP,
      *                          IP, RM)
       PRINT 1, IP(1), IP(2), X(1), F
     1 FORMAT (16H BЫЧИCЛ. ФУHKЦИИ, I2, 3X,
      *        18H BЫЧИCЛ. ГPAДИEHTA, I3, 3X, 3H X = , E15.8,
      *        2X, 3H F = , E15.8)

   Подпрограмма вычисления значений функции  f(x):

       SUBROUTINE  FUNC (X, F, FE)
       DIMENSION  X(1)
       F = 100.*EXP(-X(1)) + X(1)
       RETURN
       END

   Подпрограмма вычисления градиента функции:

       SUBROUTINE  GRAD (X, G, GE, I0)
       DIMENSION  X(1), G(1), GE(1), I0(1)
       G(1) = -100.*EXP(-X(1)) + 1.
       RETURN
       END

Результаты счета:

     BЫЧИCЛЕНИЙ ФУHKЦИИ  1        BЫЧИCЛЕНИЙ ГPAДИEHTA  10

     X  =  4.6051700 + 00              F  =  5.6051700 + 00