книги / Численные методы решения некорректных задач

Р и с. 4.7

в) программой PTICR0 производится минимизация функционала невязки y(z) на отрезке z - ХА, где X G [0, 1]; результатом работы про­ граммы PTICR0 является величина шага X;

г) изменяется приближенное решение в соответствии с найденным шагом X;

д) снова вычисляются значение оператора на новом приближении, градиент в новой точке, значение невязки и квадрат нормы градиента; е) производится проверка условий окончания работы программы: если условие окончания не выполнены, то проверяется, уменьшилось ли значение функционала невязки <p(z), если не уменьшилось, то програм­ ма завершает свою работу с кодом завершения 1; в противном случае запоминается текущее значение невязки и продолжается итерационный

процесс.

3. Описание программы PTIGR2. Программа PTIGR2 служит для вы­ бора оптимальной вершины t0 соответствующего многогранника Л/с,, т.е. вершины, в которой достигается минимальное значение функционала (.g>t ):

где g —заданный фиксированный вектор. Форма обращения к программе:

CALLPTIGR2 (TOP, G, N,Cl, С2,1C)

Здесь:

ТОР —найденная оптимальная вершина (массив длины N) —результат работы программы.

141

G - градиент функционала невязки g - фиксированный задаваемый вектор длины N.

Cl, С2 — константы, значения которых Сь С2 определяют множест­ во Me­

осуществляется перебором. Координаты вершин Г*7* легко выписывают­ ся в явном виде, как это сделано в § 3 гл. III. При реализации перебора нет необходимости каждый раз явно вычислять координаты векторов

воспользуемся априорной информацией о монотонности и выпуклости этого решения. Размерности сеток по обеим переменным выберем равны­ ми 41 (п = m = 41). В качестве правой части уравнения (1) воспользуем­ ся значениями, вычисляемыми в соответствии с (3), таким образом, имеем точность задания правой части уравнения, равную нулю.

Применение для решения этой задачи программы PTIGR позволило найти приближенное решение z(s), численные значения которых даны на рис. 4.7, содержащем также обращение к программе PTIGR. Значение невязки Ф (z) = |\Az - ||2 на найденном приближенном решении соста­ вило 2,08XIО-7, что соответствует относительной погрешности примерно 0.15% (по отношению к максимуму правой части). Всего для нахождения полученного решения потребовалось 400 итераций.

§ 6. Описание программы решения некорректно поставленных задач на специальных множествах. Метод проекций сопряженных градиентов

Метод, описанный в § 2 гл. 3 для решения интегрального уравнения

(1) на множествах монотонных ограниченных функций Zlc> монотонных выпуклых Z4 и выпуклых неотрицательных функций, основанный на применении метода проекции сопряженных градиентов для минимизации функционала невязки Ф(г) = |\Аг - щ ||£ 2 на соответствующем множест­

ве, реализован в виде программ PTILR, PTILR1, PTILR0, PTILR3, PTILR5, PTILR6, PTILR7, PTILRA, PTILRB.

1. Описание программы PTILR. Головная программа - PTILR. Форма обращения к ней:

CALL PTILR (АК, Uф, А, В, С, D, N, М, Z, AN2, DL, ITER, * ШАХ, С2,1C, О, NU, IERR)

142

Z- входной и выходной параметр. При обращении к программе PTILR

вZ должно быть задано начальное приближение (массив из N чисел), начальное приближение должно быть допустимой точкой, т.е. монотонным, выпуклым или монотонным и выпуклым одновременно, неотрицательным,

а также в случае множества Z i c (1C = 1) оно должно быть ограничено сверху константой С2. После окончания работы программы PTILR массив Z будет содержать найденное приближенное решение уравнения (1).

AN2 —выходной параметр —содержит после окончания работы про­ граммы значение функционала невязки Ф(г) = \\Az - и& ||2 на найденном приближенном решении.

функционал Ф (z) до этого значения.

ITER — выходной параметр — число сделанных программой PTILR больших циклов (число изменений активных ограничений).

ШАХ —входной параметр —максимально допустимое число больших циклов.

С2 —входной параметр —константа, значение которой С2, ограничи­

пуклых вверх функций (ZI ) .

U —рабочий массив длиной не менее пт + 5п2 + 5 п + max{/i, т) + 1, если 1C = 1, то длина рабочего массива не менее пт + 5м2 + 8м + + тах{и, т) + 4.

NU - входной параметр —длина рабочего массива U. IERR —выходной параметр —код завершения:

IERR = 0 - найдена точка минимума функционала невязки Ф(г) на соответствующем множестве (выход на шаг 6 в алгоритме § 2 гл. III); IERR = 1 - достигнуто значение функционала невязки меньшее задан­

ного значения погрешности DL;

IERR = 2 —найдена точка минимума функционала невязки Ф (z) (нор­ ма градиента функционала равна машинному нулю);

143

IERR = 3 — сделано ШАХ больших циклов (изменений множества активных ограничений);

IERR = 64 —не хватает длины рабочего массива U, значения всех вы­ ходных параметров не определены;

IERR = 65 —ошибка в задании множества корректности, задано 1C либо меньшее 1, либо большее 3, либо при 1C = 1 задано значение С2 25 0.0. Зна­ чения всех выходных параметров в этом случае не определены.

Вызываемые подпрограммы: PTILR1, PTILR2, PTILR3, PTILR4, PTILR5, PTILR6, PTILR7, PTILRA, PTILRB, PTICR0, PTICR1, PTICR2, PTICI2, PTICR3, PTICR5, PTICR6, PTILR0.

Краткое описание работы программы PTILR. Операторы с 9 по 26 формируют начала рабочих массивов и проверяют достаточность длины рабочего массива U программы для размещения рабочих массивов. Обра­ щением к программе PTICR0 формируется матрица линейного оператора, действующего из R" в Rw, аппроксимирующего интегральный оператор уравнения (1). Затем программа PTILRA преобразует функционал не­ вязки <p(z) = II Az - us II2 к виду <p(z) = (z, Qz) + (<d, z). + e, т.е. фор­ мирует матрицу Q = A*A и вектор d = 2A*us. Программа PTILRB форми­ рует матрицу F и вектор g ограничений Fz < g в соответствии со значением параметра 1C, определяющего множество корректности.

Основа метода реализована в виде программы PTILR1, служащей для минимизации функционала невязки (z, Qz) + (d, z) + e при наличии ограничений вида Fz < g. Невязка AN2, выдаваемая программой PTILR1, вычисляется в пространстве R'” , и поэтому после обращения к PTILR1 она нормируется в соответствии с метрикой L 2 [с, d] .

Вход PTILRE, служащий для продолжения счета по программе, позво­ лит избежать при повторном обращении к алгоритму формирования матриц оператора и ограничений и преобразования функционала невязки к виду (z, Qz) + (d, z) + e. При этом вместо обращения к программе PTILR1 происходит обращение к ее входу PTILR2.

2. Описание программы P11LR1. Поскольку основой реализации ал­ горитма является, как уже отмечалось, программа PTILR1, представляю­ щая самостоятельный интерес, то приведем краткое ее описание. Програм­ ма PTILR1 служит для минимизации квадратичной функции \\Az - и ||^ w = = (z, Qz) + (d, z) + e при наличии ограничений вида Fz <g. Минимизация проводится методом проекции сопряженных градиентов, как это описано в § 2 гл. III.

Форма обращения:

CALL PTILR1 (М, MF, MASK, A, AJ, F, Р, U0, PI, G, GR,

*W,PI, WORK, IMAX, ITER, Q ,X, D, N, NN, AN2, DL,

*ANGRD, IEND)

Здесь:

M —число точек сетки, на которой задается правая часть и (в случае WORK = .FALSE, можно не задавать).

MF - целая рабочая переменная —число активных ограничений.

MASK - рабочий массив длиной NN —число ограничений; содержит признаки активности ограничений.

144

А - матрица М * N оператора А (в случае WORK = .FALSE, можно не задавать).

AJ - рабочий массив длиной не менее NN * N.

F - матрица ограничений размерности NN * N.

Р - рабочий массив длиной не менее NN * N.

Uф - правая часть уравнения (1) (в случае WORK = .FALSE, можно не задавать).

PI - рабочий массив длиной не менее N * N. G - вектор ограничений g длиной NN.

GR —рабочий массив длиной не менее N.

W - рабочий массив длиной не менее max{NN, М, N}. Р1 —рабочий массив длиной не менее N.

WORK - логическая переменная, определяющая режим вычисления невязок:

WORK = .TRUE. - заданы входные параметры М, А, Ц0, DL; на каждом

М, A, U0, DL при обращении можно не задавать).

ШАХ —максимально допустимое число больших циклов. ITER - число сделанных больших циклов.

Q —матрица Qфункционала невязки \р(z) = (z, Qz) + (d, z) + е.

X - начальное приближение; после завершения программы массив X длины N содержит найденное приближенное решение.

D - вектор d функционала \р(z) = (z, Qz) + (d9z) +e длины N. N - размерность пространства решений.

NN - число ограничений.

AN2 - содержит после окончания работы программы значение функ­ ционала невязки <p(z) = \\Az - и ||^w на найденном приближенном реше­ нии в случае задания параметра WORK = .TRUE.

если квадрат нормы градиента становится менее этого значения, то про­ грамма прекращает свою работу.

IEND - код завершения программы:

IEND * 0 —найдена точка минимума <p(z) на рассматриваемом мно­ жестве (выход на шаг 6 алгоритма § 2 гл. III);

IEND = 1 — достигнуто значение функционала ^(z), меньшее DL; IEND = 2 - достигнуто значение нормы градиента, меньшее ANGRD; IEND = 3 —сделано ШАХ больших циклов.

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

Обращаем внимание на то, что параметры DL и AN2, характеризующие невязку, имеют для программ PTILR и PTILR1 различный смысл.

145

146

Рис. 4.8

Использование программы PTILR1 возможно также при решении за­ дачи (1) при наличии информации о положении, например, экстремумов искомой функции. Для этого достаточно в программе PTILR использо­ вать другую программу PTILRB формирования матрицы F и вектора огра­ ничений g. При этом следует следить, чтобы в любой допустимой точке строки матрицы активных ограничений были линейно независимы. В про­ стейшем случае изменения сводятся к перемене знаков у соответствую­ щих строк матрицы ограничений. Приведем для примера программу PTILRB для случая задачи на множестве функций, неубывающих при a <лг0 и невозрастающих при х0 < &. Предполагается, что точка Хо принадлежит сетке, на которой ищется решение и задается в COMMONблоке своим номером J0.

SUBROUTTINE PTILRB (F, G, N, NN, ITASK, С, IERR)

COMMON J0

С* ПАРАМЕТРЫ ITASK, С ФИКТИВНЫЕ DIMENSION F(NN, N), G(NN)

REAL * 8 F, G, C

INTEGER N, NN, ITASK, IERR, L, I, J0

С* ФОРМИРУЕМ NN - ЧИСЛО ОГРАНИЧЕНИЙ

NN = N - 1

С* ЗАНУЛЯЕМ МАТРИЦУ

CALL PTICR2(F, 0. О, N * NN)

С* ФОРМИРУЕМ ВЕКТОР ОГРАНИЧЕНИЙ CALL PTICR2(G, 0.0, NN)

147

12CONTINUE IERR = 0

RETURN END.

Вкачестве тестового расчета предлагается использовать расчет с по­

ми 41. Правая часть выбирается в соответствии с формулой (3). На рис. 4.8 приведена программа-обращение, реализующая решение этой задачи с использованием априорной информации о монотонном убывании искомой функции и ее ограниченности сверху константой С2 = 1. Там же приведены результаты - приближенное решение z (s).

Отметим, что решение задачи (1) с помощью программы PTILR тре­ бует неоднократного обращения матрицы активных ограничений, что требует значительных затрат машинного времени. Кроме того, задание ’’неудачного” начального приближения (на ребре малой размерности) может привести к необходимости накопления активных ограничений, требующего фиктивных шагов, на каждом из которых обращается матри­ ца активных ограничений. Тем самым для задач большой размерности при­ менение программы PTILR может оказаться неэффективным.

§ 7. Описание программы решения некорректно поставленных задач на специальных множествах. Метод сопряженных градиентов

спроецированием на множество векторов

снеотрицательными компонентами

ван в виде программ PTISR, PTISR1, PTISR2, PTISR3, PTISR4.

1. Описание программы PTISR. Головная программа - PTISR. Форма обращения к ней:

CALL PTISR(AK,U0, А, В, С, D, N, М, Z, AN2, ITER, * DL, IMAX, \С, U, NU, IERR)

148

Здесь:

АК(Х, S) —подпрограмма-функция вычисления ядра К(х, s) уравне­

А, В, С, D —входные параметры —соответственно величины d, Ъ, с, d

обязательно допустимой точкой, т.е. монотонным, выпуклым или моно­ тонным и выпуклым одновременно, а также неотрицательным. После окончания работы программы PTISR массив Z будет содержать найден­ ное приближенное решение.

нала равна машинному нулю или найден минимум в текущем подпростран­ стве и никакое ограничение нельзя исключить из множества активных) ;

149

шее применение метода проекции сопряженных градиентов бесполезно:

этом случае значениям выходных параметров доверять нельзя. Вызываемыеподпрограммы: PTISR 1, PTISR2, PTISR3, PTISR4, PTICRO, PTICR0. PT1CR1, PTICR2, PTICI2, PTICR3, PTICR4, PTICR5, PTICR6, PTICR7, PTICR8, PTICR9.

Кратко опишем работу подпрограммы PTISR. Операторы с 9 но 19 формируют начала рабочих массивов и проверяют достаточность длины рабо­ чего массива U программы PTISR. Затем обращением к программе PTISRjft формируется матрица А линейного оператора, аппроксимирующего ин­ тегральный оператор уравнения (1). Вслед за этим программа PTISR2 преобразует матрицу А и начальное приближение Z в соответствии с заме­ ной переменных, переводящей задачу минимизации в множество II* (см. § 3 гл. III). Программа PTISR1 осуществляет минимизацию функционала

невязки <р(х) = \\Ах - и ||„„ на множестве П+ векторов с неотрицатель-

к

ными компонентами методом проекции сопряженных градиентов. Заме­ тим, что в программе PTISR невязка вычисляется в аналоге интегральной метрики Ь2[ с, d ] , а в программе PTISR1 эта величина вычисляется в обыч­ ной евклидовой метрике я-мерного пространства R". После завершения минимизации экстремальная точка х подвергается обратному преобразо­ ванию z = Тх из П + в множество, определяемое значением параметра 1C. Это преобразование осуществляется программой PTISR3.

Вход PTISRE служит для повторного обращения к программе и для продолжения счета. В этом случае нет необходимости заново формиро­ вать матрицу А и осуществлять ее преобразование и преобразование на­ чального приближения. Матрица А в преобразованном виде сохраняется в первых N * М элементах рабочего массива U. Для сохранения текущего приближения х G II+ используется 26-й оператор, пересылающий его в следующие N элементов массива U. Таким образом, для продолжения счета достаточно сохранить и восстановить первые N * (М + 1) элементов рабочего массива U.

2.Описание программы PTISR1. Основой реализации алгоритма явля­

ется программа PTISR 1, служащая в общем случае для минимизации сглаживающего функционала

Ma [z] = \\Az-uW\m + a ( ||z ||* n + p ||z '||^ - i )

150