Фрактал «магнитный маятник»
ОГЛАВЛЕНИЕ
Здесь мы рассмотрим реализацию фрактала «магнитный маятник»
• Скачать демонстрационный проект - 173 Кб
• Скачать исходники - 237 Кб
Введение
Статья под названием "Эксперименты над хаосом" из немецкой версии Scientific American датируется 1994 годом и среди прочего описывает модель, показывающую беспорядочное движение маятника под влиянием гравитации и трех магнитов. Была написана программа, реализующая эту модель. Программа выдает изображение фрактала с высоким разрешением
Вычисление занимает немало времени, примерно 4-5 часов, что не странно для размера изображения 1000 x 1000 пикселей с использованием относительно быстрого одноядерного процессора. Приложение не дает быстрых результатов в реальном времени. Ниже кратко описано то, что есть в статье:
Что вам понадобится:
• Интерес к теории хаоса.
• Быстрый процессор.
• Аппаратная поддержка OpenGL. (Нет, не трехмерного, были переработаны двухмерные процедуры из экранной заставки wator.)
• Достаточно времени.
Что вы получите:
• Игрушку для создания красивых картинок и для игры
• Конфигурировать модели с любым числом источников с помощью файлов INI
• Можно прерывать и продолжать вычисления
Вы можете повторно использовать код из данной статьи в своих собственных проектах. Если вам нужно это, то используемые в данном проекте вспомогательные классы могут быть весьма полезны. Есть генератор строковых меток и другие полезные вспомогательные функции, есть класс вектора фиксированного размера, использующий приемы метапрограммирования шаблонов, и есть класс для чтения файлов INI. Очень полезные вещи, но в статье не рассматриваются.
Эффект бабочки
Программа демонстрирует эффект бабочки. Для тех, кто незнаком с этим выражением, приведено краткое объяснение из Википедии:
Эффект бабочки – выражение, заключающее в себя более техническое понятие о критичной зависимости от начальных условий в теории хаоса. Малые изменения начального состояния динамической системы могут вызвать большие изменения в долгосрочном поведении системы.
Более подробно смотрите в оригинале: Эффект бабочки в Википедии.
Обзор модели
Данная статья описывает численную модель, демонстрирующую, как малые изменения в начальных условиях моделирования приводят к большим изменениям результатов. Результатом является непредсказуемость результата моделирования, так как даже малейшее изменение в условиях может сильно повлиять на результат.
Начнем с деталей. Классическая модель предполагает наличие магнитного маятника, притягиваемого тремя магнитами, каждый из которых имеет разный цвет. Магниты расположены под маятником на круге с центром в точке крепления маятника. Магниты достаточно сильные, чтобы притягивать маятник так, чтобы он не остановился в среднем положении. Следующий рисунок показывает вид модели сверху. Цветные круги обозначают магниты, пересечение в центре – точка крепления маятника. Моделирование будет вычислять траекторию, выбираемую маятником под влиянием всех трех источников с гравитацией и трением. Из-за потери энергии в результате трения маятник рано или поздно остановится над одним из магнитов. Затем начальная точка окрашивается в цвет этого же самого магнита. Выполнение этого для всех пикселей даст красивую карту, показывающую узор, состоящий из красных, зеленых и синих пикселей.
Рисунок 1 (слева): Отрисовка экспериментальной модели ; (справа): Малые изменения в начальных положениях приводят к большим изменениям в движении маятника
Результаты показаны на рисунках 1 и 3. Белое пятно в центре рисунка 3 вызвано остановкой маятника в центре. Это обусловлено параметрами, использованными в этом примере. Так как малоокрашенные изображения некрасивы, надо было добавить больше цвета. Единственная информация, доступная для каждого пикселя, - цвет магнита, верно? Поиск в интернете показывает, что большинство кодов, реализующих эту модель, останавливаются тут, но чуть больше кода может существенно улучшить результаты. Помните, для каждого пикселя вычисляется полная траектория маятника. Следовательно, перевести часть информации из траектории в цветовую информацию – хорошая идея. Самое банальное – длина траектории. Поэтому использование длины траектории для определения яркости пикселя кажется логичным шагом (это уже было сделано в оригинальной версии, опубликованной в Scientific American). Это можно сделать с помощью функций пересчета цвета, принимающих длину траектории и максимальную длину траектории в качестве параметров. Пересчет применяется к цвету источника. Вообще, приложение может использовать любую функцию, но на деле полезны следующие три:
Рисунок 2: Обзор с помощью разных цветовых функций
Как видно, результат весьма интересный, если, конечно, вы интересуетесь теорией хаоса (если же нет, то любопытно, почему вы все еще читаете?). Начальные точки, дающие более длинные траектории, показаны в более темных цветах, придавая дополнительную сложность изображению. Приложение позволяет определять пользовательские формулы для пересчета цвета. Эти формулы будут интерпретироваться с помощью muParser.
Рисунок 3 (слева направо): Структура модели; Цвет определяется индексом магнита; Цвет определяется индексом магнита и длиной траектории
Упрощения
Вся эта трехмерная физика реализована так просто, что об этом неловко говорить. Имеет место мухлеж!
Маятник является упрощенной двумерной версией, предполагающей, что сила, тянущая маятник назад к центру, подчиняется закону Гука (пропорциональна расстоянию). Это упрощение, избавляющее от необходимости вычислять углы поворота, перекрёстные произведения и все вещи, которые понадобились бы в противном случае. Реализация физически верной модели потребовала бы слишком много добавочной работы. Главной целью было получить картинку для стены, а физически верная версия преобразовалась бы в шар, а не в плоскость. Поскольку шар нельзя повесить на стену, была выбрана двухмерная версия. Разумеется, двухмерная версия верна лишь для маленьких удлинений.
Рисунок 4: Примеры кривых сила против расстояния для силы Гука и для магнитной силы
Предполагается, что магниты вызывают силу, пропорциональную обратному квадрату расстояния до маятника. В принципе, это сродни закону тяготения или закону Кулона. Все эти законы очень похожи, но здесь речь идет о (гипотетических) магнитных монополях, а не массах или зарядах. Это предположение соответствует тому, что все делают при моделировании маятника и магнитов. В реальности магниты являются диполями. Диполь вызывает силы, пропорциональные 1/r?, а не 1/r?. Расчет силы не учитывает это, хотя можно было бы смоделировать диполь с помощью двух источников-монополей. Предполагается, что маятник сделан из железа, причем не учитываются вихревые токи, которые были бы индуцированы в реальности.
Основные уравнения
Движение маятника вычисляется путем двойного интегрирования по ускорениям, действующим на маятник. Обычно бы речь шла не об ускорениях, а о силах. Согласно первому закону Ньютона, сила, необходимая для перемещения тела, равна произведению массы на ускорение. Это уравнение было решено для ускорения.
Так как начальные условия дают начальное положение и начальную скорость (считающуюся нулевой), надо лишь вычислить ускорения. Для простоты масса считается равной одной единице массы. Как правило, моделирование особо не заботится о физических единицах. Это не проблема, поскольку использование реальных единиц лишь наложило бы коэффициенты пересчета на параметры. Следующие уравнения перечисляют все ускорения, важные для моделирования:
Код для прослеживания траектории маятника
При наличии начального положения маятника и начальной скорости маятника остается лишь найти подходящую схему интегрирования и следовать по траектории маятника. Для этого моделирования был использован алгоритм интегрирования Бимена. Применение данной схемы не требует много кода, и он весьма точен. Представленный в виде псевдокода алгоритм выглядит так:
пока прослеживается траектория маятника
положение += скорость
ускорение = 0
для всех источников силы
ускорение += ускорение, вызванное источником
если (маятник рядом с источником, и скорость маленькая) то
магнит останова = индекс источника
прервать
конец ветвления if
конец цикла for
ускорение -= ускорение, вызванное трением
скорость += ускорение
длина траектории += длина(скорость)
сохранить магнит останова
сохранить длину траектории
конец цикла while
Реализация кода для прослеживания траектории, выбранной маятником, на C++ выглядит так:
for (int ct=0; ct<m_nMaxSteps && bRunning; ++ct)
{
// вычислить новое положение
pos[0] += vel[0]*dt + sqr(dt) * (2.0/3.0 *
(*acc)[0] - 1.0 / 6.0 * (*acc_p)[0]);
pos[1] += vel[1]*dt + sqr(dt) * (2.0/3.0 *
(*acc)[1] - 1.0 / 6.0 * (*acc_p)[1]);
(*acc_n) = 0.0; // сбросить ускорение
// Вычислить силу, речь идет о силах, пропорциональных
// расстоянию или обратному квадрату расстояния
for (std::size_t i=0; i<src_num; ++i)
{
const Source &src( m_vSources[i] );
r = pos - src.pos;
if (src.type==Source::EType::tpLIN)
{
//---------------------------------------
// Закон Гука: _
// _ r _
// m * a = - k * |r| * --- = -k * r
// |r|
//
(*acc_n)[0] -= src.mult * r[0];
(*acc_n)[1] -= src.mult * r[1];
}
else
{
//---------------------------------------
// Магнитные силы: _
// _ r
// m * a = k * -----
// |r?|
//
double dist( sqrt( sqr(src.pos[0] - pos[0]) +
sqr(src.pos[1] - pos[1]) +
sqr(m_fHeight) ) );
(*acc_n)[0] -= (src.mult / (dist*dist*dist)) * r[0];
(*acc_n)[1] -= (src.mult / (dist*dist*dist)) * r[1];
}
// Проверить условие прерывания
if (ct>m_nMinSteps && abs(r)<src.size && abs(vel)<m_fAbortVel)
{
bRunning = false;
stop_mag = (int)i;
break;
}
} // для всех источников
//--------------------------------------------------------------
// 3.) Трение пропорционально скорости
(*acc_n)[0] -= vel[0] * m_fFriction;
(*acc_n)[1] -= vel[1] * m_fFriction;
//--------------------------------------------------------------
// 4.) Интегрирование Бимена для скоростей
vel[0] += dt * ( 1.0/3.0 * (*acc_n)[0] + 5.0/6.0 *
(*acc)[0] - 1.0/6.0 * (*acc_p)[0] );
vel[1] += dt * ( 1.0/3.0 * (*acc_n)[1] + 5.0/6.0 *
(*acc)[1] - 1.0/6.0 * (*acc_p)[1] );
//--------------------------------------------------------------
// 5.) перевернуть буфер ускорений
tmp = acc_p;
acc_p = acc;
acc = acc_n;
acc_n = tmp;
};
Использование приложения
Графический пользовательский интерфейс
Объяснение графического пользовательского интерфейса простое. Его нет! Все параметры модели читаются из файла INI, передаваемого в качестве единственного параметра программы. Приложение открывает единственное окно. Сразу начинается расчет. Так как ждать результата расчета скучно, была встроена интерактивность. Всегда, когда вы двигаете мышью, вычисляется траектория, начинающаяся с текущего положения мыши, и рисуется в окне приложения. Это игрушка, но помогает получить представление об итоговом изображении путем изучения определенных начальных точек и влияния маленьких изменений. Если вы вышли на данную статью, ища эффект бабочки, используйте указанную функцию и наблюдайте за изменением траектории.
Чтобы узнать, что приложение вычисляет в текущий момент, нажмите правую кнопку мыши, чтобы увидеть траектории для пикселей, вычисляющихся сейчас.
Рисунок 5: Окно приложения с производящимся расчетом
Приложение использует многоядерные процессоры путем порождения одного потока вычислений на каждое ядро и привязывая поток к этому ядру. Данный подход использует Windows API, а не OpenMP, что было бы проще, но требовало бы Visual Studio 2005.
Формат файла конфигурации
Формат файла конфигурации точно такой же, как и для всех файлов INI Windows. Надо знать разделы, ключи и значения файла INI.
|
Раздел |
Ключ |
Описание |
|
[Поле] |
COLS |
Число столбцов, используемых для дискретизации поля. Этот параметр задает ширину выходного изображения. |
|
HEIGHT |
Число строк, используемых для дискретизации поля. Этот параметр задает высоту выходного изображения. |
|
|
SIM_WIDTH |
Ширина моделируемого поля в единицах длины (метр). Установлен вCOLS,если не задан. |
|
|
SIM_HEIGHT |
Высота моделируемого поля в единицах длины (метр). Установлен вROWS,если не задан. |
|
|
WIN_WIDTH |
Ширина окна вывода в пикселях. Используйте этот параметр для корректировки ширины окна вывода при расчете изображений, превышающих размеры экрана. Установлен вCOLS,если не задан. |
|
|
WIN_HEIGHT |
Высота окна вывода в пикселях. Используйте этот параметр для корректировки высоты окна вывода при расчете изображений, превышающих размеры экрана. Установлен вROWS,если не задан. |
|
|
[Моделирование] |
THREADS |
Определяет число потоков, порождаемых для расчета. Если не задано, это число равняется числу процессоров, сообщенному системой. Все потоки выполняются на разных ядрах, давая увеличение производительности в многоядерных системах (около 40% на двухъядерном процессоре). |
|
FRICTION |
Коэффициенттрения. Сила трения пропорциональна этому значению за вычетом скорости маятника. |
|
|
PEND_HEIGHT |
Высота плоскости маятника над плоскостью магнита в пикселях. |
|
|
DELTA_T |
Размер шагов интегрирования по времени в единицах времени. Чем меньше, тем точнее; это параметрhв формуле интегрирования Бимена. |
|
|
MIN_STEPS |
Минимальное число шагов, которое должна иметь каждая траектория перед проверкой условия прерывания. |
|
|
MAX_STEPS |
Максимальное число шагов, которое может иметь траектория. При достижении этого числа моделирование останавливается, даже если маятник не остановился над магнитом. |
|
|
ABORT_VEL |
Если скорость маятника падает ниже этой величины, он считается остановившимся. |
|
|
COLOR_SCHEME |
Уравнение, определяющее пересчет цвета пикселя. Пересчет цвета зависит от длины траектории и от максимальной длины траектории. Это математическое выражение должно содержать переменныеlenиmax_len(смотрите рисунок 2). Пример:1/(exp(0.000001*(len*len))) |
|
|
BATCH_MODE |
Установите его в единицу для активации пакетного режима. В пакетном режиме приложение само завершает работу после окончания расчета. Это можно использовать для создания анимаций с применением сценариев командного процессора. |
|
|
[Источник...] |
TYPE |
Тип источника. ИспользуйтеINV_SQRилиLINEAR. Сила, вызываемая источником, пропорциональна расстоянию до маятника или обратному квадрату расстояния. Первый источник аналогичен закону Гука, а второй источник аналогичен закону Кулона(хотя это не закон Кулона, потому что речь идет о магнитах, а не электрических зарядах). |
|
COLOR |
Цвет источника. |
|
|
RAD |
Положения источников представлены в полярных координатах. Это расстояние источника от центра модели. |
|
|
THETA |
Положения источников представлены в полярных координатах. Это угол источника от центра модели. |
|
|
MULT |
Коэффициент пересчета силы. |
|
|
SIZE |
Размер источника. Условия прерывания проверяются, только если маятник ближе к источнику, чем многие пиксели. |
Примеры файлов конфигурации находятся в подкаталоге данных. Если хотите создать собственные файлы конфигурации, начните с изменения одного из существующих файлов.
Галерея
Наконец, ниже представлено несколько изображений, вычисленных с помощью приложения. Изображения были получены с использованием разных наборов параметров относительно силы и количества магнитов, а также положения точки крепления маятника.
Рисунок 6: Галерея, показывающая изображения на основе разных наборов параметров