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')
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. $\endgroup$