Runge Kutta and Binary Search to analyze system of ODEs
3 次查看(过去 30 天)
显示 更早的评论
The objective of this code is to solve a system of ODE's representing the motion of a raindrop with inclusion of drag, air density, and mass loss. Within the for loop I used the Runge Kutta method to solve for the position, velocity, and diameter^2 of the drop with respect to time.
Within the while loop lies my issue: I have to find the maximal initial diameter (found in the d vector) such that the drop evaporates before hitting the ground. That is to say, the final position will be greater than zero, for j = 0 (j = d^2). To do this, I am trying to use a binary search algorithm to converge to the diameter I need. I have tested this code with numerous variations (changing n, h, y(1)) in hopes of isolating where my problem is, to no avail. Any critiques, help, suggestions would be greatly appreciated.
clear all
close all
clc
format longE
load('Atmos1.mat');
c = 2.12;
b = 2.28; %Defining constants
a = 0.013;
sigma = 0.07;
K = 6.67e-4;
mu = 1.89e-5;
rho_w = 1000;
g = 9.8;
error = 1e-5;
h = 0.1; %Arbitrary step size
n = 0:h:4000; %Analyzing up to t = some value
d = 0.0005:0.0001:0.005; %Array of initial diameter values
t = zeros(1,1);
v = zeros(1,1);
y = zeros(1,1);
j = zeros(1,1);
rho = zeros(1,length(t));
v(1) = 0; %Initial values for the IVP
y(1) = 400; %I would like to have y(1) = 4500, but for testing purposes it is an arbitrary smaller value
t(1) = 0;
%j(1) = d(1)^2;
low = 1;
high = length(d) - 1;
mid = round((high + low)/2)
d0 = d(mid);
j(1) = d0^2; %j = diameter^2
while (low <= high) %start of binary search
for i = 1:1:length(n) -1 %start of runge kutta method
t(i) = n(i);
i;
if (y(i)<0)
y(i)=0;
break
elseif(j(i) < 0)
j(i) = 0;
break;
elseif(isnan(y(i)) == 1 || isinf(y(i)) == 1) %this often happens, I would love to know why so I can erradicate this condition
y(i) = y(i-1);
break
end
x0 = Atmos1.Hm(60); %linear interpolation equation to define rho(y) function based on given dataset
p = @(x) (x - x0)/(Atmos1.Hm(2) - Atmos1.Hm(1));
f = @(x) p(x)*(p(x)-1)/2;
g = f(y(i))*(p(y(i) - 2))/3;
rho(i) = Atmos1.rhosi(60) + p(y(i))*(Atmos1.rhosi(61) - Atmos1.rhosi(60))...
- f(y(i))*(Atmos1.rhosi(62) -2*Atmos1.rhosi(61)+ Atmos1.rhosi(60))+ g*(Atmos1.rhosi(63) - 3*Atmos1.rhosi(62) + 3*Atmos1.rhosi(61) - Atmos1.rhosi(60));
m = rho_w*(1/6)*pi*(sqrt(j(i)))^3; %mass of raindrop at time i
F1 = @(t,y,v,j) v ; %ODE 1: dy/dt = v
F2 = @(t,y,v,j) (1/m)*3*pi*sqrt(j)*mu*abs(v)... %ODE 2: dv/dt = Fd/m - g
*(1+0.16*(rho(i)*abs(v)*sqrt(j)/mu)^(2/3))*(1 + ((a*(rho(i)*v^2)*sqrt(j)/sigma) + b)^c - a*b^c) - 9.8;
F3 = @(t,y,v,j) -K*(d0)^2; %ODE 3: dj/dt = -Kdo^2
%Runge Kutta coefficients%
k1 = F1(t(i),y(i),v(i),j(i));
p1 = F2(t(i),y(i),v(i),j(i));
s1 = F3(t(i),y(i),v(i),j(i));
k2 = F1(t(i)+h/2, y(i)+(h/2)*k1, v(i)+(h/2)*p1, j(i)+(h/2)*s1);
p2 = F2(t(i)+h/2, y(i)+(h/2)*k1, v(i)+(h/2)*p1, j(i)+(h/2)*s1);
s2 = F3(t(i)+h/2, y(i)+(h/2)*k1, v(i)+(h/2)*p1, j(i)+(h/2)*s1);
k3 = F1(t(i)+h/2, y(i)+(h/2)*k2, v(i)+(h/2)*p2, j(i)+(h/2)*s2);
p3 = F2(t(i)+h/2, y(i)+(h/2)*k2, v(i)+(h/2)*p2, j(i)+(h/2)*s2);
s3 = F3(t(i)+h/2, y(i)+(h/2)*k2, v(i)+(h/2)*p2, j(i)+(h/2)*s2);
k4 = F1(t(i)+h, y(i)+(h)*k3, v(i)+(h)*p3, j(i) + h*s3);
p4 = F2(t(i)+h, y(i)+(h)*k3, v(i)+(h)*p3, j(i) + h*s3);
s4 = F3(t(i)+h, y(i)+(h)*k3, v(i)+(h)*p3, j(i) + h*s3);
y(i+1) = y(i) + h * (k1+2*k2+2*k3+k4)/6 ;
v(i+1) = v(i) + h * (p1+2*p2+2*p3+p4)/6 ;
j(i+1) = j(i) + h * (s1+2*s2+2*s3+s4)/6 ;
%t(i+1) = n(i);
end
figure(1) %Plot to properly see what is happening with each iteration. Goal: Find maximal initial diameter s.t.
plot([j(1) j(i)],[y(1) y(i)]); % the line of motion hits the y axis above 0. Meaning j = 0, y > 0.
title('Position Vs. Diameter^2');
xlabel('d^2 [m^2]');
ylabel('Position [m]');
hold on
mid = round((high + low)/2)
val1 = j(i) %simply just to see what the values are
val2 = y(i)
%Binary search logic%
if ( y(i) < error && j(i) == 0)
fprintf('case 1'); %Just to see what condition is met
break;
elseif (y(i) > error && j(i) >= 0)
fprintf('case 2');
low = mid
elseif (y(i) < error && j(i) > 0)
fprintf('case 3');
high = mid
end
d0 = d(mid)
j(1) = d0^2; %new initial diameter to start the RK method again
end
hold off
0 个评论
采纳的回答
Alan Stevens
2021-7-2
At the moment you have
dj/dt = -K*d0^2;
which means rate of reduction of radius is constant.
Since d0^2 is a constant this means
j(t) = d0^2*(1-K*t)
So j will go negative when t = 1/K.
Perhaps you should have
dj/dt = -K*j^2;
which means rate of reduction of radius is proportional to surface area.
This has the solution
j(t) = d0^2/(1+d0^2*K*t)
This approaches zero asymptotically at infinity.
Neither of these depend on the atmospheric conditions! Perhaps a different model is needed.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!