[Math] Ball motion with air resistance coupled differential equation for fourth-order Runge-Kutta

MATLABordinary differential equationsrunge-kutta-methods

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
theta = 30*pi/180; % rad
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')

Best Answer

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');

Figure 1

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.

Related Question