Когда-то давно я занимался решением подобной задачи: хотел построить график зависимости силы притяжения от расстояния между двумя магнитами, магниты смоделировал просто как кольца с током.
1. Для вычисления силы использовал известную формулу, сила равна градиенту потенциальной энергии. Тут я сделал небольшое упрощение, рассматривал только силу, приложенную к центру магнитов - и не рассматривал то, что магниты могут пытаться повернуть друг друга. Но это тоже легко вычисляется.
2. Энергия магнитного поля, в свою очередь, это пространственный интеграл произведения B*H = mu*B*B по dV. Пространство ограничиваем какими-то разумными пределами, я использовал ~200-1000 диаметров кольца.
3. Магнитная индукция аддитивна, то есть B в "точке" dV складывается из B1 и B2 - полей одного и другого кольца. Кольца одинаковые, поэтому достаточно предварительно вычислить поле вокруг одного кольца, а потом для взятия интеграла складывать уже заранее вычисленные поля.
4. Магнитное поле кольца осесимметрично, поэтому достаточно вычислить поле всего лишь на одной четверти плоскости, построенной на радиусе кольца и его оси.
5. Для вычисления этого поля - формула Био-Савара-Лапласа, и циркулярный интеграл по контуру кольца. Тут уже наверное никак особо не соптимизируешь, разве что считать только половину кольца.
Таким образом, по сути задача сводится к: циркулярный интеграл Био-Савара-Лапласа по половине кольца, посчитанный в каждой точке четвертьплоскости, и интегрирование этого поля по объёму с нужными смещениями и поворотами.
В трёх местах есть проблемы:
1. Поле вблизи кольца. Есть хороший метод в учебниках физики (к сожалению, не помню где) - рассматривать не бесконечно тонкое кольцо, тогда эта проблема легко решается.
2. Грануляция сетки для более быстрого и точного расчёта. Стандартный способ - около кольца бьём плоскость на много сегментов, а дальше от кольца - сегменты всё больше и больше.
3. Несовпадение сеток двух колец, если кольца повёрнуты друг относительно друга - я бы попробовал алгоритм Брезенхема для поиска ближайшего сегмента или даже интерполяции. Повторюсь, повороты я не реализовывал.
Это было сделано лет пять назад на ноутбуке с Core2Duo 8500, на Delphi с применением SSE2, в 5 потоков. В зависимости от грануляции расчёт занимал до суток. В итоге была получена обратно-кубическая зависимость, как и ожидалось.
Если требуется определить поведение системы - устремляете систему к минимуму энергии вариацией положений/поворотов магнитов, я бы действовал методом Ньютона. Благо, пересчёт системы с новыми положениями не требует перерасчёта магнитного поля кольца, а лишь требует расчёта пространственного интеграла.
Что бы я сделал сейчас:
1. перенёс векторные вычисления (циркулярный интеграл) на CUDA, пространственный интеграл - наверное, тоже.
2. проверил, правомерна ли аппроксимация магнита кольцом с током - просто разместил бы много маленьких колец внутри объёма магнита (и конечно, в одном направлении, как домены), и сравнил поле такой имитации с полем одного кольца.
3. сделал полноценные графики градиентов энергии по всем степеням свободы. Думаю, нашлось бы много интересного.
Простой
Сложный
