Построение графиков в консольных приложениях с помощью MathGL

Введение

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

Рисунок 1 - Построено с помощью GNU Octave и gnuplot

чем в виде массива данных

Рисунок 2 - Данные использованные при построении рисунка 1

А в аэродинамических расчетах без визуального представления вообще не обойтись

 

Рисунок 3 - Турбулентное обтекание мотоцикла (стандартный расчет из OpenFOAM)

Данные, по которым строились рисунки 3, занимают на жестком диске порядка 120 МБ в сжатом виде.

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

1. Приложение обладает графическим интерфейсом и способно само визуализировать обрабатываемые данные в одном из своих окон. Такой подход, наверное, является единственно верным при разработке приложений для бизнеса, когда потребитель хочет видеть все сразу в реальном времени. В области вычислительной механике приложения, осуществляющие и расчет и визуализацию, также известны, например, это всевозможные CAE-системы (SolidWorks, ANSYS). 
2. Приложение консольной и выводит данные в файл, данные из файла обрабатываются некоторым визуализатором типа gnuplot или ParaView. Такой подход применяется, например, при работе с пакетом OpenFOAM. OpenFOAM предоставляет набор консольных решателей, а расчетные данные отрисовываются через ParaView.
3. Приложение консольное, а за визуализацию отвечает некоторая прилинкованная библиотека.

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

Библиотеки визуализации

На просторах интернета я нашел несколько библиотек, реализующих визуализацию: MathGL, gnuplot, VTK.

MathGL является разработкой нижегородским программиста Алексея Балакина. Именно об использовании MathGL я буду писать в данной статье.

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

VTK (Visualization Toolkit) - мощный инструмент визуализации, используемый как в технике, так и в медицине (см. галерею на официальном сайте). Рисунок 3 построен в пакете ParaView, использующем VTK. Об использовании VTK буду писать в следующих статьях.

Использование MathGL

Для использования данной библиотеке необходимо подключить заголовочный файл

#include <mgl2/mgl.h>

и прилинковать библиотеку libmgl.so. Команда сборки из консоли будет выглядеть примерно так

g++ -O3 main.cpp -lmgl -o plotting 

Использование MathGL очень незамысловато, все можно осуществить посредством двух классов: mglData и mglGraph. Первый является оберткой для пользовательских данных, второй отвечает за вывод графики. Далее рассмотрим использование MathGL на примере нелинейного маятника.

Нелинейный маятник

Нелинейный маятник является классической механической колебательной системой. Эта задача интегрируется в эллиптических функциях, а достаточно подробный анализ системы можно найти в книге Г.М. Заславского и Р.З. Сагдеева "Введение в нелинейную физику".

Данная система является консервативной, поэтому имеется интеграл энергии

где φ - угол отклонения, ω - угловая скорость.

Колебания маятника описываются следующей системой уравнений

Построение фазового портрета системы

Поскольку движение рассматрвиаемой системы происходит всегда на некотором заданном уровне энергии, то фазовые траектории представляют собой изолинии энергии. В качестве входных данных при построении фазового портрета будем использовать сетку по φ
double * phi;
сетку по ω
double * omega;
двумерный массив энергии
double ** E;
в котором первый индекс отвечает φ, а второй - ω, размер сетки определяется переменными
int Nphi; int Nomega;

В первую очередь необходимо сформировать объекты-контейнеры типа mglData
mglData X(Nphi), Y(Nomega), En(Nphi, Nomega);

Здесь объект X задает сетку по φ, объект Y задает сетку по ω, объект En задает распределение энергии на заданной сетке. Конструктор класса mglData принимает на вход количество значений в координатных направлениях x, y, z. Значения незаданных размеров по умолчанию равны 1. Скопировать данные в контейнеры X, Y, En, можно, например, через циклы (есть еще иные способы копирования, которые я рассматривать не буду)
for(int i = 0; i < Nphi; i++) {
X.SetVal(phi[i], i);
}

for(int j = 0; j < Nomega; j++) {
Y.SetVal(omega[j], j);
}

for(int i = 0; i < Nphi; i++) {
for(int j = 0; j < Nomega; j++) {
En.SetVal(E[i][j], i, j);
}
}

Далее необходимо определить уровни энергии, для которых будут изображать траектории системы,
mglData v(21);
min_E = En.Minimal();
max_E = En.Maximal();
v.Fill(min_E, max_E);
Методы mglData::Minimal и mglData::Maximal определяют наименьшее и наибольшее значения в контейнере mglData.

После того как все данные размещены в контенерах mglData можно приступать к настройке объекта mglGraph
mglGraph gr(0, 1200, 800);
Здесь первый параметр задает тип вывода, 0 соответствует выводу в файл, второй и третий параметры задают размер изображения в пикселях.

MathGl позволяет формировать подписи осей в долях некоторой величины. Вызов
gr.SetTicks('x', M_PI, 0, NAN, "\\pi");
указывает объекту gr выводить подписи горизонтальной оси ('x') в долях числа π. Хочу обратить внимание, что последний параметр SetTicks принимает LaTeX-овский код символа π.

Вызов
gr.SetRanges(phi[0], phi[Nphi - 1], omega[0], omega[Nomega - 1], min_E, max_E);
указывает диапазон значений для каждого координатного направления. Пятый и шестой параметры используются объектом gr для определения цветов изолиний.

Вызовы
gr.Box();
gr.Axis();
gr.Label('x', "\\varphi", 0);
gr.Label('y', "\\omega", 0);
отвечают за отрисовку рамки, вывод осей и подписей осей. Второй параметр Label может принимать LaTeX-овский код.

Вызов
gr.Cont(v, X, Y, En, "_");
формирует само изображение.

Последее что нужно сделать - это сохранить изображение в файл
gr.WritePNG("E_MathGL.png");
Результат выполнения программы показан на следующем рисунке

 Рисунок 4 - Фазовый портрет системы (изолинии энергии)

Построение простых графиков

Уравнения нелинейного маятника можно достаточно легко проинтегрировать численно. Я сделал это с помощью метода Рунге-Кутты 4 порядка точности. Эволюцию динамических переменных φ и ω во времени я изобразил в виде графиков.

Если входные данные представить в виде массива (размером nt) моментов времени

double * time;

двумерного массива (размером nt x 2), хранящего φ и ω

double ** y;

то код отрисовки будет следующим

mglData T(nt), Phi(nt), Omega(nt);

for(int i = 0; i < nt; i++) {

T.SetVal(time[i], i);

Phi.SetVal(y[i][0], i);

Omega.SetVal(y[i][1],i);

}

double min_dPhi = Omega.Minimal();

double max_dPhi = Omega.Maximal();

mglGraph gr_Phi(0, 1200, 800);

gr_Phi.Clf();

gr_Phi.SetTicks('y', M_PI, 0, NAN, "\\pi");

gr_Phi.SetRanges(0, time[nt - 1], -M_PI, M_PI);

gr_Phi.Box();

gr_Phi.Axis();

gr_Phi.Label('x', "Time", 0);

gr_Phi.Label('y', "\\varphi", 0);

gr_Phi.Plot(T, Phi);

gr_Phi.WritePNG("Phi_MathGL.png");

mglGraph gr_dPhi(0, 1200, 800);

gr_dPhi.Clf();

gr_dPhi.SetRanges(0, time[nt - 1], floor(min_dPhi), ceil(max_dPhi));

gr_dPhi.Box();

gr_dPhi.Axis();

gr_dPhi.Label('x', "Time", 0);

gr_dPhi.Label('y', "\\omega", 0);

gr_dPhi.Plot(T, Omega);

gr_dPhi.WritePNG("dPhi_MathGL.png");

Результат выполнения программы изображен на следующих рисунках 

Рисунок 5 - Отклонение маятника (слева), угловая скорость маятника (справа)

В качестве заключения

На мой вкус библиотека MathGL архиполезный и простой в использовани инструмент. Для построения графика требуется небольшое количество строк кода. Библиотека понимает LaTeX-овские коды. Особенно мне понравился вывод метод оси в долях числа π.

MathGL категорически рекомендуется к использованию. 

Исходные коды доступны по ссылке

В продолжение затронутой темы будут рассмотрены gnuplot и VTK.