Unexpected Imaginary Numbers in Runge Kutta Method Code
4 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
Categories
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!