Unexpected Imaginary Numbers in Runge Kutta Method Code
4 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Michal Amar
le 27 Juin 2021
Commenté : Walter Roberson
le 28 Juin 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 commentaires
Réponse acceptée
Walter Roberson
le 27 Juin 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 commentaires
Walter Roberson
le 28 Juin 2021
Is R wind resistance? if so then it acts against the velocity. So -sign(v) * value calculated based on abs(v)
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!