ODE45 Related Question

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

Accepted Answer

Cris LaPierre
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)

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!