How to treat with the error "Undefined function or variable 'M'." where M is mentioned in a Matlab ode function, and it's also a function of defined variables in the Matalb code?

I have Matlab code and the function that it calls it.

where inside the ode function file:

M=sqrt(Mo-f*t);

where Mo is defined in the code, f is defined in the file of ode function, and t is the time change through ode function.

When I search for the values of M after succeful run on th Commans Window, Matlab gives the error: Undefined function or variable 'M'.

Thanks in advance for each support.

Answer 1

[Jan]

Each function has its own worksapce, whereby "workspace" means the set of locally known variables. If you want to access a value of a variable, which is used inside a function, you have to export this variable as output argument.

[Kyrillos Youssef]

If the variable that I'm asking about is R in the next code, how to make it as an output argument?

%Primary Conditions 
v0 = 10; 
initial_vx=v0;      %initial value for v_x  
initial_vy=0;       %initial value for v_y
initial_x11=0;      %initial value for s_x
initial_y22=0;      %initial value for s_y
xyv0 = [initial_x11, initial_y22, initial_vx, initial_vy];
tspan=[0 4.5];  
[t,xyv]=ode45(@rate,tspan,xyv0);
s_x = xyv(:,1);
s_y = xyv(:,2);
v_x = xyv(:,3);
v_y = xyv(:,4);
plot(t,v_x,'o'),grid  
xlabel('Time (S)') 
ylabel('Vx (m/s)') 
figure 
plot(t,v_y,'o'),grid
xlabel('Time (S)') 
ylabel('Vy (m/s)')
figure
plot(s_x,s_y,'o'),grid
xlabel('Sx (m)') 
ylabel('Sy (m)')
set ( gca, 'ydir', 'r')
function dxyvdt = rate(t,xyv)
       
        v_x = xyv(3);
        v_y = xyv(4);
        rho_p=0.9986*10^3;      %density of object
        rho_f=1.2077;           %density of air 
        v=sqrt(v_x^2+v_y^2);    %object velocity
        g=9.8;
        k1=74*10^-12;           %evaporation constant depends on R_H and T_p0 
        R0=100*10^-6;           %initial object radius  
        vis=1.488*10^-5;         %Kinematic viscosity of air 
        R = sqrt(R0^2 - 2*k1*t);
        Reynolds_no=2*R*v/vis;  % Re depends on current value of R
        c_d = 21.12/Reynolds_no + 6.3/sqrt(Reynolds_no) + 0.25;  %drag coefficient
     
        K = 3*c_d*rho_f/(8*R*rho_p); 
        
        dxyvdt = [v_x;
                  v_y;
                  -K*v_x*v;
                  (1-rho_f/rho_p)*g - K*v_y*v];
        
end

[Jan]

One of these methods are working:

[t, xyv] = ode45(@rate, tspan, xyv0);
k1 = 74e-12;           %evaporation constant depends on R_H and T_p0 
R0 = 100e-6;           %initial object radius  
R  = sqrt(R0^2 - 2 * k1 * t);

or:

[t, xyv] = ode45(@rate,tspan,xyv0);

% Request internal parameters:
[dxyvdt, R] = rate(t.', xyv.');

end

function [dxyvdt, R] = rate(t,xyv)
v_x   = xyv(3, :);        % Accept row vectors also
v_y   = xyv(4, :);
rho_p = 0.9986e3;         %density of object
rho_f = 1.2077;           %density of air 
v     = sqrt(v_x.^2 + v_y.^2);    %object velocity
g     = 9.8;
k1    = 74e-12;           %evaporation constant depends on R_H and T_p0 
R0    = 100e-6;           %initial object radius  
vis   = 1.488e-5;         %Kinematic viscosity of air 
R     = sqrt(R0^2 - 2*k1*t);
Reynolds_no = 2 * R .* v ./ vis;  % Re depends on current value of R
c_d   = 21.12 ./ Reynolds_no + 6.3 ./ sqrt(Reynolds_no) + 0.25;  %drag coefficient
K     = 3 * c_d .* rho_f ./ (8 * R * rho_p);

dxyvdt = [v_x; ...
          v_y; ...
          -K .* v_x .* v; ...
          (1 - rho_f / rho_p) * g - K .* v_y .* v];
end

[Kyrillos Youssef]

All my appreciation for your effort dear friend @Jan. Thanks!

If I did for loop for the value of R0 like:

R0 = [50, 100, 150, 200]*10^-6;  %initial object radius  
tspan=0:0.1:4.5; 
for i = 1:numel(R0)
 
[t,xyv]=ode45(@(t,xyv) ratefor(t,xyv,R0(i)),tspan,xyv0);

s_x(:,i) = xyv(:,1);
s_y(:,i) = xyv(:,2);
v_x(:,i) = xyv(:,3);
v_y(:,i) = xyv(:,4);
end
k1 = 74e-12;           %evaporation constant depends on R_H and T_p0 
R  = sqrt(R0(i)^2 - 2 * k1 * t)

and it successed, But when I demand the vector of R, Matlab shows only what is correponding to the final value of R0 (which is 200). How can I get the vector of R which is corresponding to all the values of R0 (50, 100, 150, 200)?

[Jan]

R  = sqrt(R0 .^ 2 - 2 * k1 * t);  % Auto-expand! >= R2016b

R0 is a row vector and t a column vector. Then the output R is a matrix.

[Kyrillos Youssef]

Great!

[Walter Roberson]

Note that it is common for people to expect to be able to get "all" the values of one of the local variables calculated by an ode function that is being executed under the control of ode45() or similar.

Something like that is possible to get, but it is seldom the right thing to ask for. The ode*() functions evaluate at a whole series of points "near" a proposed point; odePQ evaluates at (P+1) points in order to figure out what the output value should be, and another (Q-P) extra points to cross-check that the prediction works out. If the prediction turns out to be too far out of range, then the step fails and is withdrawn and a new smaller step size is tried. This can happen a number of times in steep areas, so dozens of intermediate locations could have been evaluated at in order for the numeric routines to figure out a single good prediction. Also, points outside the boundaries set by event functions could have been evaluated at.

Furthermore, the ode*() functions do not necessarily report the results of any location they actually evaluated at. Especially when a vector of times was given in the tspan, the ode*() functions make make predictions based upon the current step and tolerances and report those rather than at exact values.

It isn't that it isn't possible to get at all of the values that were evaluated at; it is that it is seldom useful to get at all of those values.

The second approch Jan used, of taking the reported outputs and evaluating the ode function at the reported outputs, turns out to be much more useful in practice.

[Kyrillos Youssef]

Many thamks for your detailed explanation dear @Walter Roberson