ODE45 Related Question
3 views (last 30 days)
Show older comments
Ryan Griffith
on 26 Feb 2021
Edited: Ryan Griffith
on 1 Mar 2021
I am having trouble implenting ode45 in a way that covers a time span of 66 minutes for two equations. The equations are nearly identical except the second is a slighly reduced form of the first. The difference lies in essentially two terms within this function. The first is (P_B - v*t) and the second is (r*v/3). I need (P_B - v*t) to reduce to P_F from the interval 6 to 66 minutes and (r*v/3) to go to zero from the same interval. What can I do to combine these two equations such that they are evaluated from their respective intervals (the non-reduced form from 0 to 6 minutes and the reduced form from 6 to 66 minutes)?
clc; clear; close all;
alpha = 0.0125; %Ostwald N2 solubility
D = 2*10^-8; %Diffusion coefficient (cm^2/sec)
h = 3*10^-4; %Bubble thickness (dyne/cm)
H = 2.5*10^8; %Bulk modulus (dyne/cm^2)
P_B = 14.7*6.8947*10^4; %Initial ambient pressure (dyne/cm^2)
P_F = 4.3*6.8947*10^4; %Final ambient pressure (dyne/cm^2)
gamma = 30; %Surface tension (dyne/cm^2)
V_tiss = 1; %Tissue volume (cm^3)
M = H/V_tiss; %Modulus of elasticity (dyne/(cm^2*cm^3))
P_N2o = 11.6*6.8947*10^4; %N2 initial tension (dyne/cm^2)
P_met = 1.76*10^5; %Metabolic gas tension (dyne/cm^2)
r = 0.0003; %Nucleation radius (cm)
t_a = 6*60; %Ascent time
v = (P_B - P_F)/t_a; %Ascent rate
P_N2 = P_N2o-P_N2o*(1-exp(-.00192*120));
tspan = 0:.1:6;
f = P_B - v*tspan;
plot(tspan, f)
tspan = [0 6*60];
r0 = r;
[t,r] = ode45(@(t,r)((-alpha*D/h)*(P_B-v*t+ 2*(gamma/r) + (4/3)*pi*M*r^3 ...
-P_N2 - P_met) + r*v/3)/(P_B-v*t + (4/3)*gamma/r + (8/3)*pi*M*r^3), tspan, r0);
BGI = r/r0;
plot(t/60,BGI)
grid on
xlabel('Time After Ascent (min)')
ylabel('BGI')
ylim([0 16])
tspan = [6*60 66*60];
r0 = r(41);
[t,r] = ode45(@(t,r)((-alpha*D/h)*(P_F + 2*(gamma/r) + (4/3)*pi*M*r^3 ...
-P_N2 - P_met))/(P_F + (4/3)*gamma/r + (8/3)*pi*M*r^3), tspan, r0);
BGI = r/0.0003;
plot(t/60,BGI)
grid on
xlabel('Time After Ascent (min)')
ylabel('BGI')
ylim([0 16])
0 Comments
Accepted Answer
Cris LaPierre
on 27 Feb 2021
Edited: Cris LaPierre
on 1 Mar 2021
I'm no expert, but I think if you create your own odefun, you can use an if statement to change the equation you use based on the value of t. Only one timepoint is evaluated at a time.
I used a nested function approach to avoid having to change your code too much. This gives the odefun access to the constants. I do not know, however, if the results are what you want.
sim
% created nested functions
function sim()
alpha = 0.0125; %Ostwald N2 solubility
D = 2*10^-8; %Diffusion coefficient (cm^2/sec)
h = 3*10^-4; %Bubble thickness (dyne/cm)
H = 2.5*10^8; %Bulk modulus (dyne/cm^2)
P_B = 14.7*6.8947*10^4; %Initial ambient pressure (dyne/cm^2)
P_F = 4.3*6.8947*10^4; %Final ambient pressure (dyne/cm^2)
gamma = 30; %Surface tension (dyne/cm^2)
V_tiss = 1; %Tissue volume (cm^3)
M = H/V_tiss; %Modulus of elasticity (dyne/(cm^2*cm^3))
P_N2o = 11.6*6.8947*10^4; %N2 initial tension (dyne/cm^2)
P_met = 1.76*10^5; %Metabolic gas tension (dyne/cm^2)
r = 0.0003; %Nucleation radius (cm)
t_a = 6*60; %Ascent time
v = (P_B - P_F)/t_a; %Ascent rate
P_N2 = P_N2o-P_N2o*(1-exp(-.00192*120));
tspan = 0:.1:6;
f = P_B - v*tspan;
plot(tspan, f)
tspan = [0 66*60];
r0 = r;
[t,r] = ode45(@(t,r) BGIfun(t,r), tspan, r0);
BGI = r/r0;
plot(t/60,BGI)
grid on
xlabel('Time After Ascent (min)')
ylabel('BGI')
ylim([0 16])
% Nested function
function dydt= BGIfun(t,r)
% select equation to use based on time value t
if t <= 6*60
dydt = ((-alpha*D/h)*(P_B-v*t+ 2*(gamma/r) + (4/3)*pi*M*r^3 ...
-P_N2 - P_met) + r*v/3)/(P_B-v*t + (4/3)*gamma/r + (8/3)*pi*M*r^3);
else
dydt = ((-alpha*D/h)*(P_F + 2*(gamma/r) + (4/3)*pi*M*r^3 ...
-P_N2 - P_met))/(P_F + (4/3)*gamma/r + (8/3)*pi*M*r^3);
end
end
end
More Answers (0)
See Also
Categories
Find more on Mathematics and Optimization in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!