- #1

Bluestribute

- 194

- 0

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.

You are using an out of date browser. It may not display this or other websites correctly.

You should upgrade or use an alternative browser.

You should upgrade or use an alternative browser.

- Thread starter Bluestribute
- Start date

- #1

Bluestribute

- 194

- 0

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.

- #2

DrClaude

Mentor

- 7,943

- 4,602

That's a numerical methods, so I have moved the thread to Programming.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).

Start by setting out the differential equations you need to solve.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).

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).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.

- #3

rcgldr

Homework Helper

- 8,791

- 580

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).set of coupled first-order ODEs.

Code:

```
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)
```

Last edited:

- #4

Bluestribute

- 194

- 0

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?

- #5

rcgldr

Homework Helper

- 8,791

- 580

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.So in the example above, are you using the K values only on calculating acceleration?

Last edited:

- #6

Bluestribute

- 194

- 0

- #7

rcgldr

Homework Helper

- 8,791

- 580

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:thinkif I can get that relationship, I could figure out the implementation . . . Ithink

Code:

```
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.

Last edited:

Share:

- Last Post

- Replies
- 8

- Views
- 502

- Replies
- 11

- Views
- 631

- Replies
- 0

- Views
- 339

- Last Post

- Replies
- 32

- Views
- 834

- Replies
- 8

- Views
- 575

C/C++
Operating system and C++

- Last Post

- Replies
- 22

- Views
- 804

- Replies
- 36

- Views
- 1K

- Replies
- 14

- Views
- 515

- Last Post

- Replies
- 8

- Views
- 724

- Replies
- 11

- Views
- 509