Solving Spring System with RK4: Step-by-Step Guide

[Bluestribute]

I'm not sure if this is the appropriate place, but I've just been searching, with no luck, on how to solve a spring system with RK4 (hence the math section, since it's a Runge Kutta question). How would you set it up to change position based on spring properties? I haven't really found any consistent things either, except that x' = v and v' = a (which is standard).

What would one iteration look like? As far as having position and velocity as your inputs, and outputting an updated position after the iteration. I know how to do RK4, just not so great at implementing it.

 

[DrClaude]

I'm not sure if this is the appropriate place, but I've just been searching, with no luck, on how to solve a spring system with RK4 (hence the math section, since it's a Runge Kutta question).

That's a numerical methods, so I have moved the thread to Programming.

How would you set it up to change position based on spring properties? I haven't really found any consistent things either, except that x' = v and v' = a (which is standard).

Start by setting out the differential equations you need to solve.

What would one iteration look like? As far as having position and velocity as your inputs, and outputting an updated position after the iteration. I know how to do RK4, just not so great at implementing it.

Once you have done the above, applying the RK4 algorithm is trivial (meaning that there is nothing special, it will look the same as it does for any set of coupled first-order ODEs).

 

[rcgldr]

set of coupled first-order ODEs.

That's probably the issue, knowing how to use RK4 with coupled ODEs. Example sequence, ak1 -> ak4 are calculated accelerations, vk1 -> vk4 are calculated velocities, p[ i ] is current position. Note that vk2 is based off ak1, vk3 is based off ak2, ... ; trying to use vk2 based off ak2, vk3 based off ak3, ... , will throw off the RK4 coupled ODE calculations. Similarly, ak2 is based off vk1, ak3 is based off vk2 (in case the code was changed to calculate vk's before ak's).

    ak1 = a(p[i])
    vk1 = v[i]
    ak2 = a(p[i] + vk1 (Δt/2))
    vk2 = v[i] + ak1 (Δt/2)
    ak3 = a(p[i] + vk2 (Δt/2))
    vk3 = v[i] + ak2 (Δt/2)
    ak4 = a(p[i] + vk3 (Δt))
    vk4 = v[i] + ak3 (Δt))

    v[i+1] = v[i] + (Δt/6) (ak1 + 2 ak2 + 2 ak3 + ak4)
    p[i+1] = p[i] + (Δt/6) (vk1 + 2 vk2 + 2 vk3 + vk4)

 

[Bluestribute]

So in the example above, are you using the K values only on calculating acceleration? Then do you solve it like normal, or would you repeat that process for X and V, where V would be solved using K values? Cause the way I know it is that the K values are in the form h*f (where f is the multivar function and you have different additives for each K1-4), and you're only doing that with accelerations . . .

Or would I just setup a spring system as:

v' = v - x (or whatever it is, though I didn't see much of what it should be . . . ) and do it similarly? Where I would solve first vk1 by plugging in v,x, and then solve xk1 by . . . uh . . . is it the same as above, the "k1"*Δt/2?

 

[rcgldr]

So in the example above, are you using the K values only on calculating acceleration?

As mentioned in my example, both accelerations and velocities are calculated as K values: ak1->ak4 are accelerations, vk1->vk4 are velocities, calculated during a time step Δt. Then the accelerations ak1->ak4 are used to update velocity (v[ i+1 ]) and velocities vk1->vk4 are used to update position (p[ i+1 ]). I'm assuming that acceleration is a function of position a(p[ i ]). In some cases, like re-entry of a space craft, acceleration would be a function of velocity (aerodynamic drag) and position (density of atmosphere, force of gravity). For an example spring like case, if a(p[ i ]) = -p[ i ], you get some variation of a sine wave, depending on initial condition, like v[0] = 1, p[0] = 0.

 

[Bluestribute]

So when you say acceleration is a function of position, are you saying F = ma = -kΔx (english language to math language is not my strong suit! Haha)? So would ak1 just be -kΔx/m? I think if I can get that relationship, I could figure out the implementation . . . I think

 

[rcgldr]

So when you say acceleration is a function of position, are you saying F = ma = -kΔx (english language to math language is not my strong suit! Haha)? So would ak1 just be -kΔx/m? I think if I can get that relationship, I could figure out the implementation . . . I think

Almost, for your example a(x) = -kx / m, based on position x, not change in position Δx. So ak1 = - k x[ i ] / m. x[0] would be the initial position, and v[0] would be the initial velocity. Note I corrected the formulas for ak4 and vk4 (they use Δt not Δt/2) Rewriting the example loop fragment using x instead of p:

    ak1 = a(x[i])
    vk1 = v[i]
    ak2 = a(x[i] + vk1 (Δt/2))
    vk2 = v[i] + ak1 (Δt/2)
    ak3 = a(x[i] + vk2 (Δt/2))
    vk3 = v[i] + ak2 (Δt/2)
    ak4 = a(x[i] + vk3 (Δt))
    vk4 = v[i] + ak3 (Δt))

    v[i+1] = v[i] + (Δt/6) (ak1 + 2 ak2 + 2 ak3 + ak4)
    x[i+1] = x[i] + (Δt/6) (vk1 + 2 vk2 + 2 vk3 + vk4)

ak1 and vk1 are the initial acceleration and velocity at the start of the current time step. x[ i ] + vk1 (Δt/2) is the first estimated position at the mid-point of the current time step and so a(x[ i ] + vk1 (Δt/2)) would be the first estimated acceleration at the mid-point of the current time step. ak2 and vk2 are then estimated acceleration and velocity at the mid-point of the current time step. ak3 and vk3 are improved estimates of acceleration and velocity at the mid-point of the current time step, using vk2 and ak2 as feedback. ak4 and vk4 are estimates of acceleration and velocity at the end of the current time step, using vk3 and ak3 as feedback.