Unexpected Imaginary Numbers in Runge Kutta Method Code

I have the Runge Kutta method to solve a system of 3 ODEs. While testing this code, the y vector (and sometimes v) yields either imaginary numbers, infinity, or NaN. I am not sure where my error is, as I am fairly certain the RK equations are correct. Any help on the matter would be greatyl appreciated!

% Defining functions %
c = 2.12;
b = 2.28;
a = 0.013;
sigma = 0.07;
k = 6.67e-4;
mu = 1.89e-5;
d0 = 0.01;

m = @(j) 1000*(4/3)*(sqrt(j)/2)^3;    % j = d^2
R = @(t,y,v,j,r) r*v*sqrt(j)/mu;      %r is a temporary variable to later be replaced with the rho(y) function
Ct = @(t,y,v,j,r) 1+ 0.16*R(t,y,v,j,r)^(2/3);
W = @(t,y,v,j,r) r*(v^2)*sqrt(j)/sigma;
Cd =  @(t,y,v,j,r) 1 + a*((W(t,y,v,j,r) + b)^c) - a*b^c;
Fd =  @(t,y,v,j,r) 3*pi*sqrt(j)*mu*v*Ct(t,y,v,j,r)*Cd(t,y,v,j,r);

F1 = @(t,y,v,j) v ;                         %ODE 1: dy/dt = v 
F2 = @(t,y,v,j,r) Fd(t,y,v,j,r)/m(j) - 9.8; %ODE 2: dv/dt = Fd/m - g
F3 = @(t,y,v,j) -k*(d0)^2;                  %ODE 3: dj/dt = -kdo^2

h = 1;                          %Arbitrary step size
t = 0:h:100;                    %Analyzing up to t = some value (Droplet should hit the ground at t = 1500 apx.
v = zeros(1,length(t));
y = zeros(1,length(t));
j = zeros(1,length(t));
rho = zeros(1,length(t));
v(1) = 0;                       %Initial values for the IVP
y(1) = 2000;
j(1) = d0^2;
distance = zeros(1,229);

for i = 1:1:length(t) -1

    if (y(i) < 0)
        y(i) = 0;
        break;
    else 
        for k = 1:1:46      % Defining Rho(y) function via interpolation about the point y(i) with given dataset
        distance(k) = abs(y(i) - Atmos1.Hm(k));
        end                 %Note: This part of the code works without a problem, 
        for k = 1:1:46      %the error occurs when y(i) is NaN or INF or simply imaginary, which shouldn't happen
            if (distance(k) == min(distance))
            break
            end
        end
        
    x0 = Atmos1.Hm(k);      %Setting up the linear interpolation equation to define rho(y) function
    p = abs(y(i) - x0)/(Atmos1.Hm(2) - Atmos1.Hm(1)); 
    q =  p*(p-1)/2;
    rho(i) = Atmos1.rhosi(k) + p*(Atmos1.rhosi(k + 1) - Atmos1.rhosi(k)) - q*(Atmos1.rhosi(k + 2) -2*Atmos1.rhosi(k + 1)+ Atmos1.rhosi(k));
    end
   
  % Runge Kutta Method to solve system of ODEs % 
k1 =  h*F1(t(i),y(i),v(i),j(i))
p1 =  h*F2(t(i),y(i),v(i),j(i),rho(i))
s1 =  h*F3(t(i),y(i),v(i),j(i))

k2 =  h*F1(t(i)+h/2, y(i)+(1/2)*k1, v(i)+(1/2)*p1, j(i)+(1/2)*s1)
p2 =  h*F2(t(i)+h/2, y(i)+(1/2)*k1, v(i)+(1/2)*p1, j(i)+(1/2)*s1, rho(i))
s2 =  h*F3(t(i)+h/2, y(i)+(1/2)*k1, v(i)+(1/2)*p1, j(i)+(1/2)*s1)

k3 =  h*F1(t(i)+h/2, y(i)+(1/2)*k2, v(i)+(1/2)*p2, j(i)+(1/2)*s2)
p3 =  h*F2(t(i)+h/2, y(i)+(1/2)*k2, v(i)+(1/2)*p2, j(i)+(1/2)*s2, rho(i))
s3 =  h*F3(t(i)+h/2, y(i)+(1/2)*k2, v(i)+(1/2)*p2, j(i)+(1/2)*s2)

k4 =  h*F1(t(i)+h, y(i)+k3, v(i)+p3, j(i) + s3)
p4 =  h*F2(t(i)+h, y(i)+k3, v(i)+p3, j(i) + s3, rho(i))
s4 =  h*F3(t(i)+h, y(i)+k3, v(i)+p3, j(i) + s3)

y(i+1) = y(i) + (k1+2*k2+2*k3+k4)/6 
v(i+1) = v(i) + (p1+2*p2+2*p3+p4)/6
j(i+1) = j(i) + (s1+2*s2+2*s3+s4)/6
end

[Walter Roberson]

Atmos1.Hm is not defined in the code, so we cannot test.

[Michal Amar]

You are correct, my mistake. I just uploaded it to the question. Thank you!

Answer 1

[Walter Roberson]

v starts out as 0, and due to the acceleration of -9.8 in

F2 = @(t,y,v,j,r) Fd(t,y,v,j,r)/m(j) - 9.8; %ODE 2: dv/dt = Fd/m - g

then

p1 =  h*F2(t(i),y(i),v(i),j(i),rho(i))

goes negative and then

p2 =  h*F2(t(i)+h/2, y(i)+(1/2)*k1, v(i)+(1/2)*p1, j(i)+(1/2)*s1, rho(i))

the 0 initial value plus -9.8/2 gives a negative third parameter for the F2 call (where it is known as v)

F2 passes that negative value to Fd

Fd =  @(t,y,v,j,r) 3*pi*sqrt(j)*mu*v*Ct(t,y,v,j,r)*Cd(t,y,v,j,r);

which passes it to Ct

Ct = @(t,y,v,j,r) 1+ 0.16*R(t,y,v,j,r)^(2/3);

Ct passes it to R

R = @(t,y,v,j,r) r*v*sqrt(j)/mu;      %r is a temporary variable to later be replaced with the rho(y) function

R receives that negative v, and so calculates a negative result

Ct = @(t,y,v,j,r) 1+ 0.16*R(t,y,v,j,r)^(2/3);

and that negative result is raised to the (2/3) . But negative raised to a fraction is imaginary.

[Michal Amar]

Wow. I truly appreciate this, your explanation is really helpful.

I am not sure how to go about solving this problem, though. V sill always be negative as the projectile is falling in the negative direction. I am thinking that changing my axis system may be very complicated. What do you suggest?

[Walter Roberson]

Is R wind resistance? if so then it acts against the velocity. So -sign(v) * value calculated based on abs(v)