VPF::[Матан, Delphi] Метод наискорейшего спуска

ne_tru_e

Дата 27.5.2008, 14:49 (ссылка)

(нет голосов)

Загрузка ...

Приношу жертвы MS

Профиль
Группа: Участник
Сообщений: 94
Регистрация: 20.1.2008

Репутация: нет
Всего: 4

Вопрос тем, кто разбирается в вышке.
На сайте http://www.kti.ru/data/84/mp_nsp.html есть программа на Pascal, по нужной мне теме.
Я пробовал её запускать как на Turbo Pascal, так и на Delphi, результат один - floating point overflow, то есть результаты вычислений явно выходят за границы возможных значений.
Хотя на сайте результаты есть.
Вот что выводит у меня:

--***----- ТЕКУЩИЕ ЗНАЧЕНИЯ-----***--
abs>eps
13,9512812401139
abs>eps
13,9512812401139
x[1]=3,232205
x[2]=0,023727
x[3]=-5,166087
z=13,951281
abs>eps
0,880716055388026
abs>eps
0,880716055388026
x[1]=1,189384
x[2]=2,747488
x[3]=-4,558104
z=0,880716

В Delphi все вещественные типы я поменял на Extended, все равно не помогло. Однако результаты с сайта откуда-то же могли взяться?
Я не знаю, может там намеренно сделали ошибку.
Помогите пожалуйста исправить программу на Pascal. Я тоже работаю в этом направлении. Теория есть на сайте по ссылке внизу.
Программу на Си ещё не проверял.

volvo877

Дата 27.5.2008, 16:46 (ссылка)

(нет голосов)

Загрузка ...

Эксперт

Профиль
Группа: Комодератор
Сообщений: 2073
Регистрация: 15.11.2004

Репутация: 34
Всего: 116

Цитата(ne_tru_e @ 27.5.2008, 14:49

)

может там намеренно сделали ошибку.

Использовали переводчик с С на Паскаль... У них - переводчик глючит... Вот результат работы моего транслятора с небольшой ручной правкой:

Код


type
  float = double;
const
  k: integer = 0;

function f(a: array of float): float;
begin
  inc(k);
  f := sqr(a[0]-1)+sqr(a[1]-3)+4*sqr(a[2]+5);
end;

function sgn(a, b: float): integer;
begin
  if (a - b) >= 0 then sgn := 1
  else sgn := -1;
end;


const
  max_n = 75;
var
  x, y, g, d, l, ff: array[0 .. pred(max_n)] of float;
  n: integer;


function pr(x: array of float): float;
var
  p: float;
  i: integer;
begin
  p := 0;
  g[0]:=2*(x[0]-1);
  g[1]:=2*(x[1]-3);
  g[2]:=8*(x[2]+5);
  for i := 0 to pred(n) do p := p + sqr(g[i]);
  pr := sqrt(p);
end;

var
  pf: text;
  h, eps, z, g0, g2, dn: float;
  vrem, v: float;
  i, j, m, s1, s2, s3: integer;
begin
  assign(pf, 'result.dan'); rewrite(pf);
  write('Введите число переменных:'); readln(n);
  writeln(pf,'Введите число переменных:',n);
  writeln('Введите начальную точку:');
  writeln(pf,'Введите начальную точку:');
  for i := 0 to pred(n) do begin
    read(x[i]);
    write(pf, ' ', x[i], ' ');
  end;

  writeln;
  write('Введите шаг h:'); readln(h);
  writeln(pf);
  writeln(pf,'Введите шаг h:', h);
  write('Введите точность eps:'); readln(eps);
  writeln(pf,'Введите точность eps:', eps);
  writeln('----***----Текущие значения----***----');
  writeln(pf,'----***----Текущие значения----***----');

  for i := 0 to pred(n) do begin
    y[i] := x[i];
  end;
  z := f(x);
  g0 := pr(x);

  while g0 > eps do begin
    for i := 0 to pred(n) do d[i] := -g[i] / g0;
    l[0] := 0;
    ff[0] := z;
    repeat
      m := 0; l[2] := h;
      for i := 0 to pred(n) do x[i] := y[i]+l[2]*d[i];
      z := f(x);
      ff[2] := z;
      g0 := pr(x);
      g2 := 0;
      for i := 0 to pred(n) do g2 := g2 + g[i]*d[i];
      if (ff[2]>=ff[0]) or (g2>=0) then m := 0
      else begin
        { Удвоить длину шага, чтобы накрыть min }
        h := h * 2;
        m := 1;
      end;
    until m = 0;

    l[1] := h/2;
    for i := 0 to pred(n) do x[i] := y[i]+l[1]*d[i];
    z := f(x);
    ff[1] := z;

    { Выполнить первую квадратичную интерполяцию }
    l[3]:=h*(ff[1]-0.75*ff[0]-0.25*ff[2]);
    l[3]:=h*(ff[1]-0.75*ff[0]-0.25*ff[2]);
    l[3]:= l[3] / (2*ff[1]-ff[0]-ff[2]);
    if l[3]<0 then writeln('Âíèìàíèå!');
    for i := 0 to pred(n) do x[i] := y[i]+l[3]*d[i];
    z := f(x);
    ff[3] := z;
    { Имеем 4 значения L и 4 значения функции, упорядочим их в порядке убывания}
    repeat
      m := 0;
      for i := 0 to pred(n) do
        for j := i + 1 to n do
          if ff[i] > ff[j] then begin
            vrem:=l[i]; l[i]:=l[j]; l[j]:=vrem;
            vrem:=ff[i]; ff[i]:=ff[j]; ff[j]:=vrem;
          end;

      { Закончить поиск в данном направлении если точность достигнута }
      if abs(l[0]-l[1])<eps*50 then m := 0
      else begin
        s1:=sgn(l[1],l[0]);
        s2:=sgn(l[2],l[0]);
        s3:=sgn(l[3],l[0]);
        if (s1 = s2) and (s1=-s3) then begin
          l[2]:=l[3]; ff[2]:=ff[3];
        end;
        dn:=(l[1]-l[2])*ff[0]+(l[2]-l[0])*ff[2]+(l[0]-l[1])*ff[1];
        v:=(ff[0]-ff[1])/(2*dn);
        v := v * (l[1]-l[2])*(l[2]-l[0]);
        l[3]:=(l[0]+l[1])/2+v;
        for i := 0 to pred(n) do x[i]:=y[i]+l[3]*d[i];
        z:=f(x);
        ff[3]:=z;
        m:=1;
      end
    until m = 0;

    for i := 0 to pred(n) do begin
      x[i]:=y[i]+l[0]*d[i];
      y[i]:=x[i];
      write('X[', i+1,']=', x[i], ' ');
      write(pf,'X[', i+1,']=', x[i], ' ');
    end;
    z:=f(x);
    g0:=pr(x);
    writeln; writeln(' z=', z);
    writeln(pf); writeln(pf, ' z=', z);
    h := h / 2;
  end;

  for i := 0 to pred(n) do begin
    writeln('X[',i+1,']=',x[i]);
    writeln(pf,'X[',i+1,']=',x[i]);
  end;

  writeln(' Минимум функции F( X1,X2,...,Xn)=', z);
  writeln(pf,' Минимум функции F( X1,X2,...,Xn)=', z);
  close(pf);
end.

Результат работы:

Код


Введите число переменных:3
Введите начальную точку:
4 -1 2

Введите шаг h:4
Введите точность eps:0.00001
----***----Текущие значения----***----
X[1]= 3.232204998418222E+000 X[2]= 2.372666877570376E-002 X[3]=-5.16608668142992
8E+000
 z= 1.395128124011389E+001
X[1]= 1.189383908236841E+000 X[2]= 2.747488122350878E+000 X[3]=-4.55810421411404
0E+000
 z= 8.807160553880276E-001
X[1]= 1.140914568862085E+000 X[2]= 2.812113908183887E+000 X[3]=-5.01048471494509
6E+000
 z= 5.559781620543981E-002
X[1]= 1.011955421566354E+000 X[2]= 2.984059437911528E+000 X[3]=-4.97210401634517
5E+000
 z= 3.509777240806756E-003
X[1]= 1.008895650592871E+000 X[2]= 2.988139132542839E+000 X[3]=-5.00066187876435
1E+000
 z= 2.215651103015730E-004
X[1]= 1.000754721486953E+000 X[2]= 2.998993704684063E+000 X[3]=-4.99823898319711
2E+000
 z= 1.398695550596274E-005
X[1]= 1.000561564358530E+000 X[2]= 2.999251247521960E+000 X[3]=-5.00004178306239
0E+000
 z= 8.829680993531402E-007
X[1]= 1.000047644034941E+000 X[2]= 2.999936474620078E+000 X[3]=-4.99988883058513
5E+000
 z= 5.573998316516983E-008
X[1]= 1.000035450417647E+000 X[2]= 2.999952732776471E+000 X[3]=-5.00000263767988
4E+000
 z= 3.518751952197949E-009
X[1]= 1.000003007671180E+000 X[2]= 2.999995989771760E+000 X[3]=-4.99999298210057
9E+000
 z= 2.221316655759078E-010
X[1]= 1.000002237912880E+000 X[2]= 2.999997016116160E+000 X[3]=-5.00000016651137
4E+000
 z= 1.402272098465330E-011
X[1]= 1.000002237912880E+000
X[2]= 2.999997016116160E+000
X[3]=-5.000000166511374E+000
 Минимум функции F( X1,X2,...,Xn)= 1.402272098465330E-011

ne_tru_e

Дата 27.5.2008, 20:11 (ссылка)

(голосов:1)

Загрузка ...

Приношу жертвы MS

Профиль
Группа: Участник
Сообщений: 94
Регистрация: 20.1.2008

Репутация: нет
Всего: 4

Спасибо огромное за помощь.
Мог бы конечно и сам догадаться.

furystorm

Дата 2.6.2008, 12:41 (ссылка)

(нет голосов)

Загрузка ...

Новичок

Профиль
Группа: Участник
Сообщений: 2
Регистрация: 2.6.2008

Репутация: нет
Всего: нет

уважаемые эксперты, можете пожалуйства переделать этот код под visual basic ? smile

Это сообщение отредактировал(а) furystorm - 2.6.2008, 12:41

0 Пользователей читают эту тему (0 Гостей и 0 Скрытых Пользователей)
0 Пользователей:
« Предыдущая тема \| Центр помощи \| Следующая тема »