«Численные методы оптимизации функции нескольких переменных» Вариант №13

Федеральное бюджетное государственное общеобразовательное учреждение
высшего профессионального образования
Ижевский государственный технический университет им. М.Т. Калашникова
Факультет «Информатика и вычислительная техника»
Кафедра «Вычислительная техника»

Курсовая работа
по дисциплине: «Методы оптимизации»
на тему: «Численные методы оптимизации функции
нескольких переменных»
Вариант №13

Выполнил: студент гр. 8-78-2
__________

Проверил: к.т.н., доцент
__________

Ижевск, 2014
Оглавление

Введение 3
Задание 5
Безусловная оптимизация 6
Задачи одномерной безусловной минимизации 6
Задачи многомерной безусловной минимизации 8
Решение задачи аналитическим методом 9
Решение задачи методом Хука-Дживса с использованием метода дихотомии 10
Решение задачи методом наискорейшего спуска с использованием метода квадратичной аппроксимации 13
Заключение 18
Список используемой литературы 19
Приложение 20
Введение
Оптимизация - целенаправленная деятельность, заключающаяся в получении наилучших результатов при соответствующих условиях[4]. Она находит применение в науке, технике и в другой области человеческой деятельности.
Поиски оптимальных решений привели к созданию специальных математических методов и уже в 18 веке были заложены математические основы оптимизации. Однако до второй половины 20 века методы оптимизации во многих областях науки и техники применялись очень редко, поскольку практическое использование математических методов оптимизации требовало огромной вычислительной работы, которую без ЭВМ реализовать было крайне трудно, а в ряде случаев - невозможно. Особенно большие трудности возникали при решении задач оптимизации из-за большого числа параметров и их сложной взаимосвязи между собой.
При всем многообразии содержания конкретных оптимизационных задач они имеют общую форму. Все эти задачи можно классифицировать как задачи минимизации (максимизации) вещественнозначной функции f(х) n-мерного векторного аргумента x=(x1,x2,…,xn), компоненты которого удовлетворяют системе уравнений hk(x)=0, набору неравенств gj(x)≥0 а также ограничены сверху и снизу, т. е. ai ≤xi ≤bi Функцию f(х) принято называть целевой функцией, уравнения hk(x)=0 — ограничениями в виде равенств, а неравенства gj(x)≥0 — ограничениями в виде неравенств. При этом предполагается, что все фигурирующие в задаче функции являются вещественнозначными, а число ограничений конечно.
Итак, в самом общем виде задача оптимизации может быть сформулирована следующим образом:
f(x) min (max) (1)
при ограничениях
hk(x)=0, k=1, ..., K, (2)
gj(x)≥0, j=1, ..., J, (3)
ai ≤xi ≤bi , i=1,…, n (4)
В задаче (1)-(4), - целевая функция, а множество точек удовлетворяющих ограничениям (2)-(4) – это допустимое множество X.
Задача, в которой нет ограничений, т. е. J=K=0 , ai=-∞ и bi=∞, i=1,…, n называется оптимизационной задачей без ограничений или задачей безусловной оптимизации, в противном случае - задачей оптимизации с ограничениями или задачей условной оптимизации.
Задачи оптимизации классифицируются в соответствии с видом функций f, h, g и размерностью вектора X . Задачи без ограничений, в которых X представляет собой одномерный вектор, называются задачами с одной переменной и составляют простейший, но вместе с тем весьма важный подкласс оптимизационных задач. Задачи условной оптимизации, в которых функции h, g являются линейными, носят название задач с линейными ограничениями. В таких задачах целевые функции могут быть либо линейными, либо нелинейными. Задачи, которые содержат только линейные функции вектора непрерывных переменных X, называются задачами линейного программирования; в задачах целочисленного программирования компоненты вектора X должны принимать только целые значения. Задачи с нелинейной целевой функцией и линейными ограничениями иногда называют задачами нелинейного программирования с линейными ограничениями. Оптимизационные задачи такого рода можно классифицировать на основе структурных особенностей нелинейных целевых функций. Если f(x) — квадратичная функция, то мы имеем дело с задачей квадратичного программирования; если f(x) есть отношение линейных функций, то соответствующая задача носит название задачи дробно-линейного программирования. Если функции f(x), gj(x), hj(x) зависят также от случайных параметров w, где w является элементом пространства случайных параметров Ω, тогда имеем задачу стохастического программирования.
Деление оптимизационных задач на эти классы представляет значительный интерес, поскольку специфические особенности тех или иных задач играют важную роль при разработке методов их решения.
Задание
1. Построить график функции f(x,y) = x2 + 3y2 + 4x - 5y.
2. Решить задачу безусловной оптимизации следующими способами:
2.1. Найти стационарные точки аналитическим методом, исследовать функцию в этих точках;
2.2. Найти точки экстремумов функции численным методом нулевого порядка Хука-Дживса с использованием метода дихотомии;
2.3. Найти точки экстремумов функции градиентным методом наискорейшего спуска с использованием метода квадратичной аппроксимации.


Безусловная оптимизация

Задача безусловной оптимизации состоит в нахождении минимума или максимума функции в отсутствие каких-либо ограничений . Несмотря на то, что большинство практических задач оптимизации содержит ограничения, изучение методов безусловной оптимизации важно с нескольких точек зрения. Многие алгоритмы решения задачи с ограничениями предполагают сведение ее к последовательности задач безусловной оптимизации. Другой класс методов основан на поиске подходящего направления и последующей минимизации вдоль этого направления. Обоснование методов безусловной оптимизации может быть естественным образом распространено на обоснование процедур решения задач с ограничениями.
Классификация безусловной оптимизации:
Методы одномерной оптимизации:
-Метод половинного деления
-Метод "Золотого сечения"
-Метод Фибоначчи
-Метод Пауэлла
-Метод секущих
-Метод касательной
Методы многомерной оптимизации:
Методы нулевого порядка
-Метод конфигурации Хука-Дживса
-Метод деформированного многогранника
Методы первого порядка
-Метод градиентного спуска
-Метод наискорейшего спуска
-Метод наискорейшего покоординатного спуска
-Метод сопряженных градиентов
Методы второго порядка
-Метод Ньютона
-Метод Ньютона-Рафсона
-Метод Левенберга-Марквардта



Задачи одномерной безусловной минимизации

Для заданной функции решить задачу одномерного поиска f(x)→min(max), x X (X R) и найти промежуток (X R), на котором функция унимодальная.
Пусть f(x) – действительная функция одной переменной, определенная на множестве . Точка x* X называется точкой локального минимума f(x) на множестве X, если существует такая -окрестность этой точки, что для всех x из этой окрестности, т.е., если |x-x*|< выполняется условие f(x*) < f(x). Если выполняется условие f(x*) < f(x), то x* называется точкой строгого локального минимума. У функции может быть несколько локальных минимумов. Точка x* X называется точкой глобального минимума f(x) на множестве X, если для всех x X выполняется условие f(x*) ≤ f(x). Значение функции f(x*) называется минимальным значением f(x) на множестве X, Для нахождения глобального минимума необходимо найти все локальные минимумы и выбрать наименьшее значение.
В дальнейшем будем рассматривать задачу нахождения локального минимума.
Известно, что для того, чтобы точка x* была точкой локального минимума дифференцируемой функции f(x), необходимо, чтобы выполнялось равенство
f (x*) = 0 (5)
Точка x*, удовлетворяющая равенству (5), называется стационарной точкой. Если функция f(x) дважды дифференцируема, то для того, чтобы стационарная точка x* была точкой строгого локального минимума, достаточно, чтобы выполнялось неравенство
f (x*) > 0 (6)
Если дважды дифференцируемая функция f(x) задана на отрезке [a, b], то можно предложить следующий путь решения задачи нахождения глобального минимума:
1. Найти все стационарные точки на отрезке [a, b] из условия (5), т.е. найти корни уравнения f (x) = 0, принадлежащие отрезку [a, b].
2. Найденные стационарные точки исследовать на выполнение условия (6), т.е. из найденных стационарных точек выделить точки локальных минимумов, для
которых выполняется неравенство f (x) > 0.
3. Сравнить между собой значения f(x) на концах отрезка [a, b] и в точках локальных минимумов. Наименьшему из этих значений соответствует точка глобального минимума f(x) на отрезке [a, b].
Пусть функция f(x) определена на отрезке [a, b]. Функция f(x) называется унимодальной, если на этом отрезке имеется единственная точка x* локального минимума функции f(x), причем функция строго убывает при x ≤ x* и строго возрастает при x ≥ x*. Многие алгоритмы минимизации функции одной переменной построены в предположении, что функция унимодальная на некотором отрезке. Этот отрезок будем называть отрезком локализации точки x*.
Из определения унимодальной функции вытекает следующее важное свойство. Пусть f(x) унимодальная функция на отрезке [a, b] и a ≤ x1 < x2 ≤ b. Тогда
если f(x1) ≤ f(x2), то x* [a, x2];
если f(x1) > f(x2), то x* [x1, b], (7)
где x* – точка минимума f(x) на отрезке [a, b].
Иллюстрация свойства (7) представлена на рисунке 1.

Рисунок 1 – Унимодальность функции

Для эвристического выбора начального интервала неопределенности можно применить алгоритм Свенна:
Шаг 1 Задать произвольно следующие параметры: некоторую точку x0, t>0–величину шага. Положить k=0.
Шаг 2 Вычислить значение функции в трех точках x0-t, x0, x0+t
Шаг 3 Проверить условия окончания:
• если f(x0-t) ≥ f(x0) ≤ f(x0+t), то начальный интервал неопределенности найден:[a0,b0] =[x0-t,x0+t];
• если f(x0-t) ≤ f(x0) ≥ f(x0+t), то функция не является унимодальной, а требуемый интервал неопределенности не может быть найден. Вычисления при этом прекращаются (рекомендуется задать другую начальную точку x0);
• если условие окончания не выполняется, то перейти к шагу 4.
Шаг 4 Определить величину ∆:
• если f(x0-t) ≥ f(x0) ≥ f(x0+t), то ∆=t; a0=x0; x1=x0+t; k=1;
• если f(x0-t) ≤ f(x0) ≤ f(x0+t), то ∆=-t; b0=x0; x1=x0-t; k=1.
Найти следующую точку .
Шаг 5 Проверить условие убывания функции
• если f(xk+1)<f(xk) и ∆=t, то a0=xk;
если f(xk+1)<f(xk) и ∆=-t, то b0=xk.
в обоих случаях положить k=k+1 и перейти к шагу 5;
• Если , процедура завершается. При ∆=t положить b0=xk+1, а при ∆=-t положить a0=xk+1. В результате [a0,b0] - искомый начальный интервал неопределенности.




Задачи многомерной безусловной минимизации
Для заданной функции решить задачу многомерной безусловной минимизации f(x1, x2), x X (X R) и начальной точки x0.
Пусть f(x) = f(x1, x2, … , xn) – действительная функция n переменных, определенная на множестве X Rn. Точка x* X называется точкой локального минимума f(x) на множестве X, если существует такая – окрестность этой точки, что для всех x из этой окрестности, т. е., если ||x-x*||< , выполняется условие f(x*) < f(x). Если выполняется условие f(x*) < f(x), то x* называется точкой строгого локального минимума. У функции может быть несколько локальных минимумов. Точка x* X называется точкой глобального минимума f(x) на множестве X, если для всех x X выполняется условие f(x*) ≤ f(x). Значение функции f(x*) называется минимальным значением f(x) на множестве X, Для нахождения глобального минимума необходимо найти все локальные минимумы и выбрать наименьшее значение. В дальнейшем будем рассматривать задачу нахождения локального минимума.
Теорема 1(необходимое условия экстремума первого порядка)
Пусть x* Rn есть точка локального минимума (максимума) функции f(x) на множестве Rn и она дифференцируема в точке х*. Тогда градиент функции в точке х* равен нулю.
f(x*)=0;
Теорема 2 (достаточное условие)
Пусть функция f(x) в точке дважды дифференцируема, ее градиент равен нулю, а матрица Гессе является положительно определенной (отрицательно определенной), т е.
f(x*)=0; H(x*)>0 (H(x* )<0);
Тогда точка х* есть точка локального минимума (максимума)функции на множестве .
Для определения знака матрицы Гессе используется критерий Сильвестра.

Решение задачи аналитическим методом
Найдем минимум функции классическим методом, используя необходимые и достаточные условия.
Для построения графика используем MathCad (рисунок 2).
f(x,y):=x2+3y2+4x-5y


Рисунок 2 – График заданной функции

Сначала найдем стационарные точки функции внутри рассматриваемой области. Для этого вычислим ее частные производные первого порядка функции.
; .
Вычислим стационарные точки, приравняв к нулю производные:

Найдем производные второго порядка функции f(x,y):

Составим матрицу Гессе:
Классифицируем матрицу Гессе, используя критерий Сильвестра:
;
Миноры положительны. Следовательно, матрица является положительно определенной. Используя критерии проверки достаточных условий экстремума, можно сделать вывод: точка (-2; 5/6) является точкой локального минимума. Значение функции в этой точке: f(-2; 5/6)=--6,083.
Решение задачи методом Хука-Дживса с использованием метода дихотомии
Метод Хука-Дживса

Процедура поиска Хука-Дживса представляет собой комбинацию "исследующего поиска" и "поиска по образцу". [3]
Исследующий поиск ориентирован на выявление локального характера поведения целевой функции и определение направлений вдоль "оврагов". Для проведения исследующего поиска необходимо задать величину шага, которая может быть различной для разных координатных направлений и изменяться в процессе поиска. Поиск начинается в некоторой исходной точке. Если значение целевой функции в пробной точке не превышает значение функции в исходной точке, то шаг поиска рассматривается как успешный. В противном случае необходимо вернуться в предыдущую точку и сделать шаг в противоположное направление с последующей проверкой значения целевой функции. После перебора всех N координат исследующий поиск завершается. Полученную в результате точку называют "базовой точкой".
Поиск по образцу заключается в реализации единственного шага из полученной базовой точки xk вдоль прямой, соединяющей эту точку с предыдущей базовой точкой xk-1. Новая точка образца xk+1 определяется в соответствии с формулой
xk+1 = xk + (xk - xk-1) (8)
Как только движение по образцу не приводит к уменьшению целевой функции, точка xk+1 фиксируется в качестве временной базовой точки и вновь проводится исследующий поиск. Если в результате получается точка с меньшим значением целевой функции, чем в точке xk, то она рассматривается как новая базовая точка xk+1.
С другой стороны, если исследующий поиск неудачен, то необходимо вернуться в точку xk и провести исследующий поиск с целью выявления нового направления минимизации. В конечном счете возникает ситуация, когда такой поиск не приводит к успеху. В этом случае требуется уменьшить величину шага путем введения некоторого множителя и продолжить исследующий поиск. Поиск завершается, когда величина шага становится достаточно малой.

Рисунок 3 - Поиск по методу Хука-Дживса

Алгоритм

Шаг 1. Выбрать начальную базисную точку xk =p1 и шаг длиной h1 для каждой переменной x j , j = 1, 2,…, n. Вычисляется значение функции f (p1) в базисной точке p1 .
Шаг 2. Провести исследующий поиск: каждая переменная по очереди изменяется прибавлением длины шага. Т.е. вычисляется значение функции f (p 1 +h1 e 1), где e 1 – единичный вектор в направлении оси x1 . Если это приводит к уменьшению значения функции, то p1 заменяется на p 1 +h1 e1 . В противном случае вычисляется значение функции f (p 1 -h1 e1), и если ее значение уменьшилось, то p 1 заменяем на p 1 -h1 e1 . Если ни один из проделанных шагов не приводит к уменьшению значения функции, то точка p 1 остается неизменной и рассматриваются изменения в направлении оси х 2 , т. е. находится значение функции f (p 1 +h2 e2 ) и т. д. Когда будут рассмотрены все n переменные, мы будем иметь новую базисную точку p2 .
Шаг 3. Был ли исследующий поиск удачным (p 2 ≠ p 1), то перейти к шагу 5; нет - продолжить.
Шаг 4. Проверка на окончание поиска: если длина шага будет уменьшена до заданного малого значения, то прекратить поиск. Текущая точка аппроксимирует точку оптимума x*. Иначе уменьшить длину шага и перейти к шагу 2.
Шаг 5. Провести поиск по образцу: xpk+1 = xk + λ (xk - xk-1).
Шаг 6. Провести исследующий поиск, используя xpk+1 в качестве базовой точки. Пусть xk+1 - полученная в результате точка.
Шаг 7. Если f (xk+1)< f(xk), то положить xk-1= xk, xk= xk+1 и перейти к шагу 5, иначе перейти к шагу 4.

Метод дихотомии

В основе этого метода лежит свойство унимодальной функции (7), благодаря которому можно сокращать отрезок локализации точки минимума.
Пусть f(x) – непрерывная и унимодальная на отрезке [a, b] функция, принимающая во всех точках этого отрезка конечные значения. Пусть число пробных точек x1, x2, …, xn конечно, и для определения каждой точки xk можно использовать информацию о значениях функции во всех предыдущих точках x1, x2, …, xk – 1 . Положим a0 = a, b0 = b. Середина отрезка [a, b] = [a0, b0] находится в точке . Выберем две симметричные точки
x1 = , x2 = (9)
Величина , удовлетворяющая условию 0 < < b – a , является параметром метода, как правило, – малая величина.
Вычислим значения функции в выбранных точках: f(x1) и f(x2). Определим новый отрезок локализации [a1, b1] следующим образом:
если f(x1) ≤ f(x2), то a1 = a0, b1 = x2;
если f(x1) > f(x2), то a1 = x1, b1 = b0.
Далее процедура деления отрезка [a1, b1] повторяется.
Деление продолжают до тех пор, пока половина длины отрезка [an, bn] не станет меньше заданной точности решения задачи , , т. е. пока не выполнится неравенство < . Тогда за приближение x* принимают середину отрезка [an, bn], т.е. x*≈ .

Алгоритм метода дихотомии

Шаг 1. Ввести исходные данные: a, b, , .
Шаг 2. Определить x1 и x2 по формулам (9).
Шаг 3. Вычислить f(x1) и f(x2).
Шаг 4. Если f(x1) ≤ f(x2), то перейти к новому отрезку [a, b], положив b= x2. Иначе перейти к новому отрезку [a, b], положив a = x1.
Шаг 5. Если < , то требуемая точность достигнута, перейти к шагу 6, иначе – к шагу 2 для продолжения итераций.
Шаг 6. Положить x* ≈ .. Вычислить f *≈ f(x*).

Реализация алгоритма (текст программы представлен в приложении):








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

Пусть дана функция f(х), ограниченная снизу на множестве Rn и имею¬щая непрерывные частные производные во всех его точках.
Требуется найти локальный минимум функции f (х) на множестве допус¬тимых решений X = Rn , т.е. найти такую точку х* Rn, что f(x*)= f(x).
Стратегия поиска состоит в построении последовательности точек {хk}, k=0,1.....,таких, что f(хk+1) < f(хk) , k = 0,1... Точки последовательно¬сти вычисляются по правилу
xk+1=xk-tk f(xk),
где точка х° задается пользователем; величина шага tк определяется для каждого значения k из условия
(10)
Решение задачи (5) может осуществляться с использованием необходи¬мого условия минимума и с последующей проверкой достаточного условия минимума . Такой путь может быть использован либо при достаточно простой минимизируемой функции , либо при предварительной аппроксимации достаточно сложной функции полиномом P(tk) как правило, второй или третьей степени), и тогда условие замещается условием , а условие - условием
Другой путь решения задачи (10) связан с использованием численных ме¬тодов, когда ищется . Границы интервала [a,b] задаются пользователем. При этом степень близости най¬денного значения tk к оптимальному значению tk*, удовлетворяющему условиям , , зависит от задания интервала [a,b] и точности методов одномерной минимизации.
Построение последовательности {хk}, k = 0,1., заканчивается в точке хk, для которой где - заданное число, или, если k≥M, М - предельное число итераций, или при двукратном одновременном выполнении неравенств - малое положительное число. Вопрос о том, может ли точка хk рассматриваться как найденное при¬ближение искомой точки локального минимума х* решается путем дополни¬тельного исследования.
Алгоритм метода
Шаг 1 Задать x0, >0, >0 , предельное число итераций M. Найти градиент функции в произвольной точке.
Шаг 2 Положить k=0.
Шаг 3 Вычислить .
Шаг 4 Проверить выполнение критерия окончания
а) если критерий выполнен, то x*=xk.
б) если критерий не выполнен, то перейти к шагу 5.
Шаг 5 Проверить выполнение неравенства k≥M:
а) если неравенство выполнено, то.
б) если критерий не выполнен, то перейти к шагу 6.
Шаг 6 Вычислить величину шага из условия

Шаг 7 Вычислить . xk+1=xk-tk f(xk),
Шаг 8 Проверить выполнение условий
а) если оба условия выполнены при текущем значении k и k= k -1, то расчет окончен, x*=xk+1 ;
б) если хотя бы одно из условий не выполнено, то k= k +1, перейти к шагу 3.
Геометрическая интерпретация метода для n=2 приведена на рисунке 4.

Рисунок 4 – Геометрическая интерпретация метода наискорейшего спуска
Замечания
1. Утверждение гарантирует сходимость последовательности {хk} к стационарной точке х*, где f(xk)=0. Следовательно, найденная в результате применения метода точка х* нуждается в дополнительном исследовании с целью ее классификации.
2. Метод наискорейшего спуска гарантирует сходимость последовательно¬сти {хk} к точке минимума для сильно выпуклых функций.

Метод квадратичной аппроксимации

Метод квадратичной аппроксимации относится к последовательным стратегиям. Задается начальная точка и спомощью пробного шага находятся три точки так, чтобы они были были как можно ближе к искомой точке минимума. В полученных точках вычисляются значения функции. Затем строится интерполяционный полином второй степени, проходящий через имеющиеся три точки. В качестве приближения точки минимума берется точка минимума полинома. Процесс поиска заканчивается, когда полученная точка отличается от наилучшей из трех опорных точек не более чем на заданную величину.
Алгоритм метода
Шаг 1 Задать начальную точку x1 , величину шага ∆x>0, и - малые положительные числа, характеризующие точность.
Шаг 2 Вычислить x2=x1+∆x.
Шаг 3 Вычислить f(x1)=f1 и f(x2)=f2..
Шаг 4 Сравнить f(x1) с f(x2):
а) если f(x1) > f(x2), положить x3=x1+2∆x.
б) если f(x1) ≤ f(x2), положить x3=x1-∆x.
Шаг 5 Вычислить f(x3)=f3 .
Шаг 6 Найти Fmin=min{f1,f2,f3}, xmin=xi : f(xi)=Fmin .
Шаг 7 Вычислить точку минимума интерполяционного полинома, построенног по трем точкам:
,
и величину функции .
Если знаменатель в формуле для на некоторой итерации обращается в нуль, то результатом является прямая. В этом случае рекомендуется обозначит x1=xmin и перейти к шагу2.
Шаг 8 Проверить выполнение условия окончания:
.
Тогда:
а) если оба условия выполнены, процедура заканчивается и x*≈ ;
б) если хотя бы одно из условий не выполнено и [x1,x3], выбирать наилучшую точку (xmin или ) и две точки по обе стороны от нее. Обозначить Эти точки в естественном порядке и перейти к шагу 6;
в) если хотя бы одно из условий не выполнено и [x1,x3], то положить x1= и перейти к шагу 2.

Реализация алгоритма (текст программы представлен в приложении):

















Заключение
В ходе выполнения курсовой работы на практике были рассмотрены 2 алгоритма безусловной оптимизации: метод Хука-Дживса с использованием метода дихотомии и градиентный метод наискорейшего спуска с использованием метода квадратичной аппроксимации. В результате видим, что градиентный метод показывает лучшее быстродействие, в отличие от метода Хука-Дживса.

Список используемой литературы
1) Лекции по курсу «Методы оптимизации».
2) А.В. Пантелеевб, Т.А. Летова Методы оптимизации в примерах и задачах: Учеб. Пособие – 2-е изд., исправл. – М.: Высш. Шк., 2005. – 544с.: ил.
3) Васильев Ф.П. Методы оптимизации. – М.: Издательство «Факториал Пресс», 2002. – 824 с.
4) Рейзлин В.И. Численные методы оптимизации: учебное пособие / В.И. Рейзлин; Томский политехнический университет. - Томск: Изд-во Томского политехнического университета, 2011 - 105 с.

Приложение
unit MO;

interface

uses
Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
Dialogs, StdCtrls, ExtCtrls, ComCtrls;

type
TPointf = record
x: real;
y: real;
end;
TOtrezok = array[0..1]of TPointf;
TGraf = record
zoom:real;
dx :real;
dy :real;
end;
TPLine = ^TLine;
TLine = record
x :real;
y :real;
IsX: boolean;
next :TPLine;
end;
TPDotLine = ^TDotLine;
TDotLine = record
x :real;
y :real;
xn:real;
yn:real;
next :TPDotLine;
end;
TMas = Array[1..2]of real;
TForm1 = class(TForm)
Edit1: TEdit;
Edit2: TEdit;
Edit3: TEdit;
Edit4: TEdit;
Label1: TLabel;
Label2: TLabel;
Label3: TLabel;
Label4: TLabel;
Label5: TLabel;
Label6: TLabel;
Image1: TImage;
bZoom1: TButton;
bZoom2: TButton;
bLeft: TButton;
bRight: TButton;
bUp: TButton;
bDown: TButton;
Edit5: TEdit;
Edit6: TEdit;
Label7: TLabel;
Label8: TLabel;
Edit7: TEdit;
Edit8: TEdit;
Label9: TLabel;
Label10: TLabel;
Edit10: TEdit;
Label12: TLabel;
Edit11: TEdit;
Label13: TLabel;
Label14: TLabel;
Edit12: TEdit;
Label15: TLabel;
Label16: TLabel;
Edit13: TEdit;
Edit14: TEdit;
Label17: TLabel;
Button1: TButton;
Button2: TButton;
Label18: TLabel;
Edit15: TEdit;
Label11: TLabel;
Edit9: TEdit;
Label19: TLabel;
Edit16: TEdit;
Label20: TLabel;
Label21: TLabel;
procedure FormCreate(Sender: TObject);
procedure bLeftClick(Sender: TObject);
procedure bRightClick(Sender: TObject);
procedure bUpClick(Sender: TObject);
procedure bDownClick(Sender: TObject);
procedure bZoom1Click(Sender: TObject);
procedure bZoom2Click(Sender: TObject);
procedure Button1Click(Sender: TObject);
procedure Button2Click(Sender: TObject);
private
{ Private declarations }
procedure Redraw;
procedure AddToLine(var pLine:TPLine; x:real; y:real);
procedure AddToDotLine(var pDotLine:TPDotLine; x:real; y:real; xn:real; yn:real);
procedure XD; //Ìåòîä Õóêà-Äæèâñà
function f(x:real; y:real):real; //Çíà÷åíèå ôóíêöèè
function Dhtm(p1:TMas; p2:TMas; e:real):real;//Ìåòîä äèõîòîìèè
procedure GM; //ãðàèåíòíûé ìåòîä
function df(x:real; y:real):TMas; //ãðàäèåíò ôóíêöèè â òî÷êå
public
{ Public declarations }
end;

const startzoom = 50;
var
Form1: TForm1;
koef: array[0..3]of integer;
Graf: TGraf;
FirstLine: TPLine;
FirstDotLine: TPDotLine;
//M: integer;


implementation

{$R *.dfm}

procedure TForm1.FormCreate(Sender: TObject);
begin
Graf.zoom := 50;
Graf.dx := 0;
Graf.dy := 0;
XD;
//GM;
Redraw;
end;

procedure TForm1.bLeftClick(Sender: TObject);
begin
Graf.dx := Graf.dx + sqr(startzoom)/Graf.zoom;
Redraw;
end;

procedure TForm1.bRightClick(Sender: TObject);
begin
Graf.dx := Graf.dx - sqr(startzoom)/Graf.zoom;
Redraw;
end;

procedure TForm1.bDownClick(Sender: TObject);
begin
Graf.dy := Graf.dy - sqr(startzoom)/Graf.zoom;
Redraw;
end;

procedure TForm1.bUpClick(Sender: TObject);
begin
Graf.dy := Graf.dy + sqr(startzoom)/Graf.zoom;
Redraw;
end;


procedure TForm1.bZoom1Click(Sender: TObject);
begin
Graf.zoom := Graf.zoom/2;
Redraw;
end;


procedure TForm1.Button1Click(Sender: TObject);
begin
XD;
Redraw;
end;

procedure TForm1.Button2Click(Sender: TObject);
begin
GM;
Redraw;
end;

procedure TForm1.bZoom2Click(Sender: TObject);
begin
Graf.zoom := 2*Graf.zoom;
Redraw;
end;

procedure TForm1.Redraw;
var xOy: TPoint;
buf: TPoint;
dot: TPointf;
C: integer;
str: string;
i, j: integer;
g: array[0..1]of integer;
Line: TPLine;
DotLine: TPDotLine;
begin
xOy.x := round(Image1.Width/2+Graf.zoom*Graf.dx/startzoom);
xOy.y := round(Image1.Height/2+Graf.zoom*Graf.dy/startzoom);
g[0] := -round(xOy.x/startzoom);
g[1] := round(xOy.y/startzoom);
With Image1.Canvas do
begin
Brush.Color := ClWhite; //Î÷èñòêà êàíâû
FillRect(Canvas.ClipRect);
//Ñåòêà ----------------------------------------------------------------
Pen.Color := RGB(230,230,230);
Pen.Style := psDot;
for i:=g[0] to g[0]+12 do begin
MoveTo(xOy.x+startzoom*i, 0);
LineTo(xOy.x+startzoom*i, Image1.Height);
end;
for i:=g[1] downto g[1]-12 do begin
MoveTo(0, xOy.y-startzoom*i);
LineTo(Image1.Width, xOy.y-startzoom*i);
end;
//Ëèíèè óðîâíÿ ----------------------------------------------------------
Pen.Color := RGB(230,230,255);
Font.Color := Pen.Color;
Pen.Style := psSolid;
j:=1;
for C:= -6 to 100 do begin
if (C<5)or((C mod 10)=0) then begin
dot.y :=(5-sqrt(25+12*(4+j*C)))/6;
dot.x :=sqrt(abs(5*dot.Y-3*sqr(dot.y)+j*C+4))-2;
str := f(x,y)=+FloatToStr(j*C);
TextOut(xOy.X-20+round(dot.x*Graf.zoom), xOy.Y-round(dot.y*Graf.zoom), str);
MoveTo(xOy.X+round(dot.x*Graf.zoom), xOy.Y-round(dot.y*Graf.zoom));
dot.Y := dot.Y+5/Graf.zoom;
while(dot.y <=(5+sqrt(25+12*(j*C+4)))/6)do begin
dot.x :=sqrt(5*dot.Y-3*sqr(dot.y)+j*C+4)-2;
LineTo(xOy.X+round(dot.x*Graf.zoom), xOy.Y-round(dot.y*Graf.zoom));
dot.Y := dot.Y+5/Graf.zoom;
end;
dot.y :=(5+sqrt(25+12*(j*C+4)))/6;
dot.x :=sqrt(abs(5*dot.Y-3*sqr(dot.y)+j*C+4))-2;
LineTo(xOy.X+round(dot.x*Graf.zoom), xOy.Y-round(dot.y*Graf.zoom));
dot.Y := dot.Y-5/Graf.zoom;
while(dot.y >=(5-sqrt(25+12*(j*C+4)))/6)do begin
dot.x :=-sqrt(5*dot.Y-3*sqr(dot.y)+j*C+4)-2;
LineTo(xOy.x+round(dot.x*Graf.zoom), xOy.Y-round(dot.y*Graf.zoom));
dot.Y := dot.Y-5/Graf.zoom;
end;
dot.y :=(5-sqrt(25+12*(j*C+4)))/6;
dot.x :=-sqrt(abs(5*dot.Y-3*sqr(dot.y)+j*C+4))-2;
LineTo(xOy.x+round(dot.x*Graf.zoom), xOy.Y-round(dot.y*Graf.zoom));
end;
end;
Font.Color := RGB(0,0,0);
//Îñè ------------------------------------------------------------------
Pen.Color := RGB(0,0,0);
Pen.Style := psSolid;
if (xOy.x)<=50 then begin
MoveTo(50,0);
LineTo(50,Image1.Height);
end
else if (xOy.x)>=Image1.Width-50 then begin
MoveTo(Image1.Width-50,0);
LineTo(Image1.Width-50,Image1.Height);
end
else begin
MoveTo(xOy.x,0);
LineTo(xOy.x,Image1.Height);
end;
if (xOy.y)<=50 then begin
MoveTo(0,50);
LineTo(Image1.Width, 50);
end
else if (xOy.y)>=Image1.Height-50 then begin
MoveTo(0,Image1.Height-50);
LineTo(Image1.Width, Image1.Height-50);
end
else begin
MoveTo(0,xOy.y);
LineTo(Image1.Width, xOy.y);
end;
//Ïîäïèñè îñåé ---------------------------------------------------------
buf := xOy;
if (xOy.y)<=50 then
xOy.y := 50
else if (xOy.y)>=Image1.Height-50 then
xOy.y := Image1.Height-50;
for i:=g[0] to g[0]+12 do begin
MoveTo(xOy.x+startzoom*i, xOy.y-2);
LineTo(xOy.x+startzoom*i, xOy.y+2);
TextOut(xOy.x+startzoom*i-4, xOy.y+2, FloatToStr(i*startzoom/Graf.zoom));
end;
xOy.y := buf.y;
if (xOy.x)<=50 then
xOy.x := 50
else if (xOy.x)>=Image1.Width-50 then
xOy.x := Image1.Width-50;
for i:=g[1] downto g[1]-12 do
if i <>0 then begin
MoveTo(xOy.x-2, xOy.y-startzoom*i);
LineTo(xOy.x+2, xOy.y-startzoom*i);
TextOut(xOy.x-15, xOy.y-startzoom*i-4, FloatToStr(i*startzoom/Graf.zoom));
end;
xOy.x := buf.x;
//Ïîäðîáíûé ãðàôèê ýòàïîâ ìèíèìèçàöèè-----------------------------------------
DotLine:=FirstDotLine^.next;
Brush.Color := RGB(255,255,255);
Pen.Color := RGB(255,0,0);
Pen.Style := psDot;
if (DotLine<>nil)then
while(DotLine^.next<>nil) do begin
MoveTo(xOy.X+round(DotLine^.xn*Graf.zoom), xOy.y-round(DotLine^.yn*Graf.zoom));
LineTo(xOy.X+round(DotLine^.x*Graf.zoom), xOy.Y-round(DotLine^.y*Graf.zoom));
DotLine:= DotLine.next;
end;
Pen.Color:= RGB(0,255,0);
Brush.Color:= RGB(0,255,0);
Line:= FirstLine;
MoveTo(xOy.X+round(Line^.x*Graf.zoom), xOy.Y-round(Line^.y*Graf.zoom));
Ellipse(xOy.X+round(Line^.x*Graf.zoom)-3, xOy.Y-round(Line^.y*Graf.zoom)-3,
xOy.X+round(Line^.x*Graf.zoom)+3, xOy.Y-round(Line^.y*Graf.zoom)+3);
TextOut(xOy.x+round(Line^.x*Graf.zoom)+4, xOy.y-round(Line^.y*Graf.zoom)-4, x+FloatToStr(0));
MoveTo(xOy.X+round(Line^.x*Graf.zoom), xOy.Y-round(Line^.y*Graf.zoom));
Pen.Style := psSolid;
j:=0;
while(Line.next<>nil) do begin
Line:= Line.next;
LineTo(xOy.X+round(Line^.x*Graf.zoom), xOy.Y-round(Line^.y*Graf.zoom));
Ellipse(xOy.X+round(Line^.x*Graf.zoom)-3, xOy.Y-round(Line^.y*Graf.zoom)-3,
xOy.X+round(Line^.x*Graf.zoom)+3, xOy.Y-round(Line^.y*Graf.zoom)+3);
if(Line^.IsX) then begin
inc(j);
TextOut(xOy.x+round(Line^.x*Graf.zoom)+4, xOy.y-round(Line^.y*Graf.zoom)-4, x+FloatToStr(j));
MoveTo(xOy.X+round(Line^.x*Graf.zoom), xOy.Y-round(Line^.y*Graf.zoom));
end;
end;
end;
end;

procedure TForm1.AddToLine(var pLine:TPLine; x:real; y:real);
var nLine: TPLine;
begin
new(nLine);
pLine^.next:= nLine;
nLine^.x:= x;
nLine^.y:= y;
nLine^.IsX:=False;
nLine^.next:= nil;
pLine:=nLine;
end;

procedure TForm1.AddToDotLine(var pDotLine:TPDotLine; x:real; y:real; xn:real; yn:real);
var nDotLine: TPDotLine;
begin
new(nDotLine);
pDotLine^.next:= nDotLine;
nDotLine^.x:= x;
nDotLine^.y:= y;
nDotLine^.xn:= xn;
nDotLine^.yn:= yn;
nDotLine^.next:= nil;
pDotLine:=nDotLine;
end;

procedure TForm1.XD; // Ìåòîä Õóêà-Äæèâñà
var dot: array[0..1]of TMas;
p: array[1..3]of TMas;
d: array[1..2]of TMas;
step: TMas;
lmd: real;
alf: real;
e: real;
pLine: TPLine;
pDotLine: TPDotLine;
i: integer;
bool: boolean;
begin
koef[0]:= StrToInt(Edit1.Text);
koef[1]:= StrToInt(Edit2.Text);
koef[2]:= StrToInt(Edit3.Text);
koef[3]:= StrToInt(Edit4.Text);
d[1,1]:=1; d[1,2]:=0;
d[2,1]:=0; d[2,2]:=1;
dot[0][1]:=StrToFloat(Edit5.Text);
dot[0][2]:=StrToFloat(Edit6.Text);
new(FirstLine);
pLine:=FirstLine;
pLine^.x:= dot[0][1];
pLine^.y:= dot[0][2];
pLine^.IsX:= True;
pLine^.next:=nil;
new(FirstDotLine);
pDotLine:=FirstDotLine;
pDotLine^.next:=nil;
step[1]:=StrToFloat(Edit7.Text);
step[2]:=StrToFloat(Edit8.Text);
alf:=StrToFloat(Edit10.Text);
e:=StrToFloat(Edit11.Text);
p[1]:=dot[0];
while(True) do begin
bool:=True;
for i:=1 to 2 do begin
if (f(p[i,1]+step[1]*d[i,1],p[i,2]+step[2]*d[i,2]) < f(p[i,1],p[i,2])) then
begin
p[i+1,1]:= p[i,1]+step[1]*d[i,1];
p[i+1,2]:= p[i,2]+step[2]*d[i,2];
bool:= False;
AddToLine(pLine, p[i+1,1], p[i+1,2]);
if i=2 then pLine^.IsX:= true;
end
else begin
AddToDotLine(pDotLine, p[i,1]+step[1]*d[i,1], p[i,2]+step[2]*d[i,2], p[i,1], p[i,2]);
if(f(p[i,1]-step[1]*d[i,1],p[i,2]-step[2]*d[i,2])<f(p[i,1],p[i,2]))then
begin
p[i+1,1]:= p[i,1]-step[1]*d[i,1];
p[i+1,2]:= p[i,2]-step[2]*d[i,2];
bool:=False;
AddToLine(pLine, p[i+1,1], p[i+1,2]);
if i=2 then pLine^.IsX:= true;
end
else begin
p[i+1,1]:= p[i,1];
p[i+1,2]:= p[i,2];
AddToDotLine(pDotLine, p[i,1]-step[1]*d[i,1], p[i,2]-step[2]*d[i,2], p[i,1], p[i,2]);
end;
end;
end;
if(f(p[3,1],p[3,2])<f(dot[0,1],dot[0,2])) then begin
//ïîèñê
dot[1]:=p[3];
lmd:=dhtm(dot[0],dot[1],e);
p[1,1]:=dot[1,1]+lmd*(dot[1,1]-dot[0,1]);
p[1,2]:=dot[1,2]+lmd*(dot[1,2]-dot[0,2]);
dot[0]:=dot[1];
AddToLine(pLine, p[1,1],p[1,2]);
// pLine^.IsX:= bool;
end
else begin
if ((step[1]<=e)and(step[2]<=e)) then begin
Edit13.Text := FloatToStr(dot[0,1]);
Edit14.Text := FloatToStr(dot[0,2]);
Edit15.Text := FloatToStr(f(dot[0,1],dot[0,2]));
AddToLine(pLine, dot[0,1],dot[0,2]);
pLine.IsX:=True;
break;
end
else begin
for i:=1 to 2 do
if(step[i]>e) then step[i]:=step[i]/alf;
p[1]:=dot[0];
end;
end;
end;
end;

function TForm1.f(x:real; y:real):real;
begin
Result:=koef[0]*sqr(x)+koef[1]*sqr(y)+koef[2]*x+koef[3]*y;
end;

function TForm1.Dhtm(p1:TMas; p2:TMas; e:real):real;
var dlt: real;
xc:real;
ff:TMas;
bool: boolean;
t:real;
tin: TMas;
begin
dlt:=e/5;
//ìåòîä Ñâåíà-----------------------
t:=1;
tin[1]:=0;
tin[2]:=t;
bool:= True;
while(bool) do begin
ff[1]:= f(p2[1]+tin[1]*(p2[1]-p1[1]),p2[2]+tin[1]*(p2[2]-p1[2]));
ff[2]:= f(p2[1]+tin[2]*(p2[1]-p1[1]),p2[2]+tin[2]*(p2[2]-p1[2]));
if(ff[1]>ff[2]) then begin
tin[1]:=tin[2];
t:=t*2;
tin[2]:=tin[2]+t;
end
else bool:=False;
end;
//----------------------------------
while(abs(tin[2]-tin[1])>=e) do begin
xc:= (tin[1]+tin[2])/2;
ff[1]:= f(p2[1]+(xc-dlt)*(p2[1]-p1[1]),p2[2]+(xc-dlt)*(p2[2]-p1[2]));
ff[2]:= f(p2[1]+(xc+dlt)*(p2[1]-p1[1]),p2[2]+(xc+dlt)*(p2[2]-p1[2]));
if(ff[1]<ff[2])then
tin[2]:=xc+dlt
else if(ff[1]>ff[2]) then
tin[1]:=xc-dlt
else begin
tin[2]:=xc+dlt;
tin[1]:=xc-dlt;
end;
end;
Result:= (tin[1]+tin[2])/2;
end;

procedure TForm1.GM;
var e1:real;
e2:real;
M:real;
dot: array[0..1]of TMas;
x: array[1..3]of real;
dx: real;
xmin: real;
xsr: real;
x0: TMas;
ff: TMas;
fi: array[1..3]of real;
fimin: real;
fisr: real;
t: real;
pLine: TPLine;
pDotLine: TPDotLine;
k: integer;
bool, bool2:boolean;
i, j: integer;
begin
koef[0]:= StrToInt(Edit1.Text);
koef[1]:= StrToInt(Edit2.Text);
koef[2]:= StrToInt(Edit3.Text);
koef[3]:= StrToInt(Edit4.Text);
e1:=StrToFloat(Edit9.Text);
e2:=StrToFloat(Edit16.Text);
M:=StrToFloat(Edit12.Text);
new(FirstLine);
pLine:=FirstLine;
dot[0][1]:=StrToFloat(Edit5.Text);
dot[0][2]:=StrToFloat(Edit6.Text);
pLine^.x:= dot[0][1];
pLine^.y:= dot[0][2];
pLine^.IsX:= True;
pLine^.next:=nil;
new(FirstDotLine);
pDotLine:=FirstDotLine;
pDotLine^.next:=nil;
bool:=True;
k:=0;
j:=0;
while(bool) do begin
ff:=df(dot[0,1],dot[0,2]);
if(sqrt(sqr(ff[1])+sqr(ff[2]))<e1)or(k>=M) then begin
x0:=dot[0];
bool:=False;
end
else begin
//ìåòîä êâàäðàòè÷íîé àïïðîêñèìàöèè--------------------------
dx:=5;
x[1]:=0;
bool2:= True;
while(bool2) do begin
x[2]:=x[1]+dx;
fi[1]:= f(dot[0,1]-x[1]*ff[1],dot[0,2]-x[1]*ff[2]);
fi[2]:= f(dot[0,1]-x[2]*ff[1],dot[0,2]-x[2]*ff[2]);
if(fi[1]>fi[2]) then x[3]:=x[1]+2*dx
else x[3]:=x[1]-dx;
fi[3]:= f(dot[0,1]-x[3]*ff[1],dot[0,2]-x[3]*ff[2]);
repeat
fimin:=fi[1]; xmin:=x[1];
for i:=2 to 3 do
if(fimin>fi[i]) then begin
fimin:= fi[i]; xmin:=x[i];
end;
if(((x[2]-x[3])*fi[1]+(x[3]-x[1])*fi[2]+(x[1]-x[2])*fi[3])=0) then
x[1]:=xmin
else begin
xsr:=0.5*((sqr(x[2])-sqr(x[3]))*fi[1]+(sqr(x[3])-sqr(x[1]))*fi[2]
+(sqr(x[1])-sqr(x[2]))*fi[3])/
((x[2]-x[3])*fi[1]+(x[3]-x[1])*fi[2]+(x[1]-x[2])*fi[3]);
fisr:= f(dot[0,1]-xsr*ff[1],dot[0,2]-xsr*ff[2]);
end;
if(abs((fimin-fisr)/fisr)<e2)and(abs((xmin-xsr)/xsr)<e2) then begin
t:=xsr;
bool2:=False;
break;
end
else if(xsr>x[1])and(xsr<x[3]) then begin
if(fisr<fimin) then begin
x[2]:=xsr;
if(xsr>xmin) then x[1]:=xmin
else x[3]:=xmin;
end
else begin
x[2]:=xmin;
if(xsr>xmin) then x[3]:=xsr
else x[1]:=xsr;
end;
end
else begin
x[1]:=xsr;
break;
end;
until(xsr>x[1])and(xsr<x[3]);
end;
//----------------------------------------------------------
// t:=0.1;
dot[1,1]:=dot[0,1]-t*ff[1];
dot[1,2]:=dot[0,2]-t*ff[2];
if(sqrt(sqr(dot[1,1]-dot[0,1])+sqr(dot[1,2]-dot[0,2]))<e2)
and(abs(f(dot[1,1],dot[1,2])-f(dot[0,1],dot[0,2]))<e2) then
inc(j);
if(j>=2)then begin
x0:=dot[0];
bool:=False;
end
else begin
dot[0]:=dot[1];
AddToLine(pLine, dot[0,1],dot[0,2]);
pLine.IsX:=True;
inc(k)
end;
end;
end;
AddToLine(pLine, x0[1],x0[2]);
pLine.IsX:=True;
Edit13.Text := FloatToStr(x0[1]);
Edit14.Text := FloatToStr(x0[2]);
Edit15.Text := FloatToStr(f(x0[1],x0[2]));
end;

function TForm1.df(x:real; y:real):TMas;
begin
Result[1]:=2*koef[0]*x+koef[2];
Result[2]:=2*koef[1]*y+koef[3];
end;

end.