Реализация метода интегрирования Верле

Question

Всем доброго времени суток! Я хотела реализовать метод интегрирования Верле (https://en.wikipedia.org/wiki/Verlet_integration) без использования скорости, но наткнулась на проблему наличия текущей и предыдущей координат. Я ввожу класс координат с тремя экземплярами (1D, 2D и многомерный случай):

class Coord coord where  
    infixl 6 +++
    (+++) :: coord -> coord -> coord
    infixl 7 ***
    (***) :: Double -> coord -> coord
    infixl 6 -.-
    (-.-) :: coord -> coord -> coord
    (-.-) coord1 coord2 = coord1 +++ (-1.0)***coord2

instance Coord Double where
    coord1 +++ coord2 = coord1 + coord2
    coord1 *** coord2 = coord1 * coord2

instance Coord (Double, Double) where
    (coord11,coord12) +++ (coord21,coord22) = (coord11+coord21,coord12+coord22)
    num *** (coord11,coord12) = (num*coord11,num*coord12)

instance Coord [Double] where 
    list1 +++ list2 = zipWith (+) list1 list2
    num *** list = map (* num) list

А дальше делаю саму функцию интегрирования с входными параметрами: ускорение (функция), начальная координата, начальная скорость и шаг времени, а возвращаю бесконечный список положений точки после каждого шага

verlet :: Coord coord => (coord -> coord) -> coord -> coord -> Double -> [coord]
verlet a x0 v0 dt = 
    let x1 = x0 +++ (dt***v0) +++ (0.5*(dt**2))***(a x0)
        verlet' xPrev xNow = let xNext = 2***xNow -.- xPrev +++ (dt**2))***(a xNow)
                            in (xPrev,xNow) : verlet' xNow xNext
    --- in (,) : verlet'  - вот тут проблема

Если что, я самый настоящий новичок, так что, скорее всего, это решается легко, но мне не поддается, поэтому прошу помощи)

Answer 1

verlet :: Coord coord => (coord -> coord) -> coord -> coord -> Double -> [coord] Странная система... На плюсах реализуется вот как:

struct Body
{
    Vec3d pos { 0.0, 0.0, 0.0 };
    Vec3d vel { 2.0, 0.0, 0.0 }; // 2m/s along x-axis
    Vec3d acc { 0.0, 0.0, 0.0 }; // no acceleration at first
    double mass = 1.0; // 1kg
    double drag = 0.1; // rho*C*Area - simplified drag for this example

    /**
     * Update pos and vel using "Velocity Verlet" integration
     * @param dt DeltaTime / time step [eg: 0.01]
     */
    void update(double dt)
    {
        Vec3d new_pos = pos + vel*dt + acc*(dt*dt*0.5);
        Vec3d new_acc = apply_forces(); // only needed if acceleration is not constant
        Vec3d new_vel = vel + (acc+new_acc)*(dt*0.5);
        pos = new_pos;
        vel = new_vel;
        acc = new_acc;
    }

    Vec3d apply_forces() const
    {
        Vec3d grav_acc = Vec3d{0.0, 0.0, -9.81 }; // 9.81m/s^2 down in the Z-axis
        Vec3d drag_force = 0.5 * drag * (vel * abs(vel)); // D = 0.5 * (rho * C * Area * vel^2)
        Vec3d drag_acc = drag_force / mass; // a = F/m
        return grav_acc - drag_acc;
    }
};

Может, поможет)

Answer 2

У вас уже есть x0, который вам дан, и x1, который вы посчитали в соответствии с инструкциями из Википедии. Далее в соответствии с теми же инструкциями просто используйте их в качестве первых двух аргументов для verlet':

verlet :: Coord coord => (coord -> coord) -> coord -> coord -> Double -> [coord]
verlet a x0 v0 dt = 
    let x1 = x0 +++ (dt***v0) +++ (0.5*(dt**2))***(a x0)
        verlet' xPrev xNow = let xNext = 2***xNow -.- xPrev +++ (dt**2))***(a xNow)
                             in (xPrev,xNow) : verlet' xNow xNext
    in verlet' x0 x1

Дополнение:

Судя по вашему комментарию, вам на выходе нужен список из отдельных координат, а не их пар. Если так, то просто сделайте, чтобы функция verlet' возвращала именно такой список:

...
in xPrev : verlet' xNow xNext
...

БЛОГ НА HUSL

Реализация метода интегрирования Верле

Ответы (2 шт):