Runge Kutta and Binary Search to analyze system of ODEs

3 views (last 30 days)
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

Accepted Answer

Alan Stevens
Alan Stevens on 2 Jul 2021
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.

More Answers (0)

Categories

Find more on Numerical Integration and Differential Equations in Help Center and File Exchange

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!