Unexpected Imaginary Numbers in Runge Kutta Method Code

10 Ansichten (letzte 30 Tage)
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 Kommentare
Walter Roberson
Walter Roberson am 27 Jun. 2021
Atmos1.Hm is not defined in the code, so we cannot test.
Michal Amar
Michal Amar am 27 Jun. 2021
You are correct, my mistake. I just uploaded it to the question. Thank you!

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Walter Roberson
Walter Roberson am 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 Kommentare
Michal Amar
Michal Amar am 27 Jun. 2021
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
Walter Roberson am 28 Jun. 2021
Is R wind resistance? if so then it acts against the velocity. So -sign(v) * value calculated based on abs(v)

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Numerical Integration and Differential Equations finden Sie in Help Center und File Exchange

Produkte


Version

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by