Object falling with air resistance using Runge-Kutta

Question

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

Answer 1

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.