How to extract variables within a ODE function?

Hi Everyone,

I'm solving an ODE using ode45 for variables [Wt2,X] in the following script. In addistion to the answer, I also need few more variables which are calculated within the function. I need to save few variales such as 'Old_U_liq', 'real_rb', etc. Is there any way to save them in the Workspace?

Thank you for your help.

--------------------------------------------------------------------------------------------------

Wspan = [W(1) W(1001)];
X0 = [0 0];                            
[Wt2,X] = ode45(@(Wt2,X) myode45function(Wt2,X,q,rho_l,rho_g,mu_l,g,sigma,theta,S0,U0,re_second_stage,Re_d0,Fr,Eo,V_cap,d0), Wspan, X0);

function dX = myode45function(Wt,X,q,rho_l,rho_g,mu_l,g,sigma,theta,S0,U0,re_second_stage,Re_d0,Fr,Eo,V_cap,d0) 
Fm = S0*rho_g*U0^2;
real_rb = ((Wt+V_cap-pi*(re_second_stage*sind(theta))^2*X(1))*3/(4*pi))^(1/3); 
eps = 1e-15;
a = 0; 
b = theta; 
xm = (a+b)/2;
cont = 0;
fxm = f(xm,real_rb,V_cap); 
while (abs(fxm)>eps)
    fa = f(a,real_rb,V_cap);
    fxm = f(xm,real_rb,V_cap);
    if (fa*fxm > 0)
        a = xm;
    end
    if (fa*fxm < 0)
        b = xm;
    end
    xm = (a+b)/2;
    cont = cont +1;
end
alpha = xm;
drreal_dt = ((q/1000)-pi*(re_second_stage*sind(theta))^2*X(2)*(q/1000))/(4*pi*real_rb^2);
Ub_line = drreal_dt*cosd(alpha)+X(2)*(q/1000);
a = 0.5562; b = -0.04201; c = 0.2337;
coef = a*Fr^b*Eo^c;
Old_U_liq= 4*(q/1000)/(pi*(2*real_rb)^2)*1.6*coef;
Ub_corrected = Ub_line-Old_U_liq;

Fi_part = 11*rho_l/16*(q/1000)*Ub_corrected - 11*rho_l*Wt*(q/1000)*drreal_dt/(32*pi*real_rb^3)*cosd(alpha) +11*rho_l*Wt*(re_second_stage*sind(theta))^2*X(2)*(q/1000)*drreal_dt/(32*real_rb^3)*cosd(alpha);
Fi_part_division = 11*rho_l*Wt/16 -11*rho_l*Wt*(re_second_stage*sind(theta))^2/(64*real_rb^2)*cosd(alpha);
dX = zeros(2,1);
dX(1) = X(2);
dX(2) = (Fm-Fi_part)/(Fi_part_division*(q/1000)^2);
end

--------------------------------------------------------------------------------------------------

Answer 1

[Jan]

The standard method is to add them to the outputs and call the function to be integrated for the time steps afterwards:

function [dX, Old_U_liq, real_rb] = myode45function( ...    

The while-loop looks, like the funtion to be integrated is not smooth. Is this the case? If so, remember, that Matlab's ODE integrators are designed to to interate smooth functions only. Jumps confuse the step size estimator such that you get an error message, or with less luck the step sizes are chosen such small that a huge number of steps increases the accumulated rounding errors until they dominate the calculated trajectory. From the viewpoint of numerical maths, this is an expensive random number generator only and not a scientific result.