Unexpected Imaginary Numbers in Runge Kutta Method Code
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Michal Amar
am 27 Jun. 2021
Kommentiert: Walter Roberson
am 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 Kommentare
Akzeptierte Antwort
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
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)
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Numerical Integration and Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!