Метод наименьших квадратов и ОДУ

[Red_Herring]

robot80 в сообщении #1406428 писал(а):

Ну а как её покажешь? Код программы выложить?

Если будут желающие проверить, то лучше начать с того, что в точности вы программируете и в чем вы видите "не работает"

[robot80]

Хорошо!
Решал такое тестовое уравнение с задачей Коши $y' = -y\cdot\tg(x) + \frac{1}{\cos(x)}, y(\pi) = 1$. Точное его решение $y = \sin(x)-\cos(x)$.
Целевой функционал выглядел как $F(y) = \int_\pi^{1.49\pi} (y'-f)^2 dx$. Перешёл к сеточной постановке, задав сетку из 33 точек. Интеграл вычислял методом Ньютона - Котеса 4-го порядка (поэтому и 33 точки). Минимум функции многих переменных искал адаптивным методом случайного поиска (поэтому там два вложенных цикла). В качестве стартового набора значений функции брал случайный набор чисел на интервале $[-1; 1]$, кроме первой точки - она равна $y(\pi) = 1$. В самом алгоритме менялись все точки массива значений функции, начиная со второй. Программировал на PascalABS.NET. В результате получил функцию $y = \sin(x)$, которая тоже является решением ДУ, но с другим начальным условием. Вот код:

PROGRAM INVERSE2;
TYPE ARG = ARRAY [1..33] OF DOUBLE;
VAR  A,B,PI,T,HT: double;NT: INTEGER; M0: ARG; TT: TEXT;

FUNCTION NK4(H:DOUBLE; Y:ARG; N: INTEGER):DOUBLE;
VAR S:DOUBLE; I:INTEGER;
BEGIN
  S:=0.0;I:=1;
  WHILE (I <= N-4) DO
  BEGIN
    S:= S+ 7.0*(Y[I]+Y[I+4]) + 32.0*(Y[I+1] + Y[I+3]) + 12.0*Y[I+2];  
    I:= I + 4;
  END;
  NK4:= 2.0*H*S/45.0;  
END;

function MF(A, B,ALFA: DOUBLE; NT:INTEGER;F:ARG):DOUBLE;
VAR HT,PI,T,Q: DOUBLE; DF, Y: ARG;
BEGIN
PI:=3.14159265358979;
HT:= (B - A)/(NT-1);
//производная
DF[1]:= (-3.0*F[1]+4.0*F[2]-F[3])/(2.0*HT);
DF[NT]:= (F[NT-2]-4.0*F[NT-1]+3.0*F[NT])/(2.0*HT);
FOR VAR I:=2 TO NT-1 DO
BEGIN
DF[I]:= (F[I+1]-F[I-1])/(2.0*HT);
END;
// подыинтегральное выражение
FOR VAR J:=1 TO NT DO
BEGIN
  T := A + (J-1)*HT;
  Y[J] := (DF[J]+F[J] * TAN(T) - 1.0/COS(T))**2.0;  
END;
MF := NK4(HT, Y, NT);
end;
PROCEDURE MINIR2(A,B: DOUBLE; M0: ARG; NT: INTEGER);
 VAR Z1,Z2,ALF,Z0,BET,ALFA: DOUBLE; M1, M2, K, M3: ARG;  
 BEGIN
    ALF := (SQRT(5.0) + 1.0)*0.5; BET := (SQRT(5.0) - 1.0)*0.5;
    Z0 := MF(A,B,ALFA,NT,M0);     WRITELN('MU0=   ',Z0);            
    FOR VAR KK:=1 TO 400 DO 
    BEGIN
      WRITELN ('START OF STEP');
      FOR VAR I:=1 TO NT DO BEGIN K[I]:=1.0; END;
      FOR VAR I:=1 TO 4000 DO  
      BEGIN    
          FOR VAR MJ:= 2 TO NT DO //вот здесь и ниже меняются все элементы массива начиная со второго
          BEGIN
            M1[MJ]:= M0[MJ] + K[MJ]*(2.0 * RANDOM  - 1.0); 
          END;
          Z1 := MF(A,B,ALFA,NT,M1);             
          IF (Z1 <= Z0) THEN
          BEGIN
            FOR VAR MJ:= 2 TO NT DO
            BEGIN
              K[MJ] := ALF*K[MJ]; M2[MJ]:= M1[MJ] +  ALF*(M1[MJ]-M0[MJ]); 
            END;                
            Z2 := MF(A,B,ALFA,NT,M2);
            IF (Z2 < Z0) THEN 
            BEGIN 
              FOR VAR MJ:= 2 TO NT DO
              BEGIN
                K[MJ]:= ALF*K[MJ]; M0[MJ]:= M2[MJ]; 
              END;
              Z0 := Z2;
            END
            ELSE
            BEGIN 
            FOR VAR MJ:=2 TO NT DO
            BEGIN
              K[MJ] := BET*K[MJ];  M0[MJ]:= M1[MJ]; 
            END;
            Z0 := Z1; 
            END;
          END         
          ELSE
          BEGIN
            FOR VAR MJ:= 2 TO NT DO
            BEGIN
              M0[MJ]:= M0[MJ]; K[MJ] := BET*K[MJ];
            END;
            Z0 := Z0; 
          END;          
      END;
    END;
    writeln('MU_OPT =    ',Z0);
    ASSIGN(TT, 'TEXTMINI1.TXT');
    REWRITE(TT);
    FOR VAR I:=1 TO NT DO BEGIN WRITELN(TT,M0[I]);END;
    CLOSE(TT); 
    WRITELN('PRESS');    
END;
 
BEGIN
RANDOMIZE();
PI:=3.14159265358979; A:= PI; B:= 1.49*PI; NT := 33; HT := (B - A)/(NT-1); M0[1]:=1.0;
FOR VAR I :=2 TO NT DO BEGIN M0[I]:= -1.0 + 2.0*random(); END;
MINIR2(A,B,M0,NT);
READLN;
END.

[robot80]

Red_Herring, а что означают ваши слова про начальное условие:

Цитата:

Его, кстати можно и опустить

?

[Red_Herring]

robot80 в сообщении #1406463 писал(а):

ed_Herring, а что означают ваши слова про начальное условие "Его, кстати можно и опустить"

То самое и означают. Если ввести штраф за нарушение условий, то их можно и опустить при определении множества, на котором ищется минимум. В этой задаче вряд ли, но в других м.б. полезным

[robot80]

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

[vpb]

Да уж, ЧМО --- немалая область.

Я особо не специалист, но есть такая книга Ортега, Рейнболдт, Итерационные методы решения нелинейных систем уравнений со многими неизвестными. (А вдруг Вы ее не знаете ?) Хорошая книга, годная, классическая. Так вот, там первая глава чему-то такому (применение ЧМО к ОДУ и УЧП) как раз и посвящена, в качестве мотивировки. (Я, правда, эту первую главу как раз и не читал, т.к. у меня в оптимизации другие интересы).

[Red_Herring]

robot80 в сообщении #1406469 писал(а):

но почему-то это не работает...

В каком смысле "не работает"?

[robot80]

Red_Herring , вот же в моём последнем сообщении:

Цитата:

В результате получил функцию $y = \sin(x)$, которая тоже является решением ДУ, но с другим начальным условием.

. Более подробно: первое значение найденного массива равно 1 (оно и не менялось принудительно), а второе и последующие все укладываются в синусоиду. Результат этот, несмотря на "случайность" поиска, стабилен. И понять, откуда берётся эта ошибка, в чем она, я не пока не могу.

-- 23.07.2019, 09:21 --

vpb, я просмотрел заново эту книжку, вспомнил, что давно её тоже читал. И тогда и сейчас она меня не вдохновляет честно говоря. Она слишком теоретическая и абстрактная, в ней нет примеров реализации предложенной теории. А книжки по прикладным наукам (например, ЧМ) с теорией и без подтверждающих примеров это просто говорильня. Я в своё время начитался книжек Мышкиса, Федоренко, Моисеева... Основной посыл у них таккой: теория дело хорошее, но она может быть не конструктивной в принципе, либо плохо реализуема на практике.
Может быть, я и не прав, конечно, но почему-то сложилось у меня такое мнение. Это я так - без цели розжига дискуссии.

[vpb]

robot80 в сообщении #1406534 писал(а):

Она слишком теоретическая и абстрактная, в ней нет примеров реализации предложенной теории. А книжки по прикладным наукам (например, ЧМ) с теорией и без подтверждающих примеров это просто говорильня.

Может, такие книги по прикладным наукам и есть, не знаю. Но Ортега-Рейнболдт --- это не говорильня точно.

[robot80]

vpb, если речь о ЧМ, то например, Дэннис, Шнабель. "Численные методы безусловной оптимизации и решение нелинейных уравнений".
А так любая СТАРАЯ книжка по численным методам, где как правило всегда расписывается, как абстрактный алгоритм и формула реализуется на практике. С описанием подводных камней. Тогда, мне кажется, не боялись показать, что что-то не работает.

[vpb]

robot80
Дэннис-Шнабель --- тоже хорошая, годная. Я не понял, правда, привели Вы ее в качестве примера "говорильни", или наоборот "не говорильни".

Фундаментальный матан, теоретическая оптимизация, разработка конкретных численных алгоритмов, кодирование и решение частных задач --- это всё разные виды деятельности. Хотя и связанные. Трудно ожидать, чтобы они всегда были совмещены в одном человеке. На заре развития вычислительной математики такое совмещение в одном лице бывало чаще, сейчас --- реже, ибо наука стала сложнее. По моему, тут нечему удивляться.

robot80 в сообщении #1406990 писал(а):

где как правило всегда расписывается, как абстрактный алгоритм и формула реализуется на практике. С описанием подводных камней. Тогда, мне кажется, не боялись показать, что что-то не работает.

Не хочу давать советы безапелляционно, но моя (личная!) позиция такая: если не доказано надежно, теоретически, что алгоритм/программа будет работать --- можно ожидать провала в любой момент, на любом месте. Т.е. тут презумпция не та, что если не показали подводные камни, то всё будет тип-топ, а наоборот: если не доказали точно, что всё будет работать, то ничего хорошего заранее ожидать и не стоит.

[worm2]

robot80 в сообщении #1406449 писал(а):

Минимум функции многих переменных искал адаптивным методом случайного поиска

Не уверен, что этот метод эффективен для размерности 33. Удивительно, что он выдаёт хоть какое-то решение. По моему скромному мнению, тут без градиентных методов делать нечего.

[robot80]

worm2, ну что сказать, выбор мной метода случайного поиска не случаен, но это тема скорее для раздела "вопросы преподавания" или "дискуссионные", чем "помогите решить". Не знаю, нужно ли и можно ли это обсуждать здесь или нет.

Что касается моего текущего опыта применения метода случайного поиска в рамках разработки материла для предмета методы оптимизации, то он таков:
1) конечно, метод требует новых быстрых компьютеров. К нашей кафедре как раз прикреплен компьютерный класс, куда закупили новые машинки, да и студенты ходят с достаточно новыми ноутбуками. Так что можно пробовать...
2) метод вполне хорошо работает для интегральных уравнений Вольтерра и Фредгольма 2-го порядка для размерностей 33 и 53. Количество точек подбиралось так, чтобы при известном точном виде функции-решения целевой функционал давал значение порядка $10^{-6}$.
3) Для интегральных уравнений Фредгольма 1-го рода с регуляризацией Тихонова он работает удовлетворительно. В том плане, что график функции-ответа, конечно напоминает точный ответ, но не совсем. При этом увеличение размерности пространства ситуацию не делает ни лучше, ни хуже.
3) Обозначенная в стартовом топике задача существует у меня в двух реализациях. Первая продемонстрирована в виде кода выше. Это НЕРАБОЧАЯ реализация. И где ошибка, я не понимаю (я и тему поэтому завёл). Вторая реализует модификацию Red_Herring. Она работает, но не устойчиво. Причина этой неустойчивости, на мой взгляд, в необходимости вычислять производные численно. То есть для ДУ 2-го порядка такой метод вооще не будет работать.
4) Надо сказать, что если перейти от задачи Коши для ДУ 1-го порядка к эквивалентному нелинейному ИУ Вольтерра 2-го рода, то решение ищется хорошо и устойчиво для размерности 33 (производных то нет...). Минимум ищется тем же алгоритмом, что и представлен в коде.

[worm2]

Почему мне показалось, что метод случайного поиска не должен работать?
Ну вот, например (безотносительно данной задачи, просто минимизируем интеграл от квадрата разности между текущей и известной функцией, ожидая, что в результате получится эта известная функция) мы как-то свалились в "локальный минимум": функция совпадает с искомой во всех точках, кроме первой. На самом деле это не локальный минимум: в направлении изменения только первой координаты в сторону точного решения функционал будет убывать. Но в направлении изменения остальных 32-х координат функционал будет возрастать. Вот это оптимальное направление мы ищем случайным поиском. Для того, чтобы мы на него случайно наткнулись, нужно, чтобы первая координата была значительно больше остальных (чтобы её изменение скомпенсировало ухудшение в остальных 32 координатах), да ещё со знаком должны угадать. Сколько мы будем ждать, что нам случайно выпадет такой вектор?

[robot80]

Цитата:

Сколько мы будем ждать, что нам случайно выпадет такой вектор?

Не анализируя ваш текст, можно сказать сразу - ДОЛГО. Если вы посмотрите выше, то там 400 раз повторяется алгоритм из 4000 шагов, то есть 1'600'000 шагов. Я не даром сказал про быстрые компьютеры.
Я не хочу с вами спорить: метод случайного поиска ОЧЕНЬ затратный метод. И если функция унимодальна, если производную легко и дёшево вычислить, если она не породит дополнительных погрешностей, то да - надо использовать градиентные методы. Но в градиентных методах есть одна ахилесова пята - это необходимость вычислять производную Фреше. Насколько я понимаю статьи по этой проблематике, это не самая тривиальная задача. Как её вычислить, не погружаясь в функциональный анализ, не понятно.