2
$\begingroup$

I am attempting to simulate projectile flight with drag. However, with a timestep of 0.1 seconds, I am consistently getting an error of ~0.1-1%.

while p[1, iterr] >= 0 and iterr < 9999:
            p0 = p[:, iterr]
            v0 = v[:, iterr]
            pp1 = h * v0
            vp1 = h * f(v0)
            pp2 = h * (v0 + vp1)
            vp2 = h * f(v0 + vp1/2)
            pp3 = h * (v0 + vp2)
            vp3 = h * f(v0 + vp2/2)
            pp4 = h * (v0 + vp3)
            vp4 = h * f(v0 + vp3)

            v[:, iterr + 1] = v0 + (vp1 + 2 * (vp2 + vp3) + vp4)/6
            #p[:, iterr + 1] = p0 + h * (pp1 + 2 * (pp2 + pp3) + pp4)/6
            p[:, iterr + 1] = p0 + h * v[:, iterr + 1]
            iterr += 1

Here,

def f(v):
    magv = np.linalg.norm(v)
    return np.array([k * magv * v[0], k * magv * v[1] - 9.81])

Note: One line updating position is commented out because I was/am trying to figure out if the two expressions make any significant difference, and they do not.

Is there something about this simulation that I am doing wrong, or is 0.1-1% characteristic of the RK4 method?

$\endgroup$
  • $\begingroup$ What happens if you cut the time step in half? If your code is working correctly, you should have $1/16$ of the error. Do you have some sort of closed form for the solution? $\endgroup$ – fred Oct 21 '15 at 0:19
  • $\begingroup$ To measure error, I am using the code for my dragged-motion simulation with k = 0. If you notice that sets acceleration to [0, -9.81], which is ideal projectile motion acceleration. What's interesting is that changing the timestep from 0.01 to 0.09 changes error from 2% to 0.8%, for a particular velocity (the error seems to be a function of velocity as well, which is strange.) $\endgroup$ – Chronum Oct 21 '15 at 1:01
  • $\begingroup$ It is quite unusual to have a time step 9 times as large and the error reduce by more than half. To answer your question, we usually check the scaling of the error to make sure the algorithm is working as it should. So, no, 0.1-1% is not a characteristic of RK4. $\endgroup$ – fred Oct 21 '15 at 2:05
  • $\begingroup$ @fred -- RK4 is fourth order, and $(0.9)^4=0.66$ So not quite a 60% reduction, but if you're in the pre-asymptotic range, a reduction by 60% doesn't seem completely out of the picture. In any case though, reduce the time step by factors of 2 and see how the error becomes smaller. $\endgroup$ – Wolfgang Bangerth Oct 21 '15 at 2:20
  • $\begingroup$ 1.43320277 0.75028814. Percentage error for 0.05 seconds and 0.1 seconds respectively. I have a feeling there's something wrong in my code in terms of incrementing the position variable, but I can't quite figure out what. $\endgroup$ – Chronum Oct 21 '15 at 2:37
3
$\begingroup$

This is a popular question but I can't find a readily available answer. So here are some of the details.

Let us assume that you are solving the equation $$ m \dot{\mathbf{v}} = m \mathbf{g} - k \lVert \mathbf{v}\rVert \mathbf{v} $$ where $m$ is the mass of the projectile, $\mathbf{v}$ is its velocity, $\mathbf{g}$ is the acceleration due to gravity, $k$ is a drag coefficient , $\dot{\mathbf{v}}$ is the time-derivative of the velocity, and $\lVert \mathbf{v} \rVert$ is the magnitude of the velocity. Also assume that the initial velocity is $ \mathbf{v}(0) = \mathbf{v}_0 \,. $

If we define $$ \mathbf{f}(t, \mathbf{v}) := \mathbf{g} - \frac{k}{m}\,\lVert\mathbf{v}\rVert\,\mathbf{v} \,. $$ we can write our equation as $$ \dot{\mathbf{v}} = \mathbf{f}(t, \mathbf{v}) $$ The RK4 discretization with a stepsize $h = t_{n+1} - t_n$ can then be expressed as $$ \mathbf{v}_{n+1} = \mathbf{v} + \frac{h}{6}(\mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3 + \mathbf{k}_4) $$ where $$ \mathbf{k}_1 = f(t_n, \mathbf{v}_n) \,, \mathbf{k}_2 = f(t_n + \frac{h}{2}, \mathbf{v}_n + \frac{h}{2} \mathbf{k}_1) \,, \mathbf{k}_3 = f(t_n + \frac{h}{2}, \mathbf{v}_n + \frac{h}{2} \mathbf{k}_2) \,, $$ and $$ \mathbf{k}_4 = f(t_n + h, \mathbf{v}_n + h \mathbf{k}_3) \,. $$ A similar approach is used to find the position $\mathbf{x}$ by integrating $$ \dot{\mathbf{x}} = \mathbf{v} $$ with initial position $\mathbf{x}(0) = \mathbf{0}$.

For our 2D Python implementation we can define the function $f$ as

def f(t, v, g, m, k):
    kmv = np.linalg.norm(v)*(k/m)
    return (g - kmv*v)

The RK4 method can be implemented as

def RK4(tn, xn, vn, h, g, m, k):

    k1 = f(tn, vn, g, m, k)
    k2 = f(tn + h/2, vn + k1*h/2, g, m, k)
    k3 = f(tn + h/2, vn + k2*h/2, g, m, k)
    k4 = f(tn + h, vn + k3*h, g, m, k)
    vn1 = vn + (k1 + 2*k2 + 2*k3 + k4)*(h/6)

    k1x = vn
    k2x = vn + k1x*h/2
    k3x = vn + k2x*h/2
    k4x = vn + k3x*h
    xn1 = xn + (k1x + 2*k2x + 2*k3x + k4x)*(h/6)

    return vn1, xn1

Finally, the iterative process can be implemented using

 theta = math.radians(30)
 v0 = np.array([math.cos(theta), math.sin(theta)])
 x0 = np.array([0, 0])
 t0 = 0
 x = []
 v = []
 t = []
 x.append(x0)
 v.append(v0)
 t.append(t0)

 m = 1
 k = 0.01
 g = np.array([0, -9.81])
 h = 0.01

 tn = t0
 vn = v0
 xn = x0
 maxt = 1

 while (t < maxt):
   vn, xn = RK4(tn, xn, vn, h, g, m, k)
   tn = tn + h
   t.append(tn)
   x.append(xn)
   v.append(vn)

 print t[-1]
 print x[-1]
 print v[-1]

For $h = 0.01$ we get the final position $0.85746485, -4.31393019$. The stopping condition that I've used can lead to different end times depending on the value of $h$.

I have a feeling that in the original question stopping conditions are being achieved at different times.

$\endgroup$
2
$\begingroup$

You can probably just use your original code by correcting the step computation for the positions to fall within the RK4 scheme:

while p[1, iterr] >= 0 and iterr < 9999:
        p0 = p[:, iterr]
        v0 = v[:, iterr]
        pp1 = h * v0
        vp1 = h * f(v0)
        pp2 = h * (v0 + vp1/2)
        vp2 = h * f(v0 + vp1/2)
        pp3 = h * (v0 + vp2/2)
        vp3 = h * f(v0 + vp2/2)
        pp4 = h * (v0 + vp3)
        vp4 = h * f(v0 + vp3)

        v[:, iterr + 1] = v0 + (vp1 + 2 * (vp2 + vp3) + vp4)/6
        p[:, iterr + 1] = p0 + (pp1 + 2 * (pp2 + pp3) + pp4)/6
        iterr += 1

This is then a true RK4 implementation in contrast to the other 2 answers where the position update is by the Euler method or something close to it.

You should (or are already doing without documenting it here) use interpolation to get a better value for the time of the zero crossing.

$\endgroup$
  • $\begingroup$ That is indeed what I did at first. My question here is: Does it matter? Since position is a linear function of velocity (which is changing, but we've got that handled with the RK4 for velocity), do we need RK4 for both, or does it suffice to have RK4 for simply velocity, so I can make 4 approximations for the acceleration term. And yes, I am using linear interpolation for the last two points. $\endgroup$ – Chronum Oct 25 '15 at 16:02
  • $\begingroup$ Yes, it does matter, since it effects the error order of the position. And no, there are some slight differences in this code compared to the original, namely using v1=v0+vp1/2 and v2=v0+vp2/2 both in the position and velocity increment computation. You were missing the division by 2 in the position increment. $\endgroup$ – Dr. Lutz Lehmann Oct 25 '15 at 16:09
  • $\begingroup$ I just noticed (again) that Biswajit's code implements the dual-RK4 implementation for both position and velocity. And I tried it myself right now (was trying Euler for position), and the errors are indeed smaller by a sizeable factor. $\endgroup$ – Chronum Oct 25 '15 at 16:48
  • $\begingroup$ Yes, on second view, also Biswajit's code is correct, since the velocity does not depend on the position. $\endgroup$ – Dr. Lutz Lehmann Oct 25 '15 at 19:32
  • $\begingroup$ quora.com/… this is a beautiful example of SHM using RK4 with coupled approximation terms. $\endgroup$ – Chronum Oct 25 '15 at 20:15
1
$\begingroup$

Thanks in large part to Biswajit and this MIT OCW slide, page 4, and some error on my part. My initial position was set to [0 5.5] instead of [0 0] from one of my previous attempts.

Here is my final code, using both Biswajit's answer, and the MIT slides:

while p[1, iterr] >= 0 and iterr < 9999: p0 = p[:, iterr] v0 = v[:, iterr] t0 = t[0, iterr] vp1 = f(t0, v0) vp2 = f(t0 + h/2, v0 + h * vp1/2) vp3 = f(t0 + h/2, v0 + h * vp2/2) vp4 = f(t0 + h, v0 + h * vp3) v[:, iterr + 1] = v0 + h * (vp1 + 2 * (vp2 + vp3) + vp4)/6 p[:, iterr + 1] = p0 + h * v[:, iterr + 1] t[0, iterr + 1] = t0 + h iterr += 1

And,

def f(t, v): magv = np.linalg.norm(v) return np.array([k * magv * v[0], k * magv * v1 - 9.81])

Here, k is the entire drag coefficient-mass-reference-area term, and $v^2$ is the velocity dependence expressed as $$||\mathbf v||*\mathbf v$$.

If the code needs to be any clearer, I will comment out the parts. Anyone seeing this, notice that acceleration here does not depend on the time of the simulation, however, that is not always true.

Thank you to everyone who helped me look at this problem.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.