Метод парабол. Поиск экстремума. Pascal

[Dilettante]

Приветствую вас! Стоит задача реализовать метод парабол по поиску экструмума на Паскале. Вот что у меня получилось

Используется синтаксис Pascal

 

program pr7;
uses
    crt;
var
y:array[1..3] of real;       {значения функций в аппроксимирующих точках}
eps1: real;                      {точность}
a, b, c, x, fx, t: real;
fmin: real;                       {минимум функции}
Xmin, Xmiddle: real;        {Значение аргумента при минимуме и среднее значение аргумента}
N:integer;


function f(x:real):real;
begin
f:=exp((1/3)*ln(x-2*x*x+x*x*x));
end;

function midX(var a:real; var c:real; var b:real):real; {поиск среднего занчеия аргумента}
var p1, p2:real;
begin
p1:=(b-c)*(b-c)*(f(a)-f(c))-(c-a)*(c-a)*(f(b)-f(c));
p2:=(b-c)*(f(a)-f(c))+(c-a)*(f(b)-f(c));
midX:=c+0.5*(p1/p2)
end;

BEGIN
//ClrScr;
{поиск экстремума}
{шаг 1}
N:=1;
Writeln('Input accuracy eps1');
Readln(eps1);
Writeln('Input limits a, c and b (a < c < b)');
Readln(a, c, b);
y[1]:=f(a);
y[2]:=f(c);
y[3]:=f(b);
N:=2;
While N <> 4 do
begin
{шаг 2}
while N=2 do
      begin
      t:=midX(a, c, b);      
      if t <> c then x:=t;
      if t = c then x:=(a+c)/2;
      fx:=f(x);
      N:=3;
end;

{шаг 3}
while N=3 do
      begin
      writeln('x ',x);
      if x<c then
         begin
         if fx<f(c) then begin b:=c; c:=x; y[3]:=y[2]; y[2]:=fx; end;
         if fx>f(c) then begin a:=x; y[1]:=fx; end
         else begin a:=x; b:=c; c:=(x+c)/2; y[1]:=fx; y[3]:=y[2]; y[2]:=f(c); end;
      end;
      if x>c then
         begin
         if fx<f(c) then begin a:=c; c:=x; y[1]:=y[2]; y[2]:=fx; end;
         if fx>f(c) then begin b:=x; y[3]:=fx; end
         else begin a:=c; b:=x; c:=(x+c)/2; y[1]:=y[2]; y[3]:=fx; y[2]:=f(c); end;
      end;
      N:=4;
end;

{шаг 4}
while N=4 do
      begin
      if abs(b-a)<eps1 then begin Xmin:=x; fmin:=fx; end
      else N:=2;
end;

end;

Writeln('Min: (',Xmin:0:16,' , ',fmin:0:16,')');
END.
 

Собственно, не работает, зацикливается. Подскажите, вообще хоть что-то тут правильно (по алгоритму поиска экстремума)? И подскажите, кто знает, что я делаю неправильно?

[venco]

А с чего бы ей не зацикливаться. В конце цикла (while N<>4) значение N может быть только 2, т.е. из этого цикла программа никогда не выйдет. Да и из последнего цикла (while N=4) может не выйти, если условие if выполнится.
Кто вас научил такому стилю программирования - необоснованно использовать конечный автомат?
Перепишите без этих скрывающих идею алгоритма циклов while.

[Dilettante]

Я заменил циклы while на if (внешний while на repeat). Но всё же происходит зацикливание, т.е. никогда не выполняется последний if, как я понял. Тогда получается ошибка в сравнении с точностью eps, никогда b-a не стаовяться меньше точности? Просто если взять точность 10, то работает) А снормальными значениями выдаёт зацикливание. Не знаете, как этого избежать?

Используется синтаксис Pascal

 
program pr7;
uses
    crt;
var
y:array[1..3] of real;       
eps1: real;             
a, b, c, x, fx, t: real;
dx: real;                     
fmin: real;                  
Xmin, Xmiddle: real;          
N: integer;
stop: boolean;

function f(x:real):real;
begin
f:=exp((1/3)*ln(x-2*x*x+x*x*x));
end;

function midX(var a:real; var c:real; var b:real):real;
var p1, p2:real;
begin
p1:=(b-c)*(b-c)*(f(a)-f(c))-(c-a)*(c-a)*(f(b)-f(c));
p2:=(b-c)*(f(a)-f(c))+(c-a)*(f(b)-f(c));
midX:=c+0.5*(p1/p2)
end;

BEGIN
//ClrScr;
N:=1;
stop:=false;
Writeln('Input accuracy eps1');
Readln(eps1);
Writeln('Input limits a, c and b (a < c < b)');
Readln(a, c, b);
y[1]:=f(a);
y[2]:=f(c);
y[3]:=f(b);
N:=2;

Repeat
if N=2 then
      begin
      t:=midX(a, c, b);
      //Writeln('t ', t);
      if t <> c then x:=t;
      if t = c then x:=(a+c)/2;
      fx:=f(x);
      N:=3;
end;

if N=3 then
      begin      
      if x<c then
         begin
         if fx<f(c) then begin b:=c; c:=x; y[3]:=y[2]; y[2]:=fx; end;
         if fx>f(c) then begin a:=x; y[1]:=fx; end
         else begin a:=x; b:=c; c:=(x+c)/2; y[1]:=fx; y[3]:=y[2]; y[2]:=f(c); end;
      end;
      if x>c then
         begin
         if fx<f(c) then begin a:=c; c:=x; y[1]:=y[2]; y[2]:=fx; end;
         if fx>f(c) then begin b:=x; y[3]:=fx; end
         else begin a:=c; b:=x; c:=(x+c)/2; y[1]:=y[2]; y[3]:=fx; y[2]:=f(c); end;
      end;
      N:=4;
end;

if N=4 then
      begin
      if abs(b-a)<eps1 then begin Xmin:=x; fmin:=fx; stop:=true; end
      else N:=2;
end;
Until stop=true;

Writeln('Min: (',Xmin:0:16,' , ',fmin:0:16,')');
END.
 

[venco]

Ещё бы N и все if, где она встречается, убрать. Эта переменная совсем не нужна.
А по алгоритму, похоже, в коде, который обновляет a b c не хватает пары else. Если первое условие (fx<f(c)) выполнено, два других варианта выполнять не надо.

[Dilettante]

Исправил дабы вообще не путаться в операторах вот так

Используется синтаксис Pascal

 
program pr7;
uses
    crt;
var
y:array[1..3] of real;        {çíà÷åíèÿ ôóíêöèé â òî÷êàõ x1, x2, x3}
eps1: real;                   {òî÷íîñòè}
a, b, c, x, fx, t: real;
fmin: real;                   {ìèíèìóì ôóíêöèè}
Xmin, Xmiddle: real;          {Çíà÷åíèå àðãóìåíòà ïðè ìèíèìóìå è ñðåäíåå çíà÷åíèå àðãóìåíòà}
stop: boolean;


function f(x:real):real;
begin
f:=exp((1/3)*ln(x-2*x*x+x*x*x));
end;

function midX(var a:real; var c:real; var b:real):real; 
var p1, p2:real;
begin
p1:=(b-c)*(b-c)*(f(a)-f(c))-(c-a)*(c-a)*(f(b)-f(c));
p2:=(b-c)*(f(a)-f(c))+(c-a)*(f(b)-f(c));
midX:=c+0.5*(p1/p2)
end;

BEGIN
//ClrScr;

stop:=false;
Writeln('Input accuracy eps1');
Readln(eps1);
Writeln('Input limits a, c and b (a < c < b)');
Readln(a, c, b);
y[1]:=f(a);
y[2]:=f(c);
y[3]:=f(b);

Repeat
t:=midX(a, c, b);
if t <> c then x:=t;
if t = c then x:=(a+c)/2;
fx:=f(x);

if x<c then
   begin
   if fx<f(c) then begin b:=c; c:=x; y[3]:=y[2]; y[2]:=fx; end;
   if fx>f(c) then begin a:=x; y[1]:=fx; end;
   if fx=f(c) then begin a:=x; b:=c; c:=(x+c)/2; y[1]:=fx; y[3]:=y[2]; y[2]:=f(c); end;
end;
if x>c then
   begin
   if fx<f(c) then begin a:=c; c:=x; y[1]:=y[2]; y[2]:=fx; end;
   if fx>f(c) then begin b:=x; y[3]:=fx; end;
   if fx=f(c) then begin a:=c; b:=x; c:=(x+c)/2; y[1]:=y[2]; y[3]:=fx; y[2]:=f(c); end;
end;

if abs(b-a)<eps1 then begin Xmin:=x; fmin:=fx; stop:=true; end
Until stop=true;

Writeln('Min: (',Xmin:0:16,' , ',fmin:0:16,')');
END.
 

Но удивительное дело, не работает, зацикливается. Прричём если исследуемую ф-цию перевернуть (умножить на минус в смысле), найдёт экстремум, минимум, причём верно, а вот в оригинальном виде не хочет. У этой ф-ции локальный минимум получается на конце отрезка (я строил график) и он не может его найти. И максимум тоже не ищет. Или я ошибся принципиально в методе парабол, или метод не подходит (мне кажется сомнительным второй вариант). Можкт быть кто-то подскажет, где найти верный алгоритм метода парабол? А то я уже отчаялся его найти

[venco]

Вы не исправили вот это:

venco в сообщении #385805 писал(а):

А по алгоритму, похоже, в коде, который обновляет a b c не хватает пары else. Если первое условие (fx<f(c)) выполнено, два других варианта выполнять не надо.

Точнее исправили, но наоборот - убрали else там, где он нужен.
В зависимости от соотношения fx<>f(c) надо выполнить 3 варианта изменений a,b,c. Но у Вас после одного варианта может выполнится другой. Чтобы этого не было, надо использовать else.

[Dilettante]

Может я снова туплю, но вот так:

Используется синтаксис Pascal

 
if x<c then
   begin
   if fx<f(c) then begin b:=c; c:=x; y[3]:=y[2]; y[2]:=fx; end
      else if fx>f(c) then begin a:=x; y[1]:=fx; end
           else begin a:=x; b:=c; c:=(x+c)/2; y[1]:=fx; y[3]:=y[2]; y[2]:=f(c); end;
end;
if x>c then
   begin
   if fx<f(c) then begin a:=c; c:=x; y[1]:=y[2]; y[2]:=fx; end
      else if fx>f(c) then begin b:=x; y[3]:=fx; end
           else begin a:=c; b:=x; c:=(x+c)/2; y[1]:=y[2]; y[3]:=fx; y[2]:=f(c); end;
end;
 

Тоже не работает для точностей меньше нескольких десятых.

[venco]

Ну ещё и к сравнению x>c надо бы else добавить.
Формулы пока не смотрел.

-- Пт дек 10, 2010 11:18:08 --

Кстати, а зачем Вам массив y[]?

[Dilettante]

Массив - незачем, как я теперь понял, это я что-то напутал. Убрал его.
Дополнил сравнения ещё одним else, к сравнению x>c, опять всё то же самое, с адекватной точностью отказывается работать. Сейчас проверил, вывел на экран после каждой итерации b и a, b с каждой итерацией всё медленней и медленней уменьшается, точка а так и остаётся какой была введена. В итоге b-a зависает на значении больше эпсилон.

[venco]

А какие начальные значения используете?
И приведите текущий код.

[Dilettante]

Исходные данные:
eps = 0.01
a = 0.1
c = 0.2 (пробовал всякие)
b = 0.6

a и b по условию фиксированные, c беру наугад (только чтоб между а и b), может быть зря это делаю.

Используется синтаксис Pascal

 
program pr7;
uses
    crt;
var
eps1: real;                   
a, b, c, x, fx, t: real;
fmin: real;                   
Xmin, Xmiddle: real;          
stop: boolean;


function f(x:real):real;
begin
f:=exp((1/3)*ln(x-2*x*x+x*x*x));
end;

function midX(var a:real; var c:real; var b:real):real;
var p1, p2:real;
begin
p1:=(b-c)*(b-c)*(f(a)-f(c))-(c-a)*(c-a)*(f(b)-f(c));
p2:=(b-c)*(f(a)-f(c))+(c-a)*(f(b)-f(c));
midX:=c+0.5*(p1/p2)
end;

BEGIN
//ClrScr;

stop:=false;
Writeln('Input accuracy eps1');
Readln(eps1);
Writeln('Input limits a, c and b (a < c < b)');
Readln(a, c, b);
y[1]:=f(a);
y[2]:=f(c);
y[3]:=f(b);

Repeat
{øàã 2}
t:=midX(a, c, b);
//Writeln('t ', t);
if t <> c then x:=t;
if t = c then x:=(a+c)/2;
fx:=f(x);
 
{øàã 3}
if x<c then
   begin
   if fx<f(c) then begin b:=c; c:=x; end
      else if fx>f(c) then begin a:=x; end
           else begin a:=x; b:=c; c:=(x+c)/2; end;
end
else if x>c then
   begin
   if fx<f(c) then begin a:=c; c:=x; end
      else if fx>f(c) then begin b:=x; end
           else begin a:=c; b:=x; c:=(x+c)/2; end;
end;

{øàã 4}
if abs(b-a)<eps1 then begin Xmin:=x; fmin:=fx; stop:=true; end;
Readln;
Until stop=true;

Writeln('Min: (',Xmin:0:16,' , ',fmin:0:16,')');
END.
 

[venco]

Все сравнения fx и f(c) должны быть наоборот.

-- Пт дек 10, 2010 12:57:58 --

Это потому, что Ваш код ищет минимум, а у функции в приведённом диапазоне - максимум.

-- Пт дек 10, 2010 13:20:03 --

А вообще, в методе парабол не нужны никакие сложные проверки.
Там надо просто самую старую точку заменять новой.

[Dilettante]

Да, действительно, после замены знаков ищет максимум. Как-то я не догадался. Насчёт самой старой точки - что именно вы имеете в виду? Ещё вот такой вопрос: при оч. большой точности выдаёт ошибку деления на ноль, точки видимо слишком близкие беруться. Не знаете как решить эту проблему? Вроде бы где-то я встречал решение, но не помню точно, по этому ли методу.

Да, и ведь на границе отрезка в т. с x=0.1 у функции по идее получается минимум, его никак не определить по этому методу?

[venco]

Dilettante в сообщении #385929 писал(а):

Насчёт самой старой точки - что именно вы имеете в виду?

Имея 3 точки (a, b, c) строите новую (x) и делаете замену: a=b, b=c, c=x. Повторять до нужной точности. Всё. В этом и состоит метод парабол.

-- Пт дек 10, 2010 14:52:38 --

Dilettante в сообщении #385929 писал(а):

Да, и ведь на границе отрезка в т. с x=0.1 у функции по идее получается минимум, его никак не определить по этому методу?

Этот минимум не является параболическим.

[Dilettante]

Действительно, просто. Спасибо за обширную помощь!