Unexpected Imaginary Numbers in Runge Kutta Method Code
2 views (last 30 days)
Show older comments
Michal Amar
on 27 Jun 2021
Commented: Walter Roberson
on 28 Jun 2021
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
2 Comments
Accepted Answer
Walter Roberson
on 27 Jun 2021
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.
2 Comments
Walter Roberson
on 28 Jun 2021
Is R wind resistance? if so then it acts against the velocity. So -sign(v) * value calculated based on abs(v)
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!