# Ball motion with air resistance coupled differential equation for fourth-order Runge-Kutta

I've created a MATLAB function for solving coupled differential equation with the fourth-order Runge-Kutta method based on what is provided here (Simultaneous Equations of First Order). Here the function:

function [y] = rk4_c(f, g, h, x, y, z, n)
% Runge-Kutta
% Implementation of the fourth-order method for coupled equations
% h = dt
% x is the time here

for ii=1:(n-1)
k1 = h * f(x(ii), y(ii), z(ii));
l1 = h * g(x(ii), y(ii), z(ii));
k2 = h * f(x(ii) + 0.5*h, y(ii) + 0.5*k1, z(ii) + 0.5*l1);
l2 = h * g(x(ii) + 0.5*h, y(ii) + 0.5*k1, z(ii) + 0.5*l1);
k3 = h * f(x(ii) + 0.5*h, y(ii) + 0.5*k2, z(ii) + 0.5*l2);
l3 = h * g(x(ii) + 0.5*h, y(ii) + 0.5*k2, z(ii) + 0.5*l2);
k4 = h * f(x(ii) + h, y(ii) + k3, z(ii) + l3);
l4 = h * g(x(ii) + h, y(ii) + k3, z(ii) + l3);
y(ii+1) = y(ii) + (1/6)*(k1 + 2*k2 + 2*k3 + k4);
z(ii+1) = z(ii) + (1/6)*(l1 + 2*l2 + 2*l3 + l4);
end


I would like to use this function for solving motion equations for a ball with air resistance. So starting with the basic equations (I hope they are correct): $$md^2x/dt^2 = F_{D_x}$$ and $$md^2y/dt^2 = F_{D_y}-mg$$ with $$F_{D_x} = -Dvvcos(\theta)$$ and $$F_{D_y} = -Dvvsin(\theta)$$ where $D$ is a drag coefficient. From my understanding, and from what I read, both second-order equation can be converted to two coupled equations that could be used with my coupled Runge-Kutta method function. Can someone help me understand how to convert my equation to a set of coupled equations that could be used with my MATLAB function?

EDIT

I was able to use my rk4_c function, however it only works for solving the velocities. I'm now trying to see how I could solve for the positions x and y. I'm kind of learning about Runge-Kutta methods at the same time, so sorry for the questions that may sound trivials.

% Runge Kutta code to solve projectile motion with quadratic drag
% dVx/dt = -(D/m)*vx*sqrt(vx^2+vy^2)
% dVy/dt = -(D/m)*vy*sqrt(vx^2+vy^2) - g

clc
clear all

% Constant
D = 0.24; %
m = 2; % kg
g = 9.80665; % m/s^2

% Define function handles
fVx = @(t,vx,vy) -(D/m)*vx*sqrt(vx^2+vy^2);
fVy = @(t,vx,vy) -(D/m)*vy*sqrt(vx^2+vy^2) - g;

% Initial conditions
v0 = 200; % m/s
t(1) = 0;
vx(1) = v0*cos(theta);
vy(1) = v0*sin(theta);

% Step size
h = 0.01; % s
tFinal = 2;
N = ceil(tFinal/h);

% RK4 simultaneous coupled loop
for ii = 1:N
% Update time
t(ii+1) = t(ii) + h;

% Update vx and vy
k1 = h * fVx(t(ii), vx(ii), vy(ii));
l1 = h * fVy(t(ii), vx(ii), vy(ii));
k2 = h * fVx(t(ii) + 0.5*h, vx(ii) + 0.5*k1, vy(ii) + 0.5*l1);
l2 = h * fVy(t(ii) + 0.5*h, vx(ii) + 0.5*k1, vy(ii) + 0.5*l1);
k3 = h * fVx(t(ii) + 0.5*h, vx(ii) + 0.5*k2, vy(ii) + 0.5*l2);
l3 = h * fVy(t(ii) + 0.5*h, vx(ii) + 0.5*k2, vy(ii) + 0.5*l2);
k4 = h * fVx(t(ii) + h, vx(ii) + k3, vy(ii) + l3);
l4 = h * fVy(t(ii) + h, vx(ii) + k3, vy(ii) + l3);
vx(ii+1) = vx(ii) + (1/6)*(k1 + 2*k2 + 2*k3 + k4);
vy(ii+1) = vy(ii) + (1/6)*(l1 + 2*l2 + 2*l3 + l4);
end

% Plot the solution
figure(1)
plot(t, vx, t,vy)
xlabel('Time (s)')
ylabel('Velocities (m/s)')
legend('Vx', 'Vy')

• A second order equation $y''=f(t,y,y')$ can be converted to the system $y_1'=y_2,y_2'=f(t,y_1,y_2)$. (By the way, what you want to do can be done with builtins, and it will use an adaptive step size as well. So unless there is some pedagogical reason to do this, you really don't need to do it.) – Ian May 26 '16 at 17:01
• @Ian, I know I could use the ode45 function in MATLAB, however I'm just using MATLAB as a platform for quick development. The code will be written in another programming language after development. – m_power May 26 '16 at 17:23
• OK. #1: you should still use something like vectors even in your final code. There is not much reason for this to be dimension-dependent. #2: you might want to look into how ode45 actually works so that you can imitate its adaptive step size method. – Ian May 26 '16 at 17:29
• @Ian, I don't want to use an adaptive step size method. My function does use vectors. I don't understand your #1 point. – m_power May 26 '16 at 17:41
• I'm not sure why you have separate $f,g$ in the first place. It seems to me that the interface of your code should look like the interface of ode45, plus an argument for the number of steps. (I could see a justification for having arguments for the arrays that store the result, for efficiency, but IMO this is premature optimization.) – Ian May 26 '16 at 17:42

I was trying to make your code work in the Matlab idiom.

% rk4.m

function [x,y] = rk4_c(f, tspan, y0, n)
% Runge-Kutta
% Implementation of the fourth-order method for coupled equations
% x is the time here
% More or less follows simplified interface for ode45; needs #points = n
% Thanks to @David for helpful suggestions.

h = (tspan(2)-tspan(1))/(n-1)
y = zeros(n,length(y0));
y(1,:) = y0;
x = linspace(tspan(1),tspan(2),n)';
for ii=1:(n-1)
k1 = h * f(x(ii), y(ii,:)');
k2 = h * f(x(ii) + 0.5*h, y(ii,:)' + 0.5*k1);
k3 = h * f(x(ii) + 0.5*h, y(ii,:)' + 0.5*k2);
k4 = h * f(x(ii) + h, y(ii,:)' + k3);
y(ii+1,:) = y(ii,:) + (1/6)*(k1 + 2*k2 + 2*k3 + k4)';
end

% drag.m

function yprime = drag(t,y);

global D g m
yprime = zeros(length(y),1);
v = hypot(y(2),y(4));
yprime(1) = y(2); % dx/dt
yprime(2) = -(D/m)*v*y(2); %F_x/m
yprime(3) = y(4); % dy/dt
yprime(4) = -(D/m)*v*y(4)-g; %F_y/m

% projectile.m

clear all;
close all;
global D g m

rho = 1.1; % kg/m^3
A = 0.001; % m^2
C_D = 3.2; % Or whatever...
g = 9.81; % m/s^2
D = rho*A*C_D;
m = 2.0; % kg
v0 = 686; % m/s
theta0 = 30; % °
tspan = [0 20]; % s
y0 = [0 v0*cosd(theta0) 0 v0*sind(theta0)];
[t,y] = ode45(@drag,tspan,y0);
ind = find(y(:,3)<=0);
if length(ind)>1,
ind1 = ind(2);
t(ind1)=t(ind1-1)-y(ind1-1,3)*(t(ind1)-t(ind1-1))/ ...
(y(ind1,3)-y(ind1-1,3));
y(ind1,:)=y(ind1-1,:)-y(ind1-1,3)*(y(ind1,:)-y(ind1-1,:))/ ...
(y(ind1,3)-y(ind1-1,3));
t = t(1:ind1);
y = y(1:ind1,:);
end
N = 100;
[t4,y4] = rk4_c(@drag,tspan,y0,N);
ind = find(y4(:,3)<=0);
if length(ind)>1,
ind1 = ind(2);
t4(ind1)=t4(ind1-1)-y4(ind1-1,3)*(t4(ind1)-t4(ind1-1))/ ...
(y4(ind1,3)-y4(ind1-1,3));
y4(ind1,:)=y4(ind1-1,:)-y4(ind1-1,3)*(y4(ind1,:)-y4(ind1-1,:))/ ...
(y4(ind1,3)-y4(ind1-1,3));
t4 = t4(1:ind1);
y4 = y4(1:ind1,:);
end
plot(y(:,1),y(:,3),'b-',y4(:,1),y4(:,3),'r.');
xlabel('Range (m)');
ylabel('Altitude (m)');
title('Projectile Trajectory'); Looks like it works OK.

EDIT: I experimented with the 5/11 version of the code. I added the extras functions

fx = @(t,x,y,vx,vy) vx;
fy = @(t,x,y,vx,vy) vy;


And then tracked the coordinates

m2 = h * fx(t(ii) + 0.5*h, x(ii) + 0.5*m1,y(ii) + ...
0.5*n1,vx(ii) + 0.5*k1,vy(ii) + 0.5*l1);
n2 = h * fy(t(ii) + 0.5*h, x(ii) + 0.5*m1, y(ii) + ...
0.5*n1, vx(ii) + 0.5*k1, vy(ii) + 0.5*l1);
k2 = h * fVx(t(ii) + 0.5*h, x(ii) + 0.5*m1, y(ii) + ...
0.5*n1, vx(ii) + 0.5*k1, vy(ii) + 0.5*l1);
l2 = h * fVy(t(ii) + 0.5*h, x(ii) + 0.5*m1, y(ii) + ...
0.5*n1, vx(ii) + 0.5*k1, vy(ii) + 0.5*l1);


It all seemed to work out OK.

• 1) you can replace the first for loop with y(1,:)=y0; 2) ' is conjugate transpose, .' is normal matrix transpose 3) feval is not necessary, f(t,y) is sufficient here. – David May 27 '16 at 0:55
• Whoa, that's nice to know. I thought that a column vector y0 would cause problems, but as you say it doesn't. For Hermitian adjoint vs. transpose, I didn't know that either, but in this problem we are toast if the difference matters. Not requiring feval must be relatively new. It seems to me that in the olden times it was required. That change also passed testing. The modified rk4_c.m will be edited in. – user5713492 May 27 '16 at 1:07
• No, the transpose thing shouldn't matter, but it's good to know, it can cause some very strange bugs! – David May 27 '16 at 1:29
• Thanks for your solution. Would it be possible to make it works with the algorithm that I wrote (RK4 coupled)? I think the yprime function would have to be modified (2 equations instead of 4 for the system?)? – m_power May 27 '16 at 13:08
• Yes, it could work with $2$ separate derivative functions and a user-written rk4_c function, but it's nicer with one function because one numerical integration subroutine can handle any dimension of input, and you can't compare as easily with Matlab's built-in functions as was done in this example. If you need to do it this way, try it and post all your code so potential helpers can analyze it and figure out your mistakes. – user5713492 May 27 '16 at 13:47