Я боролся со следующим кодом. Это реализация F # алгоритма Форвард-Эйлера , используемая для моделирования звезд, движущихся в гравитационном поле.
let force (b1:Body) (b2:Body) =
let r = (b2.Position - b1.Position)
let rm = (float32)r.MagnitudeSquared + softeningLengthSquared
if (b1 = b2) then
VectorFloat.Zero
else
r * (b1.Mass * b2.Mass) / (Math.Sqrt((float)rm) * (float)rm)
member this.Integrate(dT, (bodies:Body[])) =
for i = 0 to bodies.Length - 1 do
for j = (i + 1) to bodies.Length - 1 do
let f = force bodies.[i] bodies.[j]
bodies.[i].Acceleration <- bodies.[i].Acceleration + (f / bodies.[i].Mass)
bodies.[j].Acceleration <- bodies.[j].Acceleration - (f / bodies.[j].Mass)
bodies.[i].Position <- bodies.[i].Position + bodies.[i].Velocity * dT
bodies.[i].Velocity <- bodies.[i].Velocity + bodies.[i].Acceleration * dT
Хотя это работает, это не совсем "функционально". Он также страдает от ужасной производительности, он в 2,5 раза медленнее, чем эквивалентный код C #. body - это массив структур типа Body.
С чем я борюсь, так это с тем, что force () является дорогой функцией, поэтому обычно вы рассчитываете ее один раз для каждой пары и полагаетесь на то, что Fij = -Fji. Но это действительно портит любой разворачивающийся цикл и т. Д.
Предложения с благодарностью приняты! Нет, это не домашнее задание ...
Спасибо
Ade
ОБНОВЛЕНО: Для пояснения Body и VectorFloat определены как структуры C #. Это потому, что программа взаимодействует между F # / C # и C ++ / CLI. В конце концов, я собираюсь загрузить код на BitBucket, но он находится в стадии разработки. У меня есть некоторые проблемы, которые нужно решить, прежде чем я смогу это сделать.
[StructLayout(LayoutKind.Sequential)]
public struct Body
{
public VectorFloat Position;
public float Size;
public uint Color;
public VectorFloat Velocity;
public VectorFloat Acceleration;
'''
}
[StructLayout(LayoutKind.Sequential)]
public partial struct VectorFloat
{
public System.Single X { get; set; }
public System.Single Y { get; set; }
public System.Single Z { get; set; }
}
Вектор определяет тип операторов, которые вы ожидаете от стандартного класса Vector. Вы, вероятно, могли бы использовать класс Vector3D из .NET Framework для этого случая (на самом деле я изучаю его переход).
ОБНОВЛЕНИЕ 2: Улучшенный код, основанный на первых двух ответах ниже:
for i = 0 to bodies.Length - 1 do
for j = (i + 1) to bodies.Length - 1 do
let r = ( bodies.[j].Position - bodies.[i].Position)
let rm = (float32)r.MagnitudeSquared + softeningLengthSquared
let f = r / (Math.Sqrt((float)rm) * (float)rm)
bodies.[i].Acceleration <- bodies.[i].Acceleration + (f * bodies.[j].Mass)
bodies.[j].Acceleration <- bodies.[j].Acceleration - (f * bodies.[i].Mass)
bodies.[i].Position <- bodies.[i].Position + bodies.[i].Velocity * dT
bodies.[i].Velocity <- bodies.[i].Velocity + bodies.[i].Acceleration * dT
Ветвь в силовой функции для покрытия случая b1 == b2 - худший нарушитель. Вам это не нужно, если sofifyingLength всегда ненулевое, даже если оно очень маленькое (Epsilon). Эта оптимизация была в коде C #, но не в версии F # (дох!).
Math.Pow (x, -1.5) кажется намного медленнее, чем 1 / (Math.Sqrt (x) * x). По сути, этот алгоритм немного странный, потому что его производительность продиктована стоимостью одного шага.
Перемещение вычисления силы в линию и избавление от некоторых разрывов также дает некоторое улучшение, но производительность действительно снижалась из-за ветвления и во власти стоимости Sqrt.
WRT с использованием классов над структурами: в некоторых случаях (CUDA и нативные реализации этого кода на C ++ и средство визуализации DX9) мне нужно получить массив тел в неуправляемый код или на графический процессор. В этих сценариях способность запоминать непрерывный блок памяти выглядит как путь. Не то, что я получил бы от массива класса Body.