Сравнительный анализ методов оптимизации
< ε,
то примем x=0,237037034153931 и y=0,407407409218273
Z(x,y)= -0,848148148148148
Сравнив все методы, мы видим, что для данной функции
лучше подходит метод деформированного симплекса, т.к. он быстрее приводит к оптимальному
решению.
3. Условная оптимизация
Задача условной оптимизации в общем случае записывается в
известном виде:
Такая задача оптимизации, кроме целевой функции, включает
дополнительные условия в виде ограничений и граничных условий.
На рисунке 12 представлена фигура, объем, которой
необходимо максимизировать при заданной площади поверхности
Рисунок 12 – Фигура для максимизации объема при заданной
площади поверхности
Найдем полную площадь поверхности данной фигуры(без
верхней поверхности):
,
найдем объем фигуры:
Эта задача представляет собой пример задачи условной
оптимизации: необходимо найти максимальный объем при заданном значении площади
поверхности.
Эту задачу можно решить двумя методами:
Метод преобразования целевой функции,
метод штрафных функций.
3.1 Метод преобразования целевой функции
Т.к. положено ограничение типа равенства, то из этого
ограничения одну переменную выразим через другую и подставим полученную зависимость
в целевую функцию и получим преобразованную целевую функцию, но без
ограничений.
V = 4/3∙a2∙h2+7/3∙h1∙a2 →
max (1)
S = 6∙a∙h1+4∙h2∙a (2)
Выразим a из (2) и подставим в (1), получим:
V = s2∙(4∙h2+7∙h1)/3∙(6∙h1+4∙h2)2
Теперь, задав начальные условия, значение площади
поверхности, и выбрав нужную точность можно решить задачу любым методом
безусловной оптимизации.
Возьмем, например, метод правильного симплекса, и зададим
начальные условия: а=1м, h1=3м, h2=4м, s=34м. Для метода симплекса выберем
точность ε=0,001.
Т.е максимальный объем V=12,7151461307724, при заданной
площади получается при h1 = 2,946875, и h2 = 3,83229490168751
3.2 Метод штрафных функций
Методы штрафных функций относятся к группе непрямых
методов решения задач нелинейного программирования:
f(x) -> min;
gi(x) 0, i 1, ..., k;
hj(x) 0, j 1, ..., m;
a x b.
Они преобразуют задачу с ограничениями в
последовательность задач безусловной оптимизации некоторых вспомогательных
функций. Последние получаются путем модификации целевой функции с помощью
функций-ограничений таким образом, чтобы ограничения в явном виде в задаче оптимизации
не фигурировали. Это обеспечивает возможность применения методов безусловной
оптимизации. В общем случае вспомогательная функция имеет вид
F(x,a) f(x) +rS(x)
Здесь f(x) - целевая функция задачи оптимизации; S(x) -
специальным образом выбранная функция штрафа,где r— множитель, значения
которого можно изменять в процессе оптимизации.. Точку безусловного минимума
функции F(x, a) будем обозначать через х(а).
Среди методов штрафных функций различают методы
внутренней и внешней точки. Согласно методам внутренней точки (иначе называемым
методами барьерных функций), исходную для поиска точку можно выбирать только
внутри допустимой области, а для методов внешней точки как внутри, так и вне
допустимой области (важно лишь, чтобы в ней функции целевая и ограничений были
бы определены).
3.2 Методы штрафных функций
Эти методы применяются для решения задач нелинейного
программирования с ограничениями-неравенствами.
В рассматриваемых методах функции Ф(x, а) подбирают
такими, чтобы их значения неограниченно возрастали при приближении к границе допустимой
области G (Рисунок 14). Иными словами, приближение к границе “штрафуется”
резким увеличением значения функции F(x, а). На границе G построен “барьер”,
препятствующий нарушению ограничении в процессе безусловной минимизации F(x,
a). Поиск минимума вспомогательной функции F(x, а) необходимо начинать с
внутренней точки области G .
Таким образом, внутренняя штрафная функция Ф(х, а) может
быть определена следующим образом:
Здесь dG -граница области G.
Рисунок 14 - Внутренняя штрафная функция
Методы внешних штрафных функций
Данные методы применяются для решения задачи оптимизации
при наличии как ограничений-неравенств, так и ограничений-равенств.
В рассматриваемых методах функции Ф(х, а) выбирают
такими, что их значения равны нулю внутри и на границе допустимой области G, а
вне ее -положительны и возрастают тем больше, чем сильнее нарушаются ограничения
(Рисунок 15). Таким образом, здесь “штрафуется” удаление от допустимой области
G.
Рисунок – 15 Внешняя штрафная функция
Внешняя штрафная функция Ф(х, а) в общем случае может
быть определена следующим образом:
Для данного курсового проекта штрафная функция для объема
данной фигуры имеет вид:
,
где - параметр штрафа, С – полная
площадь поверхности, заданная изначально, V(a,h1,h2) = 4/3∙a2∙h2+7/3∙h1∙a2,
S(a,h1,h2) = 6∙a∙h1+4∙h2∙a.
Задача была решена методом правильного трехмерного
симплекса.
Мы видим, что при увеличении значения параметра штрафа,
значение функции уменьшается (ухудшается), а при уменьшении – увеличивается
(улучшается).
4. Симплекс таблицы
Для его применения необходимо, чтобы знаки в ограничениях
были вида «меньше либо равно», а компоненты вектора b - положительны.
Алгоритм решения сводится к следующему:
Приведение системы ограничений к каноническому виду путём
введения дополнительных переменных для приведения неравенств к равенствам.
Если в исходной системе ограничений присутствовали знаки
«равно» или «больше либо равно», то в указанные ограничения добавляются искусственные
переменные, которые так же вводятся и в целевую функцию со знаками,
определяемыми типом оптимума.
Формируется симплекс-таблица.
Рассчитываются симплекс-разности.
Принимается решение об окончании либо продолжении счёта.
При необходимости выполняются итерации.
7 На каждой итерации определяется вектор, вводимый в базис,
и вектор, выводимый из базиса. Пересчитывается таблица.
Дана функция вида:
f(x)=4x1+2x2
Подберем k геометрическим способом решения так, чтобы
область допустимых значений была пятиугольником. k=7
Рисунок – 16 Область допустимых значений
Приведем запись задачи линейного программирования к
стандартной форме, введем новых переменных, все ограничения кроме ограничения
на знак представим в виде равенств, тогда эта задача примет вид.
4у1+2у2+0у3 +0у4 +0у5
=(0;0;8;12;7) – начальные
допустимые базисные решения
Имея начальный базис, составляем симплекс таблицу для
нулевой итерации.
Итерация |
Базисная переменная |
Значение |
у1 |
у2 |
у3 |
у4 |
у5 |
0 |
у3 |
8 |
1 |
2 |
1 |
0 |
0 |
У4 |
12 |
4 |
1 |
0 |
1 |
0 |
У5 |
7 |
2 |
1 |
0 |
0 |
1 |
-f |
0 |
4 |
2 |
0 |
0 |
0 |
Вводим в базис у1 , а выводим из базиса у4.
Итерация |
Базисная переменная |
Значение |
у1 |
у2 |
у3 |
у4 |
у5 |
1 |
У3 |
5 |
0 |
1,75 |
1 |
-0,25 |
0 |
У1 |
3 |
1 |
0,25 |
0 |
0,25 |
0 |
У5 |
1 |
0 |
0,5 |
0 |
-0,5 |
1 |
-f |
-12 |
0 |
1 |
0 |
-1 |
0 |
Вводим в базис у2 , а выводим из базиса у5.
Итерация |
Базисная переменная |
Значение |
у1 |
у2 |
у3 |
у4 |
у5 |
2 |
У3 |
1,5 |
0 |
0 |
1 |
-2 |
-3,5 |
У1 |
2,5 |
1 |
0 |
0 |
0,5 |
-0,5 |
У2 |
2 |
0 |
1 |
0 |
-1 |
2 |
-f |
-14 |
0 |
0 |
0 |
0 |
-2 |
Т.к. f<0, то останавливаемся на второй итерации.
Исходя из графика ОДЗ, можно определить, что оптимальным
решением является отрезок прямой , входящий в
ОДЗ, проверим: 2,5*2+2=7.
x1 = 2,5, x2 = 2 f(x)=14.
Заключение
Целью данного курсового проекта было изучение методов
оптимизации функции. Методов одномерной оптимизации: метод дихотомии, золотого
сечения; многомерной безусловной оптимизации: покоординатный циклический спуск,
метод Хука – Дживса, правильный симплекс, деформированный симплекс, а также
методов условной оптимизации Метод преобразования целевой функции, метод
штрафных функций, табличный симплекс – метод.
Список используемой литературы
1.
А.Г.Трифонов. Постановка задачи оптимизации и численные методы ее решения;
2.
Б. Банди. Методы оптимизации. Вводный курс., 1988;
3.
Мендикенов К.К. Лекции
Приложение А
using System;
using System.Collections.Generic;
using System.ComponentModel;
using System.Data;
using System.Drawing;
using System.Text;
using System.Windows.Forms;
namespace lab1
{
public partial class Form1 : Form
{
class global
{
private global() { }
public static double a = 0.64;
public static double b = 1.77;
public static double e = 0.0001;
public static double al = 0.00001;
public static double x = 0;
public static double y = 0;
public static int iter = 0;
}
public Form1()
{ InitializeComponent(); }
private void textBox1_TextChanged(object
sender, EventArgs e)
{global.e =
Convert.ToDouble(textBox1.Text); }
private void textBox2_TextChanged(object
sender, EventArgs e)
{ global.al =
Convert.ToDouble(textBox2.Text); }
public double F(double x)
{ return (Math.Pow((2.5 - x), 2) + 3.1 *
x); }
public double Z(double x, double y)
{ return (2.5 * Math.Pow(x, 2) + 2 * x *
y + 3.1 * Math.Pow(y, 2) -2* x-3*y); }
public double Dixotom()
{
global.iter = 1;
global.a =
Convert.ToDouble(textBox4.Text);
global.b =
Convert.ToDouble(textBox3.Text);
richTextBox1.Text =
richTextBox1.Text+"a="+Convert.ToString(global.a)+"; b="+Convert.ToString(global.b)+(char)13;
while (true)
{
double x1 =
(global.a+global.b)/2-global.al;
double x2 = (global.a + global.b) / 2 +
global.al;
if (F(x1) < F(x2)) global.b = x2;
else global.a = x1;
richTextBox1.Text = richTextBox1.Text +
Convert.ToString(global.iter) + ") x1=" + Convert.ToString(x1) +
"; x2=" + Convert.ToString(x2) + "; f(x1)=" +
Convert.ToString(F(x1)) + "; f(x2)=" + Convert.ToString(F(x2)) +
"; a=" + Convert.ToString(global.a) + "; b=" +
Convert.ToString(global.b) + (char)13;
global.iter++;
if (Math.Abs(global.b - global.a) <
global.e) break;
} return (global.a + global.b) / 2;
}
public double Zolot()
{ global.iter = 1;
global.a =
Convert.ToDouble(textBox4.Text);
global.b =
Convert.ToDouble(textBox3.Text);
richTextBox1.Text = richTextBox1.Text +
"a=" + Convert.ToString(global.a) + "; b=" +
Convert.ToString(global.b) + (char)13;
double x2 = global.a+0.618*(global.b -
global.a) ;
double x1 = global.a + (1-0.618) *
(global.b - global.a);
while (true)
{
if (Math.Abs(global.b - global.a) <
global.e) break;
richTextBox1.Text = richTextBox1.Text +
Convert.ToString(global.iter) + ") a=" + Convert.ToString(global.a) +
"; b=" + Convert.ToString(global.b) + "; x1=" +
Convert.ToString(x1) + "; x2=" + Convert.ToString(x2) + ";
f(x1)=" + Convert.ToString(F(x1)) + "; f(x2)=" +
Convert.ToString(F(x2)) + (char)13;
if (F(x2) > F(x1))
{ global.b = x2; x2 = x1; x1 = global.a
+ 0.372 * (global.b - global.a); }
else { global.a = x1; x1 = x2; x2 =
global.a + 0.618 * (global.b - global.a); }
global.iter++;
}
return (global.a + global.b) / 2;
}
private void button1_Click(object
sender, EventArgs e)
{ richTextBox1.Text = "";
global.al =
Convert.ToDouble(textBox2.Text);
global.e =
Convert.ToDouble(textBox1.Text);
if (radioButton1.Checked) global.x =
Dixotom();
if (radioButton2.Checked) global.x =
Zolot();
label2.Text = "Минимум: x*=" + Convert.ToString(global.x) + "; y(x*)=" +
Convert.ToString(F(global.x)) + ", число итераций:
"+Convert.ToString(global.iter-1);
}
public void Spusk(double x,double y)
{
while (true)//идем вправо
{ x = x + global.al; if (Z(x, y) >
Z(x - global.al, y)) break; global.iter++;
richTextBox2.Text = richTextBox2.Text + Convert.ToString(global.iter)+
") x=" + Convert.ToString(x) + "; y=" + Convert.ToString(y)
+ "; z(x,y)=" + Convert.ToString(Z(x, y)) + "; al=" +
Convert.ToString(global.al) + (char)13;
x = x - global.al;//возвращаемся на неудачный шаг
while (true)//идем влево
{ x = x - global.al; if (Z(x, y) >
Z(x + global.al, y)) break; global.iter++;
richTextBox2.Text = richTextBox2.Text +
Convert.ToString(global.iter) + ") x=" + Convert.ToString(x) +
"; y=" + Convert.ToString(y) + "; z(x,y)=" +
Convert.ToString(Z(x, y)) + "; al=" + Convert.ToString(global.al) +
(char)13; }
x = x + global.al;//возвращаемся на неудачный шаг
global.x=x; global.y=y;
SpuskV(x, y);
}
public void SpuskV(double x, double y)
{
while (true)//идем вверх
{ y = y + global.al; if (Z(x, y) >
Z(x, y - global.al)) break; global.iter++;
richTextBox2.Text = richTextBox2.Text +
Convert.ToString(global.iter) + ") x=" + Convert.ToString(x) +
"; y=" + Convert.ToString(y) + "; z(x,y)=" +
Convert.ToString(Z(x, y)) + "; al=" + Convert.ToString(global.al) +
(char)13; }
y = y - global.al;//возвращаемся на неудачный шаг
while (true)//идем вниз
{ y = y - global.al; if (Z(x, y) >
Z(x, y + global.al)) break; global.iter++;
richTextBox2.Text = richTextBox2.Text +
Convert.ToString(global.iter) + ") x=" + Convert.ToString(x) +
"; y=" + Convert.ToString(y) + "; z(x,y)=" +
Convert.ToString(Z(x, y)) + "; al=" + Convert.ToString(global.al) +
(char)13; }
y = y + global.al;//возвращаемся на неудачный шаг
global.x = x; global.y = y;
if (global.al/2 > global.e) {
global.al = global.al / 2; Spusk(x, y); }
}
public void Hyg(double x, double y)
{ while (true)
{int min=Vibor(x, y);
if (min == 1) { x = x + 2 * global.e; y
= y + 2 * global.e; if (Z(x - 2 * global.e, y - 2 * global.e) < Z(x, y)) break;
}
if (min == 2) { x = x - 2 * global.e; y
= y + 2 * global.e; if (Z(x + 2 * global.e, y - 2 * global.e) < Z(x, y))
break; }
if (min == 3) { x = x - 2 * global.e; y
= y - 2 * global.e; if (Z(x + 2 * global.e, y + 2 * global.e) < Z(x, y))
break; }
if (min == 4) { x = x + 2 * global.e; y
= y - 2 * global.e; if (Z(x - 2 * global.e, y + 2 * global.e) < Z(x, y))
break; }
global.iter++;
richTextBox2.Text = richTextBox2.Text +
Convert.ToString(global.iter) + ") x=" + Convert.ToString(x) +
"; y=" + Convert.ToString(y) + "; z(x,y)=" +
Convert.ToString(Z(x, y)) +(char)13; }
global.x = x; global.y = y;
}
public int Vibor(double x, double y)
{ int min = 0;
if (Z(x + global.e, y + global.e) <
Z(x, y)) min = 1;
if (Z(x + global.e, y - global.e) <
Z(x, y)) min = 2;
if (Z(x - global.e, y - global.e) <
Z(x, y)) min = 3;
if (Z(x - global.e, y + global.e) <
Z(x, y)) min = 4;
return min; }
public void Sym(double x, double y)
{ double x0 = x; double y0 = y; double
x1 = x0 + global.al; double y1 = y0;
double x2 = x0 + (global.al) / 2; double
y2 = y0 + global.al * Math.Sin(60);
richTextBox2.Text = richTextBox2.Text +
Convert.ToString(global.iter) + ") z(x0,y0)=" +
Convert.ToString(Z(x0, y0)) + " z(x1,y1)=" + Convert.ToString(Z(x1,
y1)) + " z(x2,y2)=" + Convert.ToString(Z(x2, y2)) + " al="
+ Convert.ToString(global.al) + (char)13;
while (true)
{ //поиск наименьшего
double mx0 = x0; double my0 = y0; double
mx1 = x1; double my1 = y1; double mx2 = x2; double my2 = y2;
double z1 = Z(mx0, my0); double z2 =
Z(mx1, my1); double z3 = Z(mx2, my2);
if ((z1 < z2) && (z2 < z3)
&& (z3 > z1)) { x0 = mx0; x1 = mx1; x2 = mx2; y0 = my0; y1 = my1; y2
= my2; }
if ((z1 < z2) && (z2 > z3)
&& (z3 > z1)) { x0 = mx0; x1 = mx2; x2 = mx1; y0 = my0; y1 = my2; y2
= my1; }
if ((z1 > z2) && (z2 < z3)
&& (z3 < z1)) { x0 = mx1; x1 = mx2; x2 = mx0; y0 = my1; y1 = my2; y2
= my0; }
if ((z1 > z2) && (z2 < z3)
&& (z3 > z1)) { x0 = mx1; x1 = mx0; x2 = mx2; y0 = my1; y1 = my0; y2
= my2; }
if ((z1 < z2) && (z2 > z3)
&& (z3 < z1)) { x0 = mx2; x1 = mx0; x2 = mx1; y0 = my2; y1 = my0; y2
= my1; }
if ((z1 > z2) && (z2 > z3)
&& (z3 < z1)) { x0 = mx2; x1 = mx1; x2 = mx0; y0 = my2; y1 = my1; y2
= my0; }
//проверка на выход
if (global.al <= global.e) break;
while (true)
{ //отражение относительно 3
double kx= (x0+x1)-x2; double ky = (y0 + y1) - y2;
if (Z(x2, y2) > Z(kx, ky)) { x2 = kx;
y2 = ky; global.iter++; break; }
//отражение относительно 2
kx = (x0 + x2) - x1; ky = (y0 + y2) - y1;
if (Z(x1, y1) > Z(kx, ky)) { x1 = kx;
y1 = ky; global.iter++; break; }
//отражение относительно 1
kx = (x1 + x2) - x0; ky = (y1 + y2) - y0;
if (Z(x0, y0) > Z(kx, ky)) { x0 = kx;
y0 = ky; global.iter++; break; }
//уменьшаем треугольник
global.al = global.al / 2;
x1 = (x0 + x1) / 2; y1 = (y0 + y1) / 2;
x2 = (x0 + x2) / 2; y2 = (y0 + y2) / 2;
}
richTextBox2.Text = richTextBox2.Text +
Convert.ToString(global.iter) + ") x0=" + Convert.ToString(x0) +
" x1=" + Convert.ToString(x1) + " x2=" +
Convert.ToString(x2) + "; y0=" + Convert.ToString(y0) + " y1="
+ Convert.ToString(y1) + " y2=" + Convert.ToString(y2) + "
z(x0,y0)=" + Convert.ToString(Z(x0, y0)) + " z(x1,y1)=" +
Convert.ToString(Z(x1, y1)) + " z(x2,y2)=" + Convert.ToString(Z(x2,
y2)) + " al=" + Convert.ToString(global.al) + (char)13; }
global.x = x0; global.y = y0; }
private void button2_Click(object
sender, EventArgs e)
{ global.iter = 0;
global.al =
Convert.ToDouble(textBox7.Text);
global.e =
Convert.ToDouble(textBox8.Text);
if (radioButton4.Checked) {Spusk(Convert.ToDouble(textBox6.Text),Convert.ToDouble(textBox5.Text));
}
if (radioButton3.Checked)
Hyg(Convert.ToDouble(textBox6.Text),
Convert.ToDouble(textBox5.Text));
if (radioButton5.Checked)
Sym(Convert.ToDouble(textBox6.Text), Convert.ToDouble(textBox5.Text));
label9.Text = "Минимум: (" + Convert.ToString(global.x) + ";" +
Convert.ToString(global.y) + " f(x,y)=" +
Convert.ToString(Z(global.x, global.y)) + "), число итераций: " + Convert.ToString(global.iter); } }}
Приложение Б
procedure TForm1.Button1Click(Sender:
TObject);
begin
a:=false;
l:=0.05;
al:=1; e:=0.01; gm:=2; bt:=0.5;
x:=1.16166; y:= 1.15185;
iter:=0;
xl:= x; yl:= y; xg:= xl + l; yg:= yl;
xh:= xl + l / 2; yh:= yl + l * Sin(60);
Sym(x,y);
ZZ(x,y); //считаем значение функции в найденной точке
end;
procedure TForm1.Sym(x,y:real);
label
ext,B;
begin
//поиск наименьшего
B: mx0:= xl; my0:= yl; mx1:= xg; my1:=
yg; mx2:= xh; my2:= yh;
ZZ(mx0, my0); z1:=z; ZZ(mx1, my1);
z2:=z; ZZ(mx2, my2); z3:=z;
if ((z1 < z2) and (z2 < z3) and
(z3 > z1)) then begin xl:= mx0; xg:= mx1; xh:= mx2; yl:= my0; yg:= my1; yh:=
my2; end;
if ((z1 < z2) and (z2 > z3) and
(z3 > z1)) then begin xl:= mx0; xg:= mx2; xh:= mx1; yl:= my0; yg:= my2; yh:=
my1; end;
if ((z1 > z2) and (z2 < z3) and
(z3 < z1)) then begin xl:= mx1; xg:= mx2; xh:= mx0; yl:= my1; yg:= my2; yh:=
my0; end;
if ((z1 > z2) and (z2 < z3) and
(z3 > z1)) then begin xl:= mx1; xg:= mx0; xh:= mx2; yl:= my1; yg:= my0; yh:=
my2; end;
if ((z1 < z2) and (z2 > z3) and
(z3 < z1)) then begin xl:= mx2; xg:= mx0; xh:= mx1; yl:= my2; yg:= my0; yh:=
my1; end;
if ((z1 > z2) and (z2 > z3) and
(z3 < z1)) then begin xl:= mx2; xg:= mx1; xh:= mx0; yl:= my2; yg:= my1; yh:=
my0; end;
Richedit1.Lines.Add(IntToStr(iter));
Richedit1.Lines.Add(FloatToStr(xl)+' '+FloatToStr(yl)+'
'+FloatToStr(zl));
Richedit1.Lines.Add(FloatToStr(xg)+' '+FloatToStr(yg)+'
'+FloatToStr(zg));
Richedit1.Lines.Add(FloatToStr(xh)+' '+FloatToStr(yh)+'
'+FloatToStr(zh));
Richedit1.Lines.Add('');
x0:=(xl+xg)/2; y0:=(yl+yg)/2;
xr:=(1+al)*x0-al*xh; yr:=(1+al)*y0-al*yh;
ZZ(xl,yl); zl:=z; ZZ(xg,yg); zg:=Z;
ZZ(xh,yh); zh:=z; ZZ(xr,yr); zr:=z;
if zr<zl then
{1} begin
//растяжение
xe:=gm*xr+(1-gm)*x0; ye:=gm*yr+(1-gm)*y0;
ZZ(xe,ye); ze:= z;
if ze<zl then
{2} begin
xh:=xe; yh:=ye; ZZ(xh,yh); zh:=z;
xl:=gm*xl+(1-gm)*x0;
yl:=gm*yl+(1-gm)*y0;ZZ(xl,yl); zl:=z;
xg:=gm*xg+(1-gm)*x0;
yg:=gm*yg+(1-gm)*y0; ZZ(xg,yg); zg:=z;
Shod(xl,xg,xh,yl,yg,yh);
if a=true then goto ext else begin
inc(iter); goto B; end;
{2} end
else
{3} begin
xh:=xr; yh:=yr; ZZ(xh,yh); zh:=z;
Shod(xl,xg,xh,yl,yg,yh);
if a=true then goto ext else begin
inc(iter); goto B; end;
{3} end;
{1} end;
if ((zr>zl)and(zr<=zg)) then
begin
xh:=xr; yh:=yr; ZZ(xh,yh); zh:=z;
Shod(xl,xg,xh,yl,yg,yh);
if a=true then goto ext else begin
inc(iter); goto B; end;
end;
if ((zr>zl)and(zr>zg)) then
{4E} begin
if zr<zh then begin xh:=xr; yh:=yr;
ZZ(xh,yh); zh:=z; end;
//сжатие
xc:=bt*xh+(1-bt)*x0; yc:=bt*yh+(1-bt)*y0;
ZZ(xc,yc); zc:=z;
if zc<xh then
{5Ж} begin
xh:=xc; yh:=yc; ZZ(xh,yh); zh:=z;
xl:=bt*xl+(1-bt)*x0;
yl:=bt*yl+(1-bt)*y0;ZZ(xl,yl); zl:=z;
xg:=bt*xg+(1-bt)*x0;
yg:=bt*yg+(1-bt)*y0; ZZ(xg,yg); zg:=z;
Shod(xl,xg,xh,yl,yg,yh);
if a=true then goto ext else begin
inc(iter); goto B; end;
{5} end
{6З} else begin
//уменьшение
l:=l/2;
xh:=(xl+xh)/2; xh:=(xl+xh)/2; ZZ(xh,yh);
zh:=z;
xg:=(xl+xg)/2; yg:=(yl+yg)/2; ZZ(xg,yg);
zg:=z;
{6} end;
Shod(xl,xg,xh,yl,yg,yh);
if a=true then goto ext else begin
inc(iter); goto B; end;
{4} end;
ext:
Richedit1.Lines.Add(FloatToStr(xl));
Richedit1.Lines.Add(FloatToStr(yl));
Richedit1.Lines.Add(FloatToStr(zl));
Richedit1.Lines.Add(IntToStr(iter));
end;
procedure TForm1.ZZ(x,y:real);
begin
z:= 2.5 * x*x + 2 * x * y + 3.1 * y*y
-2* x-3*y;
end;
procedure
TForm1.Shod(xl,xg,xh,yl,yg,yh:real);
var
sigma:real;
begin
ZZ(xl,yl); zl:=z; ZZ(xg,yg); zg:=Z;
ZZ(xh,yh); zh:=z;
sigma:=sqrt((sqr(zl-((zl+zg+zh)/2+1))+sqr(zl-((zl+zg+zh)/2+1))+sqr(zg-((zl+zg+zh)/2+1))+sqr(zh-((zl+zg+zh)/2+1)))/3);
if sigma<e then a:=true;end;
Приложение В
unit Unit1;
type
procedure Button1Click(Sender: TObject);
procedure ZZ(x,y,z:real);
procedure Sym(x,y,z:real);
end;
var
X, y,z:real;
z_:real; //значение функции
iter:integer; //число итераций
s,gm,x0,x1,x2,x3:real; //координаты точек
y0,y1,y2,y3:real;
z0,z1,z2,z3:real; // треугольника
al, e:real; // длина стороны симпликса(треугольника)
kx,ky,kz, zz1,zz2:real; //координаты точки k
a,b:boolean; //для цикла
implementation
procedure TForm1.Sym(x,y,z:real);
label
l;
var
a:array[1..4,1..4] of real;//z()xyz ;
1234
k,i:integer;
buf:real;
changed:boolean;
{1} begin
x0:= x; y0:= y;z0:=z;
x1:= x0 + al; y1:= y0;z1:=z;
x2:= x0 + al / 2; y2:= y0 + al *
Sin(60);z2:=z;
x3:= x0 + al/ 2; y3:=y0 ;z3:=z0 + al *
Sin(60) ;
while (true) do
//for i:=1 to 100 do
{2} begin
ZZ(x0, y0,z0); a[1,1]:=z_; a[2,1]:=x0;
a[3,1]:=y0; a[4,1]:=z0;
ZZ(x1, y1,z1); a[1,2]:=z_; a[2,2]:=x1;
a[3,2]:=y1; a[4,2]:=z1;
ZZ(x2, y2,z2); a[1,3]:=z_; a[2,3]:=x2;
a[3,3]:=y2; a[4,3]:=z2;
ZZ(x3, y3,z3); a[1,4]:=z_; a[2,4]:=x3;
a[3,4]:=y3; a[4,4]:=z3;
//сортировка
repeat
changed:=FALSE;
for k:=1 to 3 do
if a[1,k] > a[1,k+1] then
begin
buf := a[1,k]; a[1,k] := a[1,k+1]; a[1,k+1]
:= buf;
buf := a[2,k]; a[2,k] := a[2,k+1]; a[2,k+1]
:= buf;
buf := a[3,k]; a[3,k] := a[3,k+1]; a[3,k+1]
:= buf;
buf := a[4,k]; a[4,k] := a[4,k+1]; a[4,k+1]
:= buf;
changed := TRUE;
end;
until not changed;
x0:=a[2,1]; y0:=a[3,1]; z0:=a[4,1];
x1 :=a[2,2]; y1:=a[3,2]; z1:=a[4,2];
x2:=a[2,3]; y2:=a[3,3]; z2:=a[4,3];
x3 :=a[2,4]; y3:=a[3,4]; z3:=a[4,4];
ZZ(x3, y3,z3);
//проверка на выход
if (al <= e)or(z_<0) then break;
while (true) do
{3} begin
//отражение относительно 0
kx:= 2/3*(x2+x1+x3)-x0; ky:= 2/3*(y2+y1+y3)-y0; kz:=
2/3*(z2+z1+z3)-z0;
ZZ(x0, y0,z0); zz1:=z_; ZZ(kx, ky,kz); zz2:=z_;
if (z_<0) then goto l;
if (zz1 <
zz2)and(kx>0)and(ky>0)and(kz>0)and((6*kx*ky+4*kz*kx)<=s) then begin
x0:= kx; y0:= ky; z0:= kz; inc(iter); break; end;
//отражение относительно 1
kx:= 2/3*(x2+x0+x3)-x1; ky:= 2/3*(y2+y0+y3)-y1; kz:=
2/3*(z2+z0+z3)-z1;
ZZ(x1, y1,z1); zz1:=z_; ZZ(kx, ky,kz);
zz2:=z_; if (z_<0) then goto l;
if (zz1 <
zz2)and(kx>0)and(ky>0)and(kz>0)and((6*kx*ky+4*kz*kx)<=s)then begin
x1:= kx; y1:= ky; z1:= kz;inc(iter); break; end;
//отражение относительно 2
kx:= 2/3*(x0+x1+x3)-x2; ky:= 2/3*(y0+y1+y3)-y2; kz:= 2/3*(z0+z1+z3)-z2;
ZZ(x2, y2,z2); zz1:=z_; ZZ(kx, ky,kz); zz2:=z_;if
(z_<0) then goto l;
if (zz1 <
zz2)and(kx>0)and(ky>0)and(kz>0)and((6*kx*ky+4*kz*kx)<=s)then begin
x2:= kx; y2:= ky; z2:= kz;inc(iter); break; end;
//отражение относительно 3
kx:= 2/3*(x2+x1+x0)-x3; ky:= 2/3*(y2+y1+y0)-y3; kz:=
2/3*(z2+z1+z0)-z3;
ZZ(x3, y3,z3); zz1:=z_; ZZ(kx, ky,kz); zz2:=z_;if
(z_<0) then goto l;
if (zz1 <
zz2)and(kx>0)and(ky>0)and(kz>0)and((6*kx*ky+4*kz*kx)<=s) then begin
x3:= kx; y3:= ky; z3:= kz;inc(iter); break; end;
//уменьшаем треугольник
al:= al / 2;
x1:= x0 + al; y1:= y0;z1:=z0;
x2:= x0 + al / 2; y2:= y0 + al *
Sin(60);z2:=z0;
x3:= x0 + al/ 2; y3:=y0 ;z3:=z0 + al *
Sin(60);break; inc(iter);
{3} end;
{2} end;
l: x:= x3; y:= y3;
{1} end;
procedure TForm1.ZZ(x,y,z:real);
begin
z_:=(4/3*x*x*z+7/3*y*x*x)-gm*sqr(s-(6*x*y+4*z*x)); end;
procedure TForm1.Button1Click(Sender:
TObject);
begin
RichEdit1.Text :='';iter:=0;
gm:=StrToFloat(Edit1.Text);s:=StrToFloat(Edit2.Text);al:=0.05;e:=0.01;
x:=StrToFloat(Edit3.Text);y:=StrToFloat(Edit4.Text);z:=StrToFloat(Edit5.Text);
Sym(x,y,z);ZZ(x3,y3,z3); //считаем значение функции в
найденной точке
RichEdit1.Text:= RichEdit1.Text + 'Число
итераций '+ IntToStr(iter) + #13;
RichEdit1.Text:= RichEdit1.Text + 'x3=
'+ FloatToStr(x1) +'; y3= '+ FloatToStr(y1)+'z3= '+ FloatToStr(z1) + #13;
ZZ(x3,y3,z3); RichEdit1.Text:=
RichEdit1.Text +'f(x,y,z)= '+ FloatToStr(z_)
+'('+FloatToStr(6*x3*y3+4*z3*x)+')'+ #13;end;end.
|