Параллелизм и феномен ван дер Поля

Зачем нужны потоки, если есть параллелизм ВКПа? Поговорим об этом подробнее. По существу мы тем самым продолжим тему статьи[1], рассмотрев только более сложный пример, чем простые и абстрактные счетчики. Рассмотрим по ходу сначала пример, а уж потом и его реализацию на потоке.  Поехали?!

Уравнение ван дер Поля

При обсуждении статьи[2] было предложено решить уравнения ван дер Поля[3]. Дано обещание, что знакомство с  уравнением удивит . Но состояние удивления возникает, когда загадочен процесс свершения чего-то, не понятно достижение того или иного результата. Отсюда — удивление, восхищение и всякие там охи-ахи. Посмотрим, так ли это будет.

Решая уравнение ван дер Поля «в лоб», мы лишаемся шанса на успех. По подсчетам Вячеслава Петухова, предложившего задачу, грозит безрадостная перспектива «8 лет и 8 миллионов лет на расчет одной секунды процесса». Но, предположим, мы дождемся. Все же 8 лет это не миллион, хотя срок тоже приличный. Но нужно ли решение, полученное через час, да даже через минуту, когда счет идет на доли секунд? Вопрос, конечно, риторический.

Вид того как лихо и буквально мгновенно считается уравнение в SimInTech, если кого-то может ввести в состояние удивления, то только не нас. Поскольку, если есть результат, то существует алгоритм его получения, а если известен алгоритм, то его можно почти гарантированно реализовать. Таковы были соображения при решении взяться за данную задачу. Проблема виделась в одном — в алгоритме. Дальше, как говорится, дело техники. Но кто предполагал, что именно с ним и будет загвоздка? Так что, порой, не все так может сложиться просто, как это может по началу показаться.

То, что задача не простая, следует не cтолько из страшилок Вячеслава, сколько из материала книги[4], в которой подобные задачи рассмотрены с математической дотошностью. Более того, уравнение ван дер Поля или похожее на него включено в дополнительные материалы к книге. Фактически это же уравнение привел и Вячеслав, не указав, правда, источник. Но так уж случилось, что мы на него наткнулись и это дало дополнительную надежду на успех.  

Рассмотрим уравнение подробнее. Его код на внутреннем языке программирования (ЯВУ) платформы SimInTech приведен на рис. 1а. Именно он рассмотрен Вячеславом. В дополнительных материалах к упомянутой книге в проекте под именем «Локально-неуст.prt» он имеет несколько иной вид (см. рис. 1б). Но мы не будем вникать в такие тонкости. Нас интересует поведение уравнения. А в этом смысле они достойны друг друга.

В книге Скворцова Л.М. про автоматы ни слова. Следовательно, можно полагать, что алгоритмы в ней представлены обычной моделью блок-схем. А это совсем не помеха, т.к. хорошо известно как от любой блок-схемы перейти к автомату. Но это больше к слову, т.к. это знание для реализации на потоке не понадобится. Когда перейдем к платформе ВКПа, то нужно будет лишь знать, как подключить к ней поток. Но это, к счастью,  уже пройденный этап (см. упомянутую выше статью [1]).

Ну, так что — в путь?

Почему поток?

Чем интересно уравнение ван дер Поля? С точки зрения теории автоматов  — ни чем. Простой алгоритм, простые автоматы. Большое реальное время? Тоже не проблема. Формально у автомата дискретный такт может быть сколь угодно малым лишь бы не равным нулю. Реализовать это сколь угодно малое значение? Так это дело не теории, а практики.

От чего зависит реальное время вычисления уравнения ван дер Поля? При численном решении от шага интегрирования и, как следствие, от общего числа шагов алгоритма. Зная реальное время отдельного шага, умножая его на число шагов или складывая шаги, если их времена разные, получаем реальное время расчета. Чтобы его уменьшить, нужно просто уменьшить число шагов и/или их длительность. Далее опробуем и то и другое.

В простейшем случае шагу интегрирования назначается фиксированное значение. Общее число шагов тогда определяется просто — делением заданного модельного времени расчета на значение шага. Но что делать, если шаг очень маленький и их огромное число? А это, как раз, наш случай. 

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

Реальное время расчета и модельное время слабо связанные между собой вещи, т.к. реальное время зависит от среды исполнения алгоритма, а модельное время определяется алгоритмом расчета. Средой может быть что угодно — SimInTech, MATLAB,  ВКПа и другое вплоть до просто потока. Нам больше интересно именно последнее. Хотя попробуем и остальное, исключая МАТЛАБ. При этом, чем быстрее среда исполнения, тем быстрее будет получен результат. В пересчете, конечно, на один и тот же алгоритм.

Реальный шаг автоматного алгоритма в ВКПа — это дискретный такт автоматного пространства. По сравнению с тактовой частотой процессора он, прямо скажем, не очень большой, хотя от нее зависит. В потоке длительность шага определяется временем итерации программного цикла, а потому ближе к быстродействию процессора. Особенно, если поток автономен. И уж совсем рядом, если используется язык ассемблера.

Из сказанного следует, если существует проблема с временем расчета, то поток выглядит явно предпочтительнее. Удобно и то, что уравнение ван дер Поля можно рассматривать как автономное решение. Такой взгляд на него позволяет не только просто, удобно, но и эффективно реализовать его. Ниже мы на этом и сосредоточимся.

Потоки и С++

Создадим в Qt Creator простой проект типа «Приложение Qt Widgets» (ссылка на архив проекта будет приведена в конце статьи). Это будет один в один проект, рассмотренный  нами ранее в статье[5]( https://habr.com/ru/articles/821283/). Добавим в него класс уравнения ван дер Поля. В качестве его предка возьмем отработанный класс потока-счетчика. Код нового класса приведен на листинге 1, а диалог главного окна приложения с внесенными в него изменениями на  рис. 2.

Создание нового класса по существу сводится к перегрузке метода run() базового класса. Сам код класса не намного сложнее, чем решения на ЯВУ в SimInTech (см. рис. 1), т.к. реально к уравнению имеет отношение буквально несколько строк. Они приведены отдельно на листинге 2. Все настолько просто, что было бы достаточно простого Си. Но объектные технологии имеют свои преимущества. Например, мы воспользовались принципом наследования ООП.

class ThVanDerPol: public ThCounter { public:     ThVanDerPol(long *pIntMax, MainWindow *parent = nullptr);     ~ThVanDerPol();      double  pVarStep{0};     double  pVarY1{0}, pVarY2{0};     double  pVarTimeOneClock{0};     double  dY1{0}; double dY2{0};     bool    pVarIfOpt{false};     bool    bIfRealTime{false};     bool    bIfEndTime{false};      double  pVarKxY1{1};             //  Коэффициент усиления - y1     double  pVarKxY2{110};           //  Коэффициент усиления - y2     double  endtime{100};            //  Конечное время расчета     double  hmin{1e-100};            //  Минимальный шаг     double  hmax{0.01};              //  Максимальный шаг      void run();     double Integral(double dT, double dKx, double dX, double dConstant, double *dPreviousX, double *dIntegrator);     double Parabolic(double a0, double a1, double a2, double dx);     double OptStep(double dX, double dStep);      QElapsedTimer qeOneClock;     QElapsedTimer qeRealTimeCalculation;        // real time calculation     double  dRealTimeCalculation{0};     double  dStepTime{0};      int     nCycles{0};     double  dPrevX{0}; double dPrevStep{0}; double dPrevCoefStep{0};     double  dModelTime{0};     double  dPreviousX_Y1{0}; double dIntegrator_Y1{0};     double  dPreviousX_Y2{0}; double dIntegrator_Y2{0};     double  dInY1{0}; double dInY2{0}; };  ThVanDerPol::ThVanDerPol(long *pIntMax, MainWindow *parent): ThCounter(pIntMax, parent) { }  ThVanDerPol::~ThVanDerPol(void) {     bIfRun = false;     quit();         // корректный способ завершения потока     wait(300);     // ожидание завершени потока (неограниченное ожидание wait())                     // подробнее см. Боровский А. Qt4.7... стр.170 }  // double ThVanDerPol::Integral(double dT, double dKx, double dX, double dConstant, double *dPreviousX, double *dIntegrator) {      double dCurrentX = dX;                  // текущее значение входа     dCurrentX *=dKx;     double dS = (dT * (*dPreviousX + dCurrentX))/2; // площадь текущего элемента     *dIntegrator += dS;                     // общая площадь (значение интеграла)     double dY = *dIntegrator;               // "усилить" интеграл     dY += dConstant;                        // добавить константу     *dPreviousX = dCurrentX;                // сохранить текущее значение входа     return dY; }  double  ThVanDerPol::Parabolic(double a0, double a1, double a2, double dx) {     return a0 + a1*dx + a2 * qPow(dx, 2); }  double ThVanDerPol::OptStep(double dX, double dStep) {     return 1+qFabs((dX - dPrevX)/dStep); }  void ThVanDerPol::run() {     qeOneClock.restart();     qeRealTimeCalculation.restart();     while(bIfRun)  { // шаг при отключенной оптимизации         if(!pVarIfOpt&&!bIfRealTime) { dStepTime = hmax; } //         double dKx1 = pVarKxY1; double dKx2 = pVarKxY2;         if (!bIfRealTime) { dStepTime = pVarStep; dInY1 = dY1; dInY2 = dY2; }         else { dInY1 = dY1; dInY2 = dY2; } //  вычисление функции ван дер Поля         double dParbl = Parabolic(1.0, 0.0, -1.0, dY1);         double dMul = dParbl*(dY1 + dY2);         dY2 = Integral(dStepTime, dKx2, dMul, 2.0, &dPreviousX_Y2, &dIntegrator_Y2);         dY1 = Integral(dStepTime, dKx1, dY2, 2.0, &dPreviousX_Y1, &dIntegrator_Y1); // первый цикл расчета         if(nCycles==0) { dPrevX = dY2; } // ========= адаптивный алгоритм выбора шага интегрирования         if (pVarIfOpt && dStepTime) {             double dCurCoef = OptStep(dY2, dStepTime);             double dTime = (hmax / dCurCoef);             if (dTime< hmin)  dTime = hmin;             dStepTime = dTime;         }         else             if (!bIfRealTime) dStepTime = hmax; // ========= // считать модельное время функции         dModelTime += dStepTime;         if (dPrevStep != dStepTime) { pVarStep = dStepTime;  dPrevStep = dStepTime; }         dPrevX = dY2; // притормозить поток         QString qstr = pFThCounter->pVarStrSleepLotThreads;         if (qstr.length()>0) { int nSleep = qstr.toInt(); usleep(nSleep); } // шаг интегрироваия в реальном времени         if (bIfRealTime) {             dStepTime = double(qeOneClock.elapsed())/1000;         }         pVarY1 = dY1;         pVarY2 = dY2;         pVarTimeOneClock = dStepTime;         qeOneClock.restart();         nCycles++; // завершить вычисление функции         if (bIfEndTime && (dModelTime>endtime)) {             bIfRun = false; pFThCounter->DecrementActive();         } // отразить время работы в реальном времени         if (dModelTime>endtime) {             dRealTimeCalculation = qeRealTimeCalculation.elapsed()/1000;             qeRealTimeCalculation.restart();             dModelTime = 0;         }     } } 

//  вычисление функции ван дер Поля         double dParbl = Parabolic(1.0, 0.0, -1.0, dY1);         double dMul = dParbl*(dY1 + dY2);         dY2 = Integral(dStepTime, dKx2, dMul, 2.0, &dPreviousX_Y2, &dIntegrator_Y2);         dY1 = Integral(dStepTime, dKx1, dY2, 2.0, &dPreviousX_Y1, &dIntegrator_Y1); 

Таким образом, используя прототип, мы затратили минимум усилий, для создания потока, в который поместили код уравнения. Хотя, если быть точным, это тоже лишь «заготовка», которая позже будет адаптирована под требования платформы ВКПа.

В новом окне диалога сохранены все поля из предыдущего проекта, хотя это совсем не обязательно, и введены новые. Здесь можно задать минимальный и максимальный шаги интегрирования, коэффициент одного из интеграторов — Kx(y2), определить модельное время расчета — endtime, выбрать режим оптимизации шага интегрирования — переключатель adaptation. Диалог отображает текущие значения выходов интеграторов, модельное время расчета, текущее значение шага интегрирования — Step и реальное время расчета — real time. Расчет уравнения может выполняться в непрерывном режиме или до достижения модельного времени (переключатель endtime).

Алгоритм адаптивного выбора шага прост. Сначала определяется скорость изменения выбранной переменной. После этого шаг интегрирования уменьшается пропорционально увеличению скорости изменения данной переменной. Такой адаптивный алгоритм весьма прост и потому не самое лучшее решение. Например, на больших скоростях он будет явно промахиваться при переходе на меньший шаг. Но для демонстрации самого принципа адаптивной подстройки шага его вполне достаточно.

Тестирование показало, что при фиксированном шаге 1e-4(0.0001) коэффициент усиления Kx(y2) не превышает значения 1.1e4(11000) при почти мгновенном времени расчета модельного времени 3 сек (0.001 — 0.002 сек). Если на порядок увеличить коэффициент усиления, то необходимо на порядок уменьшить и шаг, а это в свою очередь на порядок увеличит время вычисления. Так, уменьшение шага до 1e-6 увеличивает коэффициент до 1.1e6, но и время возрастает до 0.15 сек. При дальнейшем увеличении на порядок шага время возрастает до более чем 1.5 сек. Для работы в реальном времени это уже немного многовато.

А какой результат мы получим, включив адаптивный режим. Здесь возникает интересная ситуация. При диапазоне шага от 0.01 до 1e-100, с одной стороны, мы можем без проблем дойти до коэффициента 1.1e5, повышая одновременно точность расчета, за условно приемлемое время расчета около 4-х секунд. Однако, при таком же коэффициенте усиления, но жестком шаге 1e-5, позволяющем его достичь, время будет в пределах всего лишь 0.015 сек. Так в чем тут будет выгода от адаптации? В повышении точности расчета? Как-то выглядит сомнительно.  

Вывод по результатам тестирования, похоже, один: читайте умные книжки с правильными алгоритмами адаптивного численного интегрирования. В частности книгу [4]. Нашему адаптивному алгоритму уравнение ван дер Поля явно не по силам. А достигнутые результаты пока могут больше расстраивать, чем удивлять.

Визуальные формы уравнения  ван дер Поля

Мы рассмотрели программные формы решения уравнения. Проект на С++ что-то считает, но качество расчета оценить невозможно, т.к. нет визуализации данных. А появление в поле значения выхода интегратора nan говорит лишь о том, что вычисление «ухнуло», но вот о том, что было до этого, сказать сложно. Что-то считалось, но как?

Современные решения тяготеют к визуальным формам. И дело не только в языках. Визуальные платформы «в коробке» содержат готовую визуализацию результатов, наглядные языки описания разного рода схем и графических объектов, библиотеки готовых блоков, гибкие средства управления вычислениями, визуальные средства отладки и т.д. и т.п.

Визуальную форму уравнения ван дер Поля в SimInTech и результаты работы демонстрирует рис. 1. Правда, простота схемы обманчива, т.к. содержит загадки. Главная из них — адаптивный алгоритм/метод численного интегрирования.

Подобная «скрытность» один из минусов визуального проектирования, т.к., если потребуется что-то подправить или в чем-то разобраться, то сделать это, как правило, будет сложно. Мы находимся именно в такой ситуации. Как «вскрыть» адаптивный алгоритм вычисления шага интегрирования? Автору это не известно. Может, Вячеслав, подскажет?

Ну, да ладно. Пока суть да дело создадим аналогичную схему в ВКПа. Все необходимые для этого блоки есть. Проблема возникает на этапе применения параметров уравнения к схеме. К тем из них, которые существенно влияют на поведение схемы относятся: 1) коэффициент усиления интегратора y2, 2) параметры дискретного шага и 3) метод интегрирования. В схеме на рис. 1 им соответствуют: 1) коэффициент усиления: mu=1E20; 2) шаг максимальный и минимальный: hmin=1E-100, hmax=0.01; 3) метод интегрирования: ARK21(Адаптивный1). 

Задать подобные значения параметров в ВКПа пока проблема. А где-то даже не имеет смысла, т.к. она заточена под проекты реального времени, где, например, столь малое дискретное время не реализуемо в принципе. Только поэтому в ВКПа есть обычный интегратор, но нет адаптивного. Просто потому, что шаг интегрирования уже меньше 1 мсек задать почти невозможно. И хотя дискретное время можно и здесь менять «на лету», но что это в данном случае даст, если требуется 1e-100 сек? Или вам известна ОС реального времени, где такой такт возможен?

Так неужели есть нечто, что нельзя реализовать в ВКПа? Может быть даже есть, но только не в этом случае. Т.к., если имеется решение в SimInTech, то нет препятствий к его реализации и в ВКПа. До этого момента это было именно так. Справедливости ради надо сказать, что возможно и обратное. Хотя, как показывает практика, здесь могут возникнуть проблемы, но это уже тема отдельного разговора.

Задача ван дер Поля примечательна тем, что в лоб решить ее невозможно, т.к. на это потребуется очень и очень много времени (расчет Вячеслава в комментариях к упомянутой статье). Тем не менее, выход есть. В его основе адаптивные методы интегрирования. Именно их работу демонстрирует график для шага на рис. 1.

Итак, у нас есть схема, аналогичная схеме на рис. 1. Поскольку пока сложно напрямую применить параметры исходной схемы, то подберем такие, которые были бы применимы к обеим схемам и сравнивать будем на этом уровне. Поколдовав с параметрами, остановимся на следующих: 1) постоянный шаг интегрирования на уровне 0.01, 2) метод интегрирования — Эйлера и 3) коэффициент усиления для интегратора y2, например, 30.

Если мы запустим проект в SimInTech с указанными параметрами, то получим «картинку», показанную на рис. 2. Результат работы аналогичной схемы в ВКПа представлен на рис. 3. В целом получилось, скажем так, очень похоже. Но, правда, по ходу схема и ее параметры были подкорректированы. Последовательность подобной подгонки демонстрирует рис. 4.

Рассмотрим, как это происходило. Рис. 6а показывает работу схемы с такими же параметрами. Это, безусловно, далеко не то, что хотелось бы видеть. Ускорим работу блоков параболической функции, умножения и сложения. Для этого переместим их в более быстрое пространство. Результат демонстрирует рис. 6b. Уже лучше. Но только дальнейшее уменьшение шага интегрирования (скорости автоматного пространства с интеграторами) с 10 мсек до 2 мсек приводит к попаданию «в яблочко». Это и демонстрирует рис. 5.  

В целом мы уже хорошо подготовились к решению поставленной задачи. Во-первых, мы создали аналог схемы SimInTech. Во-вторых, убедились, что он работоспособен, а результаты его работы сравнимы с SimInTech. Таким образом, не в полном объеме, но мы уже реализовали задание Вячеслава Петухова.

Тестируя схему в ВКПа, мы можем узнать границы, при которых начинаются отличия в работе проектов. Можем определить также значения параметров, при которых проект сохраняет свою работоспособность. Например, в аналогичной ситуации проект в SimInTech ведет себя не очень адекватно (см. времена на рис. 7). Хотя, возможно, в такой ситуации именно такое поведение и можно считать адекватным.

Информация к размышлению

Проанализируем действия, которые привели к результату сравнимому с результатами в SimInTech. Во-первых, мы ускорили работу «обвязки». Это вполне объяснимо, т.к. в реальной ситуации простые блоки будут быстрее блоков интегрирования. Эту идею просто реализовать в ВКПа, но  проблемно в SimInTech. В последней по умолчанию все блоки по скорости неотличимы друг от друга. Кем-то утверждается даже, что они «мгновенны», хотя ни чего мгновенного нет. Назовем модель программирования, в которой все операции условно одинаковы по времени исполнения, «плоской».

Но даже в рамках «плоской модели» SimInTech блоки, как оказывается, можно условно разделить на два класса — 1) «мгновенные» и 2) включающие по умолчанию задержку на один такт/шаг модельного времени. Примеры первых блоков — математические операции, вычисление функций и т.п. Пример вторых — интеграторы, задержки, применяемые для разрыва петель схем, и другие подобные им блоки.

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

Во-вторых, при наличии проблем уменьшение шага интегрирования приближает нас к правильному решению, т.к. снижается требования к полосе пропускания прибора — третьему ограничению, присущему всем реальным объектам (подробнее об ограничениях см. в [6]). Все же любая программа — это тоже реальный объект, обладающий определенными свойствами. В том числе и присущими ей ограничениями.

Поток , как база реализации математической модели

Выше мы рассмотрели модель решения уравнения ван дер Поля в форме автоматной модели ВКПа. Соответственно, существующие ограничения платформы распространяются на модели, созданные на ее базе. Здесь полная аналогия с тем, как ограничения процессора распространяются на программы, реализованные на его базе. Чем быстрее будет процессор, тем быстрее будут программы, исполняемые им, тем меньшие ограничения им будут присущи.

Математические модели — это объекты несколько иной природы, чем модели программ. Они остаются моделями даже будучи реализованными теми же программами. Проявляется это, например, в принципиальной разнице дискретного времени формальной модели уравнения ван дер Поля и дискретного времени программной модели, которая его реализует. Дискретное время программной модели ограничено возможностями аппаратной или программной платформы (того же ядра ВКПа), а дискретное время математической модели, например, только точностью переменной, которая его моделирует.

Мы упрощаем решение, если используем свойства программной модели для моделирования свойств математической модели. Именно это мы выше сделали, используя дискретное время ВКПа в качестве модели дискретного времени математической модели. Это удобно, но приводит к ограничениям, связанным с величиной минимального дискретного времени автоматного пространства ВКПа.

Но ограничение можно снять, если моделировать время математической модели. В этом случае отдельному такту программной модели назначается какой-то отрезок дискретного времени математической модели. При этом во многом снимаются ограничения на его минимальную величину. Другое дело, что при этом на расчет может потребоваться уйма реального времени. Собственно об этом и предупреждал Вячеслав.

С другой стороны, мы можем использовать более быструю программную модель. Поток, особенно автономный, в силу определенных обстоятельств значительно быстрее автоматов в ВКПа. И этим грех было бы не воспользоваться.

Только что собой представляет дискретный такт в рамках потока?

Самое простое решение — связать дискретный такт потока с длительностью итерации исполняемого им программного цикла. Далее нужно назначить ему любое время и уже в цикле считать общее модельное время. Заметим, что назначаемая/сопоставленная величина не обязательно должна иметь фиксированное значение. Используя ее значение, можно сократить время исполнения математической модели. Что собственно и делают адаптивные методы интегрирования.  

Выше, реализовав уравнение в форме отдельного проекта, мы пошли именно таким путем: использовали поток для реализации модели времени. И установили, что «жесткое» дискретное время величиной  1e-7 может быть вполне приемлемым значением для расчета (напомним, что в ВКПа это время было ограничено значением в 1e-3 сек). Далее уже нужно использовать «умные» методы управления дискретным временем модели.

Но перенесем созданное выше решение на потоке на платформу ВКПа и посмотрим, что из этого получится. Сделав это, мы достигнем минимум двух целей. Во-первых, решим проблему реализации функций, которые в среде ВКПа считались бы очень медленно. Во-вторых, реализованные подобным образом функции получат доступ к возможностям платформы. Например, для той же визуализации результатов. Последнее часто очень важно для оценки качества решения. В первом — простом проекте Qt, об этом можно говорить лишь условно.       

Реализация уравнения на потоке в ВКПа

Перетащить поток в ВКПа дело пройденное. Подробно эта процедура описана в рамках двух предыдущих статей автора, посвященных изучению свойств потоков (см. [1], [5]). Поэтому посмотрим сразу на результат.

А результат интересный. На усеченном варианте параметров (они приведены выше) мы получили картину, которую демонстрирует рис. 8. Она явно отличается от результата на  рис. 6. Как такое получилось, если до этого все уже было отработано в рамках простейшего проекта?  

Чтобы разобраться, попробуем замедлить поток. Это можно сделать, используя режим синхронизации с автоматным пространством ВКПа, к которому привязан поток. Сразу получаем картинку (рис. 9), сравнимую с рис. 6. Т.е. с одними и теми же параметрами и на примерно той же скорости работы мы получаем идентичные результаты.

Но что это за «загогулины» на фазовом портрете, которых нет на рис. 6? Все просто. Они отражают процессы, которые видны на импульсном графике (синий цвет на рис. 9). Такая генерация возникает, как показали эксперименты, когда не хватает точности. Повысим ее, уменьшив шаг на порядок. В результате получим картинку (рис. 10), которая уже сравнима с рис. 6:

В результате повышения точности ушли процессы генерации, пропали и «загогулины». Правда увеличилось время расчета. Но это пока еще какие-то жалкие 0.003 сек.

Вид полученной картинки дает надежду, что перенос потока не повлиял на сам расчет. Но это мы убедились на медленной скорости. Но все же что творится реально на максимальной скорости?

Осциллограф ВКПа, работая на меньшей скорости, чем поток, «хватает» мгновенные значения и потому формирует искаженный вид фазового портрета. Но у него есть режим отражения данных точками, не соединяя их между собой. Переведя визуализацию в этот режим, получим картину, которая показана на рис. 11. Теперь, соединив мысленно точки на рис. 11, сравним этот образ с рис. 10. В целом это будет одна и та же картина.

Попробуем выжать максимум из режима работы при постоянном шаге интегрирования. На рис. 11 мы задали шаг — 1e-8 при коэффициенте усиления — 10000. Получили в результате время расчета более чем 350 сек (пришлось долго ждать) и следующий фазовый портрет (рис. 12). По форме похоже, но по времени, конечно, многовато.

Давайте включим режим оптимизации при том же коэффициенте усиления. Результаты показаны на рис. 13. Время расчета — чуть более 8 сек. Ускорение более чем в 40 раз! Это уже кое-что. Но если мы увеличим коэффициент на порядок, т.е. до 100000, то время возрастет до немногим более 80 сек. Но даже это все же не 350 сек, полученных ранее при фиксированном шаге и коэффициенте усиления на порядок меньше.

Так мы фактически выжали все что можно из созданных решений.  В результате коэффициент усиления  при скорости расчета примерно 0.8 сек возрос с 30 до 1000. Далее, как показывают эксперименты, каждое увеличение коэффициента на порядок на порядок же увеличивает время расчета.

Можно ли доверять SimInTech?

Мне сложно быть объективным, но мне почему-то больше нравятся решения без всякой оптимизации. Я такому решению доверяю больше. По мне, чем проще алгоритм, тем больше ему доверия. Оптимизация — вещь, которая в разных ситуациях не только может работать по-разному, но может выдать и разные результаты. Попробую пояснить.

Возьмем сначала SimInTech и самый простой вариант с ограниченными параметрами (их значения приведены в начале статьи). Будем менять только метод интегрирования, загнав проект в область неустойчивого решения,  уменьшив шаг интегрирования. В нашем случае установим его равным 0.0147 сек. Метод интегрирования напомним — Эйлера. В этом случае возникает генерация, что демонстрирует рис. 14.

А если выбрать адаптивный метод — ARK21(Адаптивный 1)?  Ему установим диапазон шага: hmin = 1e-100 и hmax = 0.001. Получим результат, показанный на рис. 15. Генерация пропала.

Для чистоты эксперимента выберем другой метод, но тоже Эйлера,  — неявный Эйлера.  Получим результат, показанный на рис. 16. Генерации  тоже нет, но фазовый портрет выглядит достаточно «коряво».

А, ведь, мы еще не обсудили моменты появления импульса сигнала с выхода интегратора y2.

А что ВКПа? Ее возможности демонстрирует рис. 17. Дискретное время — 0.012 сек (пометка в левом нижнем углу). Одинаковыми цифрами отмечены соответствующие друг другу процессы на диаграммах. В ВКПа получилось, как минимум, красивее. Не правда ли? Ну и, как я надеюсь, точнее. Кто-то, а я знаю этого человека, с этим, скорее всего, не согласится.

Выводы

Не хочется признавать поражение, но частично — признаю. Но с чистой совестью. Так как дело совсем не в автоматах. Нужен адаптивный алгоритм. Какой он должен быть, думаю, лучше всего знает автор книги [4]  Скворцов Л.М. Знают его, надеюсь, и разработчики SimInTech. Поэтому «мячик» перебрасываю на их строну. Точнее, к Петухову Вячеславу, как инициатору темы.   

Я же взялся за задачу, которая на текущий момент за пределами моей компетенции. А если быть совсем точным, то за пределами моих интересов. Лично я за то, чтобы каждый занимался своим делом. Только тогда может получится что-то стоящее. Я буду заниматься автоматами, кто-то — математикой. Кесарю, как говорится, кесарево…

Безусловно, можно закопаться в теорию. Она есть и, отталкиваясь от нее, можно создать адаптивный алгоритм. Мы все же не на не паханном поле? Но пока я не вижу смысла тратить на это время. Своих целей я достиг. А тут нет даже элементарной мотивации. Хотя бы той, которую не так давно озвучил Вячеслав Петухов. Я даже выполнил ее требование — не просто упомянул, а написал очередной трактакт  про SimInTech. И что? Где эта самая «мотивация»? 🙂

Ну, а если серьезно. Я в очередной раз убедился, что в нашем промышленном программировании сложно, если вообще на текущий момент возможно, найти истину. Каждый считает себя истинной в последней инстанции и разводит «пальцы веером». Но где тут истина, если смотреть на приведенные выше диаграммы сигналов?

Я не тонкий знаток ни SimInTech, ни MATLAB, ни LabVIEW, ни, прости Господи, выскочившей вдруг ниоткуда, как черт из табакерки,  «Российской платформы математических вычислений и динамического моделирования» Engee. Где она таилась, если до недавнего времени о ней я  слыхом не слыхивал? А потому, имея накопленный опыт работы на других платформах, я не спешил бы «задрав штаны» переходить на SimInTech или ту же Engee. Только если уж обстоятельства заставят. А они могут заставить…

Сетовать, что кто-то где-то не следует вашим призывам, наклеивая по ходу еще и какие-то ярлыки, не стоит. Разбрасываться опытом — последнее дело. Особенно инженеру. Это торгашу по фиг чем торговать… будь это Engee, SimInTech или что-то другое. Для «заманухи» нужны более серьезные «ваши доказательства». Хотя в наше меркантильное время сторонников часто элементарно зазывают «баблом», буквально перекупая/переманивая у конкурентов…

Да, чуть не упустил. Какая связь между параллелизмом и уравнением ван дер Поля? Может показаться, что ни какой. Но посмотрим внимательнее. Это система уравнений. И соединены эти уравнения перекрестными обратными связями: y1 используется для вычисления y2 и наоборот. Помните RS-триггер, который легко входит в режим генерации? Может, именно  здесь «собака порылась», если посмотреть на генерацию на рис. 14 и рис. 17?

А если еще вспомнить уравнения аттракторов, которые мы рассматривали в цикле статей [7, 8, 9]? Структурно они совпадают с уравнением ван дер Поля. А вы — RS-триггер, RS-триггер… Опять этот триггер?! Это ж элементарно! Но элементарно ли? Вот генерит уже уравнение ван дер Поля. Неспроста это, однако. Мгновенные вещи не входят в режим генерации. Это о чем-то кому-то говорит?

А теперь — там-тара-рам! — «вишенка на торте»! Все что до этого было написано, возможно, просто бла-бла-бла!  Если представить уравнение ван дер Поля, как «черный ящик», то результат можно получить без адаптивных методов интегрирования. Вам три секунды на размышление! Раз, два, три. Придумали? Нет? Рассказываю. Берем проект с минимальными параметрами (они выше приведены), прикручиваем к выходу интегратора y2 усилитель и получаем по внешнему виду нужный «портрет». Просто, гибко…

А Вы обратили внимание, что по внешнему виду он по ходу статьи фактически один и тот же? Отличие только в значении y2. Чем больше коэффициент усиления интегратора, тем оно больше. Поэтому напрашивается  чисто инженерное решение:  получаем «портрет», который затем усиливаем по одной координате. Рис. 18 демонстрирует работу двух проектов, где верхний классический, нижний — с усилителем.

Ну, почти одно и то же? Как, Вячеслав, подойдет? 🙂

Архивы проектов:

1. Простой проект: QVanDerPol.zip — https://cloud.mail.ru/public/tGjg/gw15odLZK

2. Проект ВКПа: FsaVanderPol.zip  — https://cloud.mail.ru/public/5Fwt/KK6y9GURo

Литература

1. Все секреты многопоточности. https://habr.com/ru/articles/818903/

2. Закоулками мечты. Часть 2. https://habr.com/ru/articles/842974/

3. Кузнецов А.П., Селиверстова Е.С., Трубецков Д.И., Тюрюкина Л. В. Феномен уравнения ван дер Поля. https://sgtnd.narod.ru/papers/2014PND3.pdf

4. Скворцов Л.М. Численное решение обыкновенных дифференциальных и дифференциально-алгебраических уравнений. 2-e изд., испр. и доп. — М: ДМК Пресс, 2022. https://vk.com/wall-78467318_781

5. Засады многопоточности. https://habr.com/ru/articles/821283/

6. Глушков В.М. Введение в кибернетику. Изд-во АН Украинской ССР. К.:  1964. — 324с.

7. О программных ошибках на примере MATLAB и SimInTech. https://habr.com/ru/articles/703244/

8. Параллелизм истинный и мнимый или… и ты туда же, Рикитаке. https://habr.com/ru/articles/706192/

9. Цена ошибки. https://habr.com/ru/articles/709358/