Saving derivative values in ode45 in Matlab

Question

I'm simulating equations of motion for a (somewhat odd) system with mass-springs and double pendulum, for which I have a mass matrix and function f(x), and call ode45 to solve

M*x' = f(x,t);

I have 5 state variables, q= [ QDot, phi, phiDot, r, rDot]'; (removed Q because nothing depends on it, QDot is current.) Now, to calculate some forces, I would like to also save the calculated values of rDotDot, which ode45 calculates for each integration step, however ode45 doesn't give this back. I've searched around a bit, but the only two solutions I've found are to a) turn this into a 3rd order problem and add phiDotDot and rDotDot to the state vector. I would like to avoid this as much as possible, as it's already non-linear and this really makes matters a lot worse and blows up computation time.

b) augment the state to directly calculate the function, as described here. However, in the example he says to make add a line of zeros in the mass matrix. It makes sense, since otherwise it will integrate the derivative and not just evaluate it at the one point, but on the other hand it makes the mass matrix singular. Doesn't seem to work for me...

This seems like such a basic thing (to want the derivative values of the state vector), is there something quite obvious that I'm not thinking of? (or something not so obvious would be ok too....)

Oh, and global variables are not so great because ode45 calls the f() function several time while refining it's step, so the sizes of the global variable and the returned state vector q don't match at all.

In case someone needs it, here's the code:

%(Initialization of parameters are above this line)
   options = odeset('Mass',@massMatrix);
   [T,q] = ode45(@f, tspan,q0,options);

function dqdt = f(t,q,p)
    % q = [qDot phi phiDot r rDot]';

    dqdt = zeros(size(q));

    dqdt(1) = -R/L*q(1) - kb/L*q(3) +vs/L;
    dqdt(2) = q(3);
    dqdt(3) = kt*q(1) + mp*sin(q(2))*lp*g;
    dqdt(4) = q(5);
    dqdt(5) = mp*lp*cos(q(2))*q(3)^2 - ks*q(4) - (mb+mp)*g;
end

function M = massMatrix(~,q)
    M = [
        1 0 0 0 0;
        0 1 0 0 0;
        0 0 mp*lp^2 0 -mp*lp*sin(q(2));
        0 0 0 1 0;
        0 0 mp*lp*sin(q(2)) 0 (mb+mp)
        ];
end

Answer 1

The easiest solution is to just re-run your function on each of the values returned by ode45.

The hard solution is to try to log your DotDots to some other place (a pre-allocated matrix or even an external file). The problem is that you might end up with unwanted data points if ode45 secretly does evaluations in weird places.

Answer 2

Since you are using nested functions, you can use their scoping rules to mimic the behavior of global variables.

It's easiest to (ab)use an output function for this purpose:

function solveODE

    % ....        
    %(Initialization of parameters are above this line)

    % initialize "global" variable
    rDotDot = [];

    % Specify output function 
    options = odeset(...
        'Mass', @massMatrix,...
        'OutputFcn', @outputFcn);

    % solve ODE
    [T,q] = ode45(@f, tspan,q0,options);

    % show the rDotDots    
    rDotDot



    % derivative 
    function dqdt = f(~,q)

        % q = [qDot phi phiDot r rDot]';

        dqdt = [...
            -R/L*q(1) - kb/L*q(3) + vs/L
            q(3)
            kt*q(1) + mp*sin(q(2))*lp*g
            q(5)
            mp*lp*cos(q(2))*q(3)^2 - ks*q(4) - (mb+mp)*g
        ];

    end % q-dot function 

    % mass matrix
    function M = massMatrix(~,q)
        M = [
            1 0 0 0 0;
            0 1 0 0 0;
            0 0 mp*lp^2 0 -mp*lp*sin(q(2));
            0 0 0 1 0;
            0 0 mp*lp*sin(q(2)) 0 (mb+mp)
         ];
    end % mass matrix function


    % the output function collects values for rDotDot at the initial step 
    % and each sucessful step
    function status = outputFcn(t,q,flag)

        status = 0;

        % at initialization, and after each succesful step
        if isempty(flag) || strcmp(flag, 'init')
            deriv = f(t,q);
            rDotDot(end+1) = deriv(end);
        end

    end % output function 

end 

The output function only computes the derivatives at the initial and all successful steps, so it's basically doing the same as what Adrian Ratnapala suggested; re-run the derivative at each of the outputs of ode45; I think that would even be more elegant (+1 for Adrian).

The output function approach has the advantage that you can plot the rDotDot values while the integration is being run (don't forget a drawnow!), which can be very useful for long-running integrations.