How do I extract extra parameters from the ode45 ordinary differential equation solver function?

I have the following function which represents an ever increasing mass on a spring according to Hooke's Law:

function dydt = HookeWithRisingMass(t, y, m0, k, Cload)

%HookeWithRisingMass defines Hooke's Law but with a mass increasing at a constant rate over time.

%INPUTS:

%t = [start, finish], time interval

%y = [disp; velocity], initial conditions of the system; displacement from equilibrium and starting velocity

%k = the spring constant

%m0 = initial mass
%Cload = constant of mass increase

%OUTPUTS:

%dydt = [time, [diplacement, velocity]]

%See also

mass = m0 + Cload*t; %mass increases linearly with time

disp = y(1);    %Initial displacement

dxdt = y(2);    % set dx/dt = velocity

dvdt = -(k/mass).*disp;  % set dv/dt = acceleration

dydt = [dxdt; dvdt];  % return the derivatives

It works correctly when called by e.g. [timeout,out] = ode45(@HookeWithRisingMass, tin, y, [], m0, k, Cload).

However I really need to get a time series of values of the mass vs. time for use elsewhere. At the moment I can do this by adding the following to the HookeWithRisingMass function definition:

paff = 'C:\matlabtemp\';

fileName = num2str(t);

fileName(regexp(fileName, '\D')) = 'D';

fileName = ['mass' fileName];

save([paff fileName], 'm', 't');

I can then sequentially load the saved files and store the values of m, t in an array. This is a slow method and gets very unwieldy when the maximum value of time, t, is large. Is there a method of extracting t, m, pairs and storing them in an array that is available in the workspace of the script/function that makes the ode45 call?

Answer 1

[michio]

One way is to use OutputFcn. I have attached an example function that will save mass to the workspace when the calculation is done. Please have myOutputFcn.m saved in the same directly with HookeWithRisingMass.m, then the following script would get you mass.

m0 = 1;
k = 1;
Cload = 1;
y0 = [1,1];
tin = [0,1];
odeFcn = @(t,y) HookeWithRisingMass(t,y, m0, k, Cload);
options = odeset('OutputFcn',@(t,y,flag) myOutputFcn(t,y,flag, m0, Cload));
[timeout,out] = ode45(odeFcn, tin, y0, options);

For details on this option, please see "OutputFcn" section of odeset. The following is the content of myOutputFcn.m.

function status = myOutputFcn(t,y,flag, m0, Cload)
% OutputFcn sample
persistent mass

switch flag
    case 'init'
        mass = m0 + Cload*t; %mass increases linearly with time
    case ''
        mass = [mass, m0 + Cload*t];
    case 'done' % when it's done
        assignin('base','mass',mass); % get the data to the workspace.
end

status = 0;

[Robert]

Well, this works, thanks! All I've got to do now is understand it well enough to be able to modify it for other situations...

[Robert]

Looking at this, I should be able to replace mass = m0 + Cload*t with mass = any explicit function of time but it's going to fail if I use e.g. mass = m0 + rand(1) isn't it?

[michio]

I think mass = m0 + rand(1) and any functions should work as long as the function is defined and it uses variables defined inside myOutputFcn. In the above case, t, y, flag, m0, and Cload.

Let me know if you have any issues.

[Vincent DUFOUR]

I am sorry to re-open this subject, but once you have done it, how do you get to the variable you extracted (I mean in the workspace) ?

Thanks in advance

[Simon Silge]

That is done by the assignin('base','mass',mass)