Разное

«Если у тебя есть молот, все проблемы становятся гвоздями».

Сегодня мы рассмотрим игрушечный пример с динамической системой второго порядка и расскажем, как можно использовать свободные библиотеки для анализа подобных систем. Эта заметка, надеюсь, поможет трудолюбивым исследователям в разработке собственных инструментов прототипирования и визуализации, каковую, если речь идет об отрисовке математических абстракций, иногда называют «когнитивной». Появление этой заметки, хотя она имеет мало отношения к САПР, не должно смущать читателя. Это своеобразная дань неповторимой Горьковской школе, воздухом которой я дышал весь свой университетский срок. Моя капля пота на общий алтарь. А для чего еще нужен блог, если не для подобного рода рефлексии, особенно отвлеченной, когда в процессе ты любуешься красотами Южного Урала на озере Акакуль ;)

За основу был взят набор компонент Analysis Situs, который я все чаще использую вместо привычного Draw. Дело в том, что средство прототипирования должно «лежать в руке», и в этом отношении Draw угрюмо курит в сторонке. Если тебе нужно убрать с дороги упавшее дерево, то целесообразно использовать бензопилу, а не старый-добрый унылый топор, хотя и без него бывает не обойтись. Уверен, аналогия понятна, хотя выбор инструмента не в последнюю очередь — дело вкуса.

Рассмотрим математическую модель системы фазовой автоподстройки частоты. В безразмерном виде она задается системой обыкновенных дифференциальных уравнений второго порядка. Параметры системы положительны, а интегральные кривые располагаются на фазовом цилиндре. Опираясь на результаты уже проведенного качественного исследования, мы построим интерактивное приложение для экспериментирования с конкретной динамической системой. В учебно-методических целях можно пойти дальше и попробовать создать приложение для численного исследования «произвольных» динамических систем на два- или три-многообразиях. Нам же достаточно «поиграться» с простым примером, чтобы продемонстрировать как 3D технологии помогают в исследовании динамических процессов. Итак, дана система:

Используя связку OpenCascade и VTK, не составляет труда визуализировать фазовый цилиндр, по которому «бегают» траектории. Для этого нужно создать цилиндрическую грань OpenCascade, например, вот так:

Handle(Geom_CylindricalSurface)
  surf = new Geom_CylindricalSurface( gp_Ax3(), 10.0 );
  
Handle(Geom_RectangularTrimmedSurface)
  trimmedSurf = new Geom_RectangularTrimmedSurface(surf, -50.0, 50.0, false);
  
TopoDS_Shape shape = BRepBuilderAPI_MakeFace( trimmedSurf, Precision::Confusion() );

После визуализации средствами VTK мы имеем обыкновенный трехмерный просмотрщик, позволяющий рассматривать фазовый цилиндр с любого ракурса.

Удобство интерактивной среды состоит в том, что пользователь может произвольно выбирать параметры исследуемой системы и «выпускать траектории» путем интерактивного клика на цилиндре. Отзывчивый интерфейс специализированного приложения дает оператору возможность поставить серию вычислительных экспериментов совершенно не напрягаясь, то есть, грубо говоря, экономить мозги для действительно важного дела — теоретического качественного исследования и ОБЪЯСНЕНИЯ полученных результатов. Выбор же начальных условий (клик на цилиндре) реализуется в два этапа:

1. Нахождение пересечения луча, выпущенного из камеры, с полигональным представлением цилиндра.
2. Уточнение точки клика путем проекции точки пересечения на цилиндр.

На первом этапе используется класс vtkWorldPicker библиотеки VTK. Второй этап выполняется средствами OpenCascade следующим образом:

// Picked point (comes from VTK picker)
gp_Pnt P(x, y, z);
  
// Get stored cylindrical surface
Handle(Geom_Surface) surf = ...;
  
// Project point on surface
ShapeAnalysis_Surface sas(surf);
gp_Pnt2d UV = sas.ValueOfUV(P, 1.0e-6);
//
if ( sas.Gap() > 0.1 )
{
  std::cout << "Picked point is too far from the phase cylinder" << std::endl;
}

Сами интегральные кривые находятся общеизвестным методом Рунге-Кутта, который может запрограммировать любой студент за полчаса. Отрисовка фазовых траекторий на цилиндре дает качественный портрет движения системы. «Наматывание» траекторий на поясок цилиндра дает возможность исследователю сделать заключение о существовании предельного цикла второго рода (то есть вращательного движения). На рисунке ниже зона, где существует этот предельный цикл, показана стрелкой. Синие точки отвечают начальным условиям, выбранным пользователем интерактивно.

Любопытно отметить вот какой факт. Интегрирование методом Рунге-Кутта дает траекторию как последовательность точек в двумерном фазовом пространстве. Чтобы визуализировать ее в 3D (как это показано на рисунке), сам собой напрашивается следующий простой подход:

1. Строим полигональную линию в 2D (области определения цилиндра).
2. Вычисляем 3D-траекторию в силу цилиндра. Формула цилиндрической поверхности нам известна, поэтому проблем быть не должно.

Следует сразу предостеречь читателя от использования аппроксиматоров OpenCascade для решения этой задачи, хотя идея полиномиального приближения CONS выглядит изящно. Работать это могло бы так. Заводим адаптер вида Adaptor3d_CurveOnSurface, позволяющий вычислять точки кривой на поверхности, и аппроксимируем эту кривую явно заданным сплайном (пакет Approx). Увы, практика показывает, что сильно закрученные кривые могут плохо аппроксимироваться.

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

Handle(Geom_BSplineCurve) PolylineAsSpline(const std::vector<gp_XYZ>& trace)
{
  // Initialize properties for spline trajectories
  TColgp_Array1OfPnt poles( 1, (int) trace.size() );
  //
  for ( int k = 0; k < trace.size(); ++k )
  {
    poles(k + 1) = gp_Pnt( trace[k] );
  }
  
  const int n = poles.Upper() - 1;
  const int p = 1;
  const int m = n + p + 1;
  const int k = m + 1 - (p + 1)*2;
  
  // Knots
  TColStd_Array1OfReal knots(1, k + 2);
  knots(1) = 0;
  //
  for ( int j = 2; j <= k + 1; ++j )
    knots(j) = knots(j-1) + 1.0 / (k + 1);
  //
  knots(k + 2) = 1;
  
  // Multiplicities
  TColStd_Array1OfInteger mults(1, k + 2);
  mults(1) = 2;
  //
  for ( int j = 2; j <= k + 1; ++j )
    mults(j) = 1;
  //
  mults(k + 2) = 2;
  
  return new Geom_BSplineCurve(poles, knots, mults, 1);
}

Народ, занимающийся качественным исследованием динамических систем, знает, что такое метод точечных отображений. Идея состоит в том, чтобы исследовать поведение фазовой траектории не на всем исходном пространстве, а на некоторой редуцированной области, скажем, прямой или плоскости. Для фазовой траектории выписывается функция последования, которая может иметь или не иметь неподвижную точку. Существование неподвижной точки означает предельный цикл. Точки последования фазовой траектории тоже удобно показывать в трехмерном вьювере.

Динамика точек последования удобно визуализируется так называемой диаграммой Ламерея:

Подытоживая эти бесхитростные упражнения, мне хотелось бы сослаться на известного программиста и исследователя Брукса (Brooks, F.P. 1996. Toolsmith II. Communications of the ACM 39, 3, 61–68.). Рассуждая о том, как вычислительная машина помогает ученому в получении нового знания, он вводит понятие «амплификатора интеллекта» (intelligence amplification). Позже на него будет ссылаться автор библиотеки VTK Шродер (Schroeder, W., Martin, K., Lorensen, B. 2006. Visualization Toolkit: An Object-Oriented Approach to 3D Graphics. Kitware, Colombia.), подчеркивая, что задача научной визуализации как раз и состоит в том, чтобы усилить когнитивные способности исследователя. Эта же мысль является ключевой в монографии нашего соотечественника Зенкина, который занимается «когнитивной компьютерной графикой» (Zenkin, A. 1991. Cognitive Computer Graphics.). Штука состоит в том, чтобы разгрузить мозги ученого от рутины и дать ему надежный ассистирующий инструмент для проведения вычислительных экспериментов. Да, для построения подобных систем придется потратить некоторое время, но оно многократно окупится в процессе самой исследовательской деятельности. Программные продукты такого рода следует воспринимать как своеобразный испытательный стенд и не относиться к ним легкомысленно. Можно даже сформулировать правило:

ОДНА ПРОБЛЕМА = ОДИН ИСПЫТАТЕЛЬНЫЙ СТЕНД.

Это правило, очевидно, не гарантирует успеха, но позволяет грамотно «поставить» исследовательский процесс. Сделать его осмысленным и методичным.