# Object falling with air resistance using Runge-Kutta

I am not very familiar with differential equations, nor physics in general. I am trying to program an object falling with air resistance with the use of a numerical algorithm called Runge-Kutta. The object is going to fall from the same position every-time $x=0$, with the same start speed $v=0$ every time.

I think I then need to work with the equation: $ma=F_d-mg$, where $F_d$ is the drag and it's variables are $\frac {1}{2} pdAv^2$.

$ma=F_d-mg \rightarrow ma = \frac {1}{2} pdAv^2-mg$, where $A$ is the area, a constant.

If I solve for $a$, the acceleration I get the equation:

$a = \frac {\frac 12pdav^2}{m}-g$

This is the differential equation:

$y'' = \frac {\frac 12pda(y')^2}{m}-g$

Because I am going to solve for this using the Runge-Kutta numerical method I am going to have to rewrite this into two first order differentials equations. ¨

I thus start by saying:

$\frac{dx}{dt} = v$

and then

$\frac{dv}{dt}=\frac {\frac 12pda(v)^2}{m}-g$

Now that I have my two first order differential equations I should be fine plugging them into the Runge-Kutta method, to obtain the velocity and the position.

Here is a picture of what I have programmed. Using the initial values:

• $p = 0.5$
• $d = 1.275$
• $a = 0.003184$
• $m = 0.14529$
• $g = 9.81$

using UnityEngine;
using System.Collections;

public class ODESOLVER1 : MonoBehaviour {}

float t_n = 0;
float y_n = 0;
float h = 0.1;
float x_n = 0;
float p = 0.5;
float d = 1.275;
float a = 0.003184;
float m = 0.14529;
float g = 9.81;

void Start () {
for (int i = 0; i < 200; i++)
{
RungaKrutta4(t_n, y_n, x_n, h);
Debug.Log("t_n=" + t_n + "y_n=" + y_n + "x_n=" + x_n);
};
}

void Update () {}

private void repeatKrutta() {}

private void RungaKrutta4(float t_n, float y_n, float x_n, float h) {

float k0 = h * function1(t_n, y_n, x_n);
float l0 = h*function2(t_n, y_n, x_n);
float k1 = h * function1(t_n + (h / 2), y_n + (k0 / 2), x_n + (l0 / 2));
float l1 = h * function2(t_n + (h / 2), y_n + (k0 / 2), x_n + (l0 / 2));
float k2 = h * function1(t_n + (h / 2), y_n + (k1 / 2), x_n + (l1 / 2));
float l2 = h*function2(t_n + (h/2), y_n + (k1/2), x_n + (l1/2));
float k3 = h * function1(t_n + (h / 2), y_n + (k2 / 2), x_n + (l2 / 2));
float l3 = h*function2(t_n + (h/2), y_n + (k2/2), x_n + (l2/2));

this.t_n = t_n + h;
this.y_n = y_n + ((h / 6) * (k0 + (2 * k1) + (2 * k2) + k3));
this.x_n = x_n + ((h / 6) * (l0 + (2 * l1) + (2 * l2) + l3));
}

private float function1(float t, float y, float x) {
return (((0.5f * p * d * a * y * y)/m) - g);
}

private float function2(float t, float y, float x) {
return (y);
}

private void calculateDecimals(float value) {
string s = System.Convert.ToString(value);
}


My $y_{n+1}$ almost comes out correctly except for the fact that it becomes a much smaller value. If the speed should have been $0.0981$ with a timestep $h = 0.010$, it comes out as $0.000981.$ in the first iteration for example.

My $x_{n+1}$ does not come correctly out at all. The first, second and third values should be as follows: $-0.00049, -0.00196, -0.0441$ with a timestep of $h=0.010.$

• I'm voting to close this question as off-topic because it's either a check-my-work question or a question about debugging code, which are both off-topic. – ACuriousMind Sep 23 '16 at 12:05
• But isn't there a possibility that you could at least answer the physics part of the problem? – David Lund Sep 23 '16 at 12:09
• as a hint: check how the deviation between what you get and what you should get changes with $h$ (or more bluntly: how many powers of h should you have in a first order integration) – Bort Sep 23 '16 at 12:18
• 1) There is only one "r" in Runge-Kutta. 2) I suspect that h (which you have not shown) has a value of .01, which explains the small velocity change. 3) This in turn will produce very small velocity increments, but you have not indicated what is wrong with this, other than "does not come correctly out at all", so no one can help you there. – WhatRoughBeast Sep 23 '16 at 12:23
• What is the objective of this exercise? Are you trying to learn some basic physics, Java programming, this UnityEngine software package, or how the Runge-Kutta algorithm works? Your implementation of Runge-Kutta may be correct but looks pretty convoluted to me. Most scientists/engineers don't try to implement their own ODE solver but instead use one of the many freely-available solvers. – Bill Greene Sep 23 '16 at 14:20

I am not going to compile everything to see if this is your only mistake (you can check here for an RK4 implementation). However, one clear error is that the fourth step has a stepsize h and not h/2. So you should have

float k3 = h * function1(t_n + h , y_n + k2, x_n + l2);
float l3 = h*function2(t_n + h, y_n + k2, x_n + l2);


Just for reference, the method with all halves (like you had) is only $\mathcal{O}(h^2)$, but it's more stable than RK4.

(2 * k1) + (2 * k2)


is unnecessary,

2 * (k1 * k2)


gets rid of an operation. Also, the update step should be

this.y_n = y_n + ((1 / 6) * (k0 + (2 * k1) + (2 * k2) + k3));
this.x_n = x_n + ((1 / 6) * (l0 + (2 * l1) + (2 * l2) + l3));


i.e. get rid of the $h$'s there. Most people don't know this, but the RK methods actually respect units / dimensional analysis. Check the dimensions (where $f$ is a rate): function1(t_n, y_n, x_n) is a rate, so h * function1(t_n, y_n, x_n) has the same dimensions as y_n. Since you multiplied the h already, your k's aren't rates, but actually the same units as the y_n and x_n and so the last step is just a linear combination.

There's no reason to t_n + (h / 2) so many times. Save that value to cut down on operations. RungaKrutta4 should be RungaKutta4

These are all of the mistakes I could find while quickly skimming over your RK4, and RK4 is one of the simplest methods (you have no adaptive timestepping, stabilization, interpolation, or optimizations like within-method multithreading, explicit fused-multiply added, etc.). This is one reason why, for most people, I would recommend that they use tried and tested ODE solver packages.

• So I removed the h, fixed k3 and l3, and I figured out that 1/6 also caused me some troubled so I changed it into a floating point number instead and now it works as it should. Thanks. – David Lund Sep 23 '16 at 17:22