ODE45 Related Question

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])

Answer 1

[Cris LaPierre]

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

[Ryan Griffith]

Thank you! That works!