Rungekuttan for a system of 2 DOE functions

Hello,

I have been struggling with creating a code where the goal is to simulate a tennis serve that is desscribed by the differential equations;

m¨x = −KxV 1.5

m¨y = −mg − KyV 1.5

V = √ ˙x2 + ˙y2

I have previously only worked with one equation with RK, so this obviusly creates a problem for me. But I tried to code something in this order by using RK only for se how the velocities changes on tthe x and y-direction by using the acceleration equations;

% helpfunction for acceleration in x direction 
ax = @(Vx, Vy) (-Kx * Vtot(Vx, Vy)^1.5) / m;

% helpfunction for acceleration in y direction 
ay = @(Vx, Vy) (-g - (Ky * Vtot(Vx,Vy)^1.5) / m);

and at the end calculating the lenght in x and y dirction with the velocities*time;

    x= x + dt*Vx;

    y= y + dt*Vy;

However the results i get from the plot seems unreasonalbe. First the tennisball travels 15 meter in the height, y-direction and hits the ground only after 1.2 s. This seems pretty unreasonable. I highly suspect I´m doing something wrong in the RK code but I cannot figure out what is wrong.

Here is the full code.

clear, clc;

% Intial 
m = 0.58;
Kx = 0.01;
Ky = 0.02;
g = 9.82;

% time steps 
dt = 0.01;
t_start = 0;
t_end = 6;
ttot=0;

% Initial positions and velocities
x0 = 0;
y0 = 2.3;
Vstart= 200/3.6; %m/s
alfa = 15;
alfarad = deg2rad(alfa);
Vx = 200*cos(alfarad);
Vy = 200*sin(alfarad);

% empty matrices
positions_x = [];
hast_x = [];
positions_y = [];
hast_y=[];
speeds = [];
hast_tot=[];
times = [];

% helpfunction for the total velocity
Vtot = @(Vx, Vy) sqrt(Vx^2 + Vy^2);

% helpfunction for acceleration in x direction 
ax = @(Vx, Vy) (-Kx * Vtot(Vx, Vy)^1.5) / m;

% helpfunction for acceleration in y direction 
ay = @(Vx, Vy) (-g - (Ky * Vtot(Vx,Vy)^1.5) / m);

y=y0;
x=x0;
% Runge-Kutta code
while t_start < t_end
    positions_x = [positions_x, x];
    positions_y = [positions_y, y];
    hast_x = [hast_x, Vx ];
    hast_y = [hast_y, Vy];
    hast_tot = [hast_tot, sqrt(Vx^2+Vy^2)];
    ttot = ttot + dt;
    times = [times, ttot];

    k1_x = dt*ax(Vx, Vy);
    k1_y = dt*ay(Vx, Vy);

    k2_x = dt*ax(Vx + 0.5 * dt, Vy + k1_y * 0.5);
    k2_y = dt*ay(Vx + k1_x * 0.5, Vy + 0.5 * dt);

    k3_x = dt*ax(Vx + 0.5 * dt , Vy + k2_y * 0.5);
    k3_y = dt*ay(Vx + k2_x * 0.5, Vy + 0.5 * dt);

    k4_x = dt*ax(Vx+dt, Vy + k3_y);
    k4_y = dt*ay(Vx + k3_x, Vy + dt);
    
    Vx = Vx + (1 / 6) * (k1_x + 2 * k2_x + 2 * k3_x + k4_x);
    Vy = Vy + (1 / 6) * (k1_y + 2 * k2_y + 2 * k3_y + k4_y);
    
    x= x + dt*Vx;

    y= y + dt*Vy;

    t_start = t_start + dt;

    if y<0
        disp('vi kom till y=0');
        break
    end

end

% plot the results
subplot(5, 1, 1);
plot(times, positions_x);
xlabel('Tid');
ylabel('x-position');

subplot(5, 1, 2);
plot(times, hast_x);
xlabel('Tid');
ylabel('x-hast');

subplot(5, 1, 3);
plot(times, positions_y);
xlabel('Tid');
ylabel('y-position');

subplot(5, 1, 4);
plot(times, hast_y);
xlabel('Tid');
ylabel('y-hast');

subplot(5, 1, 5);
plot(times, hast_tot);
xlabel('Tid');
ylabel('Hastighet');

Answer 1

[William Rose]

@Emelia,

I recommend that you write a set of first order differential equations for your system.

Let the four variables be x, vx, y, vy.

dx/dt=vx

dvx/dt=-Kx*(vx^2+vy^2)^0.75/m

dy/dt=vy

dvy/dt=-Ky*(vx^2+vy^2)^0.75/m-g

If you don't want to use one of the built-in Matlab ODE solvers, such as ode45(). then the formalisim above should be useful.

If you want to use one of the built-in Matlab ODE solvers, then put the equations above into the form that ode45() likes, which is to create vector of length 4, which we will call y, and a function that returns the first derivative, which we will call dydt:

dydt(1)=y(2);

dydt(2)=-Kx*(y(2)^2+y(4)^2)^0.75/m

dydt(3)=y(4);

dydt(4)=-Ky*(y(2)^2+y(4)^2)^0.75/m-g;

That will get you ready to use ode45(). See the help for details on how to call it.

Answer 2

[James Tursa]

It would be better to have a single state vector for your solution so that you can easily compare your results to MATLAB ode45( ) function results. However, for what you have coded, note that your differential equations are not functions of time or delta time. So there should be no dt in the inputs. E.g., you have this:

k1_x = dt*ax(Vx, Vy);
k1_y = dt*ay(Vx, Vy);
k2_x = dt*ax(Vx + 0.5 * dt, Vy + k1_y * 0.5);
k2_y = dt*ay(Vx + k1_x * 0.5, Vy + 0.5 * dt);
k3_x = dt*ax(Vx + 0.5 * dt , Vy + k2_y * 0.5);
k3_y = dt*ay(Vx + k2_x * 0.5, Vy + 0.5 * dt);
k4_x = dt*ax(Vx+dt, Vy + k3_y);
k4_y = dt*ay(Vx + k3_x, Vy + dt);

Why is dt used as an input to your derivative functions? It should look something like this instead:

k1_x = dt*ax(Vx             , Vy             );
k1_y = dt*ay(Vx             , Vy             );
k2_x = dt*ax(Vx + k1_x * 0.5, Vy + k1_y * 0.5);
k2_y = dt*ay(Vx + k1_x * 0.5, Vy + k1_y * 0.5);
k3_x = dt*ax(Vx + k2_x * 0.5, Vy + k2_y * 0.5);
k3_y = dt*ay(Vx + k2_x * 0.5, Vy + k2_y * 0.5);
k4_x = dt*ax(Vx + k3_x      , Vy + k3_y      );
k4_y = dt*ay(Vx + k3_x      , Vy + k3_y      );

I've lined up the input arguments to make the code more readable.

So you use RK4 for the velocity integration. But then you use Euler 1st order method to integrate the position equations (and you use updated velocity instead of present velocity):

x= x + dt*Vx;
y= y + dt*Vy;

It is best to integrate all states with the same method. In your case, you should be simultaneously integrating the position and velocity states all using RK4 method instead of the mixed methods you are currently using.