Using a another function within ODE45

I have a function that is part of an ODE, however I have no idea how to get it to work with the ODE45 function. I basically want to pass m_t into disp_vel (at the bottom), for values of t. I have also tried using an anonymous function, however this doesn't seem to work either.

The function m_f should be linear, with grad = fluid_r, until it reaches m_f_max, at which point it is constant.

function HarmonicMotion

%Default Values
    m_c = 2;        % Container Mass
    s1 = 16;        % Spring Constant 1
    s2 = 20;        % Spring Constant 2
    c = 2;          % Damping Coef
    d = 0.6;        % Gap
    t_total = 10;    % Total Integration Time
    dt = 0.001;     % Integration step
    xv0 = [2 0];  % Init Displacement & velocity
    tvals = 0:dt:t_total; % A 1D array for the time values
    m_f_max = 4;    % Max mass of fluid
    fluid_r = 1;    % Fluid mass rate
    data = [m_c, s1, s2, c, d, m_c, fluid_r]; % Setting Array for input into functions
    
    function plot_wf(data) % plots a function with the effects of fluid included
        [t xv2] = ode45(@(t,x) calc_wf(x, data), tvals, xv0);
        figure
        plot(t, xv2(:,:),[0 t_total], [d d], 'k-.');
        xlabel('time');legend('Displacement / m', 'Velocity / ms^{-1}', 'Gap');grid;
        
        
        function [disp_vel] = calc_wf(xv, data) % Calcs values for the case a fluid;
            x = xv(1); v = xv(2); %1st Col = x, 2nd = vel;
            i = 1;
%             m_f(i) = 0;
%             for t = 0:10*dt:4;
%                 if m_f <= m_f_max;
%                     m_f(i+1) = m_f(i) + (fluid_r*(10*dt));
%                 end 
%                 i = i + 1;
%             end
            
             if x > d  %Same as for calc_nf
                 s2_x = x +d;
             else
                 s2_x = 0; 
             end
            disp_vel = [v; -(s1*x + c*v + s2*s2_x)/m_t(1,1)]; % HOW TO GET THIS TO RESPOND TO DIFFERENT VALUES OF m_t
        end
    end
end

[Walter Roberson]

            m_f = @(t)

Invalid syntax. You need to complete the anonymous function

            m_t = m_c + m_f

m_f is a function handle and must be invoked to give a result.

[Fawaz Zaman]

Sorry, I should've edited that part out (I was just testing it out), even when it was completed, it would still not work with ODE45. I've corrected this in the question, thank you.

Answer 1

[Alan Stevens]

More like this (but note the comments near the end):

%Default Values
    m_c = 2;        % Container Mass
    s1 = 16;        % Spring Constant 1
    s2 = 20;        % Spring Constant 2
    c = 2;          % Damping Coef
    d = 0.6;        % Gap
    t_total = 10;    % Total Integration Time
    dt = 0.001;     % Integration step
    xv0 = [2 0];  % Init Displacement & velocity
    tvals = 0:dt:t_total; % A 1D array for the time values
    m_f_max = 4;    % Max mass of fluid
    fluid_r = 1;    % Fluid mass rate
    data = [m_c, s1, s2, c, d, fluid_r ,m_f_max]; % Setting Array for input into functions
    

        [t, xv2] = ode45(@(t,x) calc_wf(t, x, data), tvals, xv0);
        figure
        plot(t, xv2(:,:),[0 t_total], [d d], 'k-.'), grid
        xlabel('time');
        legend('Displacement / m', 'Velocity / ms^{-1}', 'Gap')
        figure
        plot(t,min(fluid_r*t + m_c, m_f_max)), grid
        xlabel('time'),ylabel('fluid mass')
       
        
        function [disp_vel] = calc_wf(t, xv, data) % Calcs values for the case a fluid;
           
        % Extract data
            m_c = data(1); s1 = data(2); s2 = data(3); c = data(4); 
            d = data(5); fluid_r = data(6); m_f_max = data(7);
            
        % Extract displacement and velocity    
            x = xv(1); v = xv(2); %1st Col = x, 2nd = vel;
        
         % Mass at time t  
            m_t = min(fluid_r*t + m_c, m_f_max);  
            
             if x > d  %Same as for calc_nf  %%%%%%%%%% In calc_nf you had x<-d
                 s2_x = x + d;   %%%%%%%%%%%%% Do you mean this? 
                                 %%%%%%%%%%%%% Ifx and d are both positive
                                 %%%%%%%%%%%%% shouldn't you have x-d here?
             else
                 s2_x = 0; 
             end
            disp_vel = [v; -(s1*x + c*v + s2*s2_x)/m_t]; 
        end

[Walter Roberson]

Using min() inside an ode function is usually a problem. The function needs to be twice differentiable, but min() is not differentiable. The if is a problem as well.

You should be using ode event functions to terminate integration at the boundary conditions, and then restarting integration .

[Alan Stevens]

Using min() inside an ode function is sometimes, but not always, a problem. Here, for example, it isn't. Compare the following

The first graph shows the results without a restart; the second shows the results following a restart, (using the pre-restart end values as initial values for the latter stage) after the mass reaches its maximum value.

(In the if statement I've replaced s2_x = x + d; by s2_x = x - d; for both sets of calculations as this should only occur when x exceeds d. This means that there isn't really a significant step change here, it's a gradually increasing one as x gradually exceeds d.)

Nevertheless, Walter is right to point out the potential danger.

[Fawaz Zaman]

Thanks for both your help, you are right it was meant to be (x-d), small oversight on my point - serves me right for working so late :)

I altered the behaviour of m_t to be min(fluid_r*t + m_c, m_f_max + m_c);, since I don't believe that the mass of the container was accounted for previously.

As for the event functions, whilst min() is definitely the easier way to go here, I'll look into them for future reference. Do you happen to know of any good resources to practise / learn these slightly more advanced topics?

Also is there any particular reason that you're extracting the data within the function individually, rather than just using data as a global variable.

as in

        % Extract data
            m_c = data(1); s1 = data(2); s2 = data(3); c = data(4); 
            d = data(5); fluid_r = data(6); m_f_max = data(7);
            

[Walter Roberson]

Global is recommended against, as it is the slowest form of variable and can be the most difficult to debug.

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

With respect to event functions, I recommend the ballode example.

[Alan Stevens]

1. Try

doc ode event location

2. Best to avoid global variables in general. See https://matlab.fandom.com/wiki/FAQ#Are_global_variables_bad.3F for example.