Getting extra parameters from ODE45 and the mystery transpose

Hi,

I've been reading various posts on getting extra parameters from ODEs and not having much luck implementing them. Concerned that I'm trying something too complicated for my Matlab skills, I've gone back to a basic model to try and understand where I'm going wrong. This has raised a couple of questions.

First, the code, the majority of which is borrowed from others. This is the top half of a simple two degree of freedom mass-spring-damper model (without the plotting and analysis parts):

clear
clc

% define time span
t0 = 0; % s. Input signal start time
t1 = 5.115; % s. Input signal end time. 5.115 gives 1024 points; 61.435 gives 12,288 points
dt = 0.005; % time resolution
tvec = t0:dt:t1; % creates a horizontal vector between t0 and t1 that increments by dt

% define chirp input
A = 10; % mm. Input signal peak amplitude
f0 = 0.5; % Hz. Input signal start frequency
f1 = 20; % Hz. Input signal end frequency
g = (f1./f0).^(1./(t1-t0)); % Exponential growth of chirp frequency
i = A.*sin(f0.*(((g.^tvec)-1)./log(g)).*2.*pi); % Ground input displacement - exponential chirp signal
idot = A.*cos(f0.*(((g.^tvec)-1)/log(g)).*2.*pi).*(2.*pi.*f0.*g.^tvec); % Ground input velocity - exponential chirp signal

% set initial conditions
x0 = [0; 0; 0; 0];

% Sprung mass parameters
ms = 540.5; % kg
ks = 41; % N/mm
cs = 1.5; % Ns/mm
% Unsprung mass parameters
mu = 40; % kg
ku = 350; % N/mm
cu = 0.35; % Ns/mm

% Solve model
[T, x] = ode45(@(t,x) Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu), tvec, x0);

And the ODE function is:

function [dx, Fs] = Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu)

i = interp1(tvec, i, t); % this is interpolating i at t
idot = interp1(tvec, idot, t); % this is interpolating idot at t
% as Matlab is using unspecified time steps it needs a value of i for
% each t

% x(1) is sprung mass displacement
% x(2) is unsprung mass displacement
dx(1) = x(3); % sprung mass velocity. This is first column of dx
dx(2) = x(4); % unsprung mass velocity. This is second column of dx, and so on....
dx(3) = 1./ms.*(ks.*(x(2)-x(1)) + cs.*(x(4)-x(3))).*1000; % sprung mass acceleration
dx(4) = 1./mu.*((ku.*(i(1)-x(2))) + (cu.*(idot(1) - x(4))) - (ks.*(x(2)-x(1))) - (cs.*(x(4)-x(3)))).*1000; % unsprung mass acceleration

dx = dx'; % transpose results from horizontal to vertical
end

This runs fine but the first question is why the need to do the transpose at the bottom when I don't see that line in anyone else's examples? Without it, the code throws an error, and I'm aware that the ODE must return the results in columns. However, I'm confused as to why I don't see the transpose in the help pages or here on Answers.

I then adjusted the code to see if I could get an extra parameter out of the ODE - just a simple dummy parameter as an example. I inserted the following in the ODE:

Fs = ks.*(x(2)-x(1)) + cs.*(x(4)-x(3)); % sum of forces on sprung mass

And adjusted two of the ODE statements to accept the new parameter:

dx(3) = 1./ms.*(Fs).*1000; % sprung mass acceleration
dx(4) = 1./mu.*((ku.*(i(1)-x(2))) + (cu.*(idot(1) - x(4))) - Fs).*1000; % unsprung mass acceleration

I adjusted the function declaration to include the new parameter:

function [dx, Fs] = Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu)

I then asked for the parameter from the ODE:

[dx, Fs] = Two_DOF_QCM_Basic_ODE(T, x, i, idot, tvec, ms, ks, cs, mu, ku, cu);

It ran without error but I only got Fs back with a single value in it, rather than a value for each element of T. How do I get Fs back at all times in T?

Thanks, Simon.

Answer 1

[Mischa Kim]

Hi Simon, to your first question: When you assign values to dx(1), dx(2), and so on you are creating a row vector. However, you need for the ode45 function to return a column vector, hence you have to transpose. What is frequently done in the ode function is to intialize the vector, e.g.

dx = zeros(4,1); % this is a 4x1 column vector

as a column vector. Assigning values to the dx now preserves the column vector. The other approach (which I use frequently) is to assign dx as:

dx    = [x(3); ... % sprung mass velocity. This is first column of dx
         x(4); ... % unsprung mass velocity. This is second column of dx, and so on....
         1./ms.*(ks.*(x(2)-x(1)) + cs.*(x(4)-x(3))).*1000; ... % sprung mass acceleration
         1./mu.*((ku.*(i(1)-x(2))) + (cu.*(idot(1) - x(4))) - (ks.*(x(2)-x(1))) - (cs.*(x(4)-x(3)))).*1000]; % unsprung mass acceleration

which by design is a column vector. To your second question, since you have an equation for Fs you can simply calculate Fs after solving the ODE, i.e. with the return values from the ode45 call.

[Simon Aldworth]

Hi Mischa,

Thanks for your answer - so quick too!

However, with regard to my second question, the code I give is just an example. You are correct that I can calculate Fs with the returned results, but that wouldn't help me understand how to do what I'm asking. The model I am actually trying to build is quite a bit more complicated and I need to learn how to extract the extra parameters out for each of the solved time steps.

I've seen some posts that use outputfcn and others that say don't bother do it another way instead! Unfortunately I understand neither.

Regards.

[Mischa Kim]

The OutputFcn functionality may also be a good approach. In the example below I am computing a randomly made up variable Fs at every integrator time step. You can also append the list of input arguments of myOutputFcn if have additional parameters.

function myPhase()
options = odeset('OutputFcn',@(t,y,flag) myOutputFcn(t,y,flag));
[~,X]   = ode45(@EOM,[0 50],[1 1],options);
u = X(:,1);
w = X(:,2);
plot(u,w)
xlabel('u')
ylabel('w')
grid
end

function status = myOutputFcn(t,y,flag)
persistent Fs

switch flag
    case 'init'
        Fs = y(2);
    case ''
        Fs = [Fs, y(2)*t];
    case 'done' % when it's done
        assignin('base','Fs',Fs); % write data to workspace
end
status = 0;
end

function dX = EOM(~, y)
dX = zeros(2,1);
u  = y(1);
w  = y(2);
A  = 1;
B  = 1;
dX = [w*u^2 - B*u;...
      A - w - w*u^2];
end

[Simon Aldworth]

Moved from below as I placed it in an Answer box instead of a Comment box. Concerned you won't have recevied a notification.

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

Thanks Mischa.

My query differs from your suggestion by the fact that my calculation of the extra parameter I'm interested in, Fs, is within the ODE. I then need to get Fs to the output function. In your case, Fs is calculated within the output function and is only reliant on t and y.

In my case Fs uses two constants that were passed to the ODE and two variables that are calculated within the ODE.

I assume I need to either pass Fs to the output function once it's been calculated in the ODE, or pass the constants to the output function and calculate Fs there, whichever is the most efficient.

Regards.

[Mischa Kim]

Hi Simon, I not sure I understand your requirements and use case. To summarize:

1. The ODE solver (in your case this would be ode45 with function Two_DOF_QCM_Basic_ODE) calls the OutputFcn after each successful integration time step.
2. All variables (e.g. time vector t and the vector that is being integrated -- the state vector) and additional parameters that you pass to the ode function you can also pass to the output function.
3. Hence, just in terms of computing your variable Fs, it makes no difference if you do this in the ode function or the output function.
4. You cannot pull out extra variables directly from the ode function. This is what the output function is for.
5. Bottom line, you can either compute Fs in the ode function and after the integration is done, as I pointed out in my first answer. Or, you do it in the ode function and the output function.

[Simon Aldworth]

Hi Mischa,

Thanks but I should have been clearer by posting the problem I am running into trying to implement OutputFcn.

If you recall, in my first post I tried asking the ODE for Fs after it had solved the equations but I only got a single value back. I then tried to implement OutputFcn and now get errors with Fs being undefined. Presumably this is an issue with the scope of Fs.

To be precise: the code in my main model is now arranged with Fs being passed to OutputFcn:

% Set up ODE function to pass to solver
odeFcn = @(t, x) Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu);
% Set up options for output function, passing all variables (because I can)
options = odeset('OutputFcn',@(t, x, flag) myOutputFcn(t, x, flag, Fs));

% Solve model
[T, x] = ode45(odeFcn, tvec, x0, options);

Fs is calculated as before in the ODE:

Fs = ks.*(x(2,:)-x(1,:)) + cs.*(x(4,:)-x(3,:)); % sum of forces on sprung mass

% x(1) is sprung mass displacement
% x(2) is unsprung mass displacement
dx(1) = x(3,:); % sprung mass velocity. This is first column of dx
dx(2) = x(4,:); % unsprung mass velocity. This is second column of dx, and so on....
dx(3) = 1./ms.*(Fs).*1000; % sprung mass acceleration
dx(4) = 1./mu.*((ku.*(i(1,:)-x(2,:))) + (cu.*(idot(1,:) - x(4,:))) - Fs).*1000; % unsprung mass acceleration

And finally myOutputFcn receives Fs and computes the persistant variable Fs_out:

function status = myOutputFcn(t, x, flag, Fs)
% OutputFcn sample
persistent Fs_out

switch flag
    case 'init'
        Fs_out = Fs;
    case ''
        Fs_out = [Fs_out, Fs];
    case 'done' % when it's done
        assignin('base','Fs',Fs_out); % get the data to the workspace.   
end

status = 0;

However, when I run the model and it gets to the line to solve using ode45 I receive an error saying that Fs is undefined:

Undefined function or variable 'Fs'.

Error in Two_DOF_Quarter_Car_Model_Basic>@(t,x,flag)myOutputFcn(t,x,flag,Fs)

Error in ode45 (line 269)

feval(outputFcn,[t tfinal],y(outputs),'init',outputArgs{:});

So, I suppose the question is how do I pass Fs to OutputFcn ?

Many thanks.

[Mischa Kim]

According to the documentation (see OutputSel) you can only pass state vector components (the solution) from the ode solver to the output function. You cannot pass other variables, e.g. Fs. This means you need to compute Fs twice. Once in the ode function (to solve the differential equation) and again in the output function.

Here you go:

clear
clc

% define time span
t0   = 0; % s. Input signal start time
t1   = 5.115; % s. Input signal end time. 5.115 gives 1024 points; 61.435 gives 12,288 points
dt   = 0.005; % time resolution
tvec = t0:dt:t1; % creates a horizontal vector between t0 and t1 that increments by dt

% define chirp input
A    = 10; % mm. Input signal peak amplitude
f0   = 0.5; % Hz. Input signal start frequency
f1   = 20; % Hz. Input signal end frequency
g    = (f1./f0).^(1./(t1-t0)); % Exponential growth of chirp frequency
i    = A.*sin(f0.*(((g.^tvec)-1)./log(g)).*2.*pi); % Ground input displacement - exponential chirp signal
idot = A.*cos(f0.*(((g.^tvec)-1)/log(g)).*2.*pi).*(2.*pi.*f0.*g.^tvec); % Ground input velocity - exponential chirp signal

% set initial conditions
x0 = [0; 0; 0; 0];

% Sprung mass parameters
ms = 540.5; % kg
ks = 41; % N/mm
cs = 1.5; % Ns/mm
% Unsprung mass parameters
mu = 40; % kg
ku = 350; % N/mm
cu = 0.35; % Ns/mm

% Set up ODE function to pass to solver
odeFcn  = @(t,x) Two_DOF_QCM_Basic_ODE(t,x,i,idot,tvec,ms,ks,cs,mu,ku,cu);
% Set up options for output function, passing all variables (because I can)
options = odeset('OutputFcn',@(t,x,flag) myOutputFcn(t,x,flag,cs,ks));
% Solve model
[T, x]  = ode45(odeFcn,tvec,x0,options);
subplot(2,1,1)
plot(T,x)
subplot(2,1,2)
plot(T,Fs)

function dx = Two_DOF_QCM_Basic_ODE(t,x,i,idot,tvec,ms,ks,cs,mu,ku,cu)
i    = interp1(tvec, i, t); % this is interpolating i at t
idot = interp1(tvec, idot, t); % this is interpolating idot at t
% as Matlab is using unspecified time steps it needs a value of i for
% each t

Fs    = ks.*(x(2,:)-x(1,:)) + cs.*(x(4,:)-x(3,:)); % sum of forces on sprung mass
% x(1) is sprung mass displacement
% x(2) is unsprung mass displacement
dx(1) = x(3,:); % sprung mass velocity. This is first column of dx
dx(2) = x(4,:); % unsprung mass velocity. This is second column of dx, and so on....
dx(3) = 1./ms.*(Fs).*1000; % sprung mass acceleration
dx(4) = 1./mu.*((ku.*(i(1,:)-x(2,:))) + (cu.*(idot(1,:) - x(4,:))) - Fs).*1000; % unsprung mass acceleration

dx = dx'; % transpose results from horizontal to vertical
end

function status = myOutputFcn(~,x,flag,cs,ks)
persistent Fs_out
switch flag
    case 'init'
        Fs_out = ks.*(x(2,:)-x(1,:)) + cs.*(x(4,:)-x(3,:));
    case ''
        Fs_out = [Fs_out, ks.*(x(2,:)-x(1,:)) + cs.*(x(4,:)-x(3,:))];
    case 'done' % when it's done
        assignin('base','Fs',Fs_out); % get the data to the workspace.   
end
status = 0;
end

[Walter Roberson]

It is common for people to think that they would like to record all calculated intermediate values such as Fs, thinking that this will give a record of the Fs values that were used for each integration step. However, that does not work in practice.

Each "step" for ode45 involves evaluating the given function with 6 different boundary conditions including a mix of different times on the 6 different calls. If the error estimate between the 5th and 6th evaluation is satisfactory, then the step is "accepted". If the error estimate is not satisfactory, then the step is rejected and the same base point is used but with a smaller step size, again making 6 calls. A number of rejections can occur in a row, until it has found a step size that is small enough to accomodate whatever rapid changes in the function are happening.

From time to time, results are output, but not necessarily at any point that was evaluated -- for example instead of evaluating at time 2.0 exactly, it might synthesize and output based upon its calculations at 1.97 and 2.22 and 2.05, extrapolating what the readings for 2.0 would have been.

Now, it is unlikely that you want the value of intermediate variables during calculations for steps that end up being rejected --- the evaluations might even have been at locations that turn out to be outside boundaries (according to event functions.)

And you probably do not want the value of intermediate variables during the 6 different calculations that go into making up one cycle (though sometimes you do). Most of the time, what you want is the value of the intermediate variables associated with the locations actually output... but as I indicated earlier, ode*() might not have evaluated at that exact location!

Because of these factors, most of the time it is better to take the output locations and invoke a function that calculates just the intermediate values for those locations. It could be a function that you would otherwise invoke from inside the ode function.

[Simon Aldworth]

@Mischa Kim, many thanks - that's working fine; and thanks for the lessons. I have also managed to get the alternative method working, whereby the solver is called to complete the integration and then ode is called a second time to provide addtional output at the solved time steps. Thanks again.

@Walter Roberson. Thank you for your explanation; I understand the concept that the integration time steps do not necessarily align to the output steps. This is why my original post was written the way it was. However, so that I fully understand your reasons for posting: are you saying that something we are doing with my code is wrong? Should I be doing something differently? Regards.

Simon.

Answer 2

[Stephen23]

Using the OutputFcn is a complex way to get the Fs values.

The simpler approach is to run the ODE solver normally, and then run the function with the values returned by the ODE solver, which automatically calculates Fs values which correspond exactly to the T and x values returned by the ODE solver. This just takes two lines of code, as shown below:

% define time span
t0 = 0; % s. Input signal start time
t1 = 5.115; % s. Input signal end time. 5.115 gives 1024 points; 61.435 gives 12,288 points
dt = 0.005; % time resolution
tvec = t0:dt:t1; % creates a horizontal vector between t0 and t1 that increments by dt

% define chirp input
A = 10; % mm. Input signal peak amplitude
f0 = 0.5; % Hz. Input signal start frequency
f1 = 20; % Hz. Input signal end frequency
g = (f1./f0).^(1./(t1-t0)); % Exponential growth of chirp frequency
i = A.*sin(f0.*(((g.^tvec)-1)./log(g)).*2.*pi); % Ground input displacement - exponential chirp signal
idot = A.*cos(f0.*(((g.^tvec)-1)/log(g)).*2.*pi).*(2.*pi.*f0.*g.^tvec); % Ground input velocity - exponential chirp signal

% set initial conditions
x0 = [0; 0; 0; 0];

% Sprung mass parameters
ms = 540.5; % kg
ks = 41; % N/mm
cs = 1.5; % Ns/mm
% Unsprung mass parameters
mu = 40; % kg
ku = 350; % N/mm
cu = 0.35; % Ns/mm

% Solve model
fun = @(t,x) Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu);

[T, x] = ode45(fun, tvec, x0) % solve

T = 1024×1

0 0.0050 0.0100 0.0150 0.0200 0.0250 0.0300 0.0350 0.0400 0.0450

x = 1024×4

0 0 0 0 0.0000 0.0085 0.0244 4.3093 0.0004 0.0502 0.1487 12.8568 0.0018 0.1397 0.4273 22.9849 0.0050 0.2786 0.8805 32.2528 0.0109 0.4578 1.4973 38.9086 0.0202 0.6619 2.2440 42.1296 0.0334 0.8734 3.0753 41.9685 0.0510 1.0774 3.9473 39.2069 0.0729 1.2631 4.8229 35.0102

[~,Fs] = cellfun(fun,num2cell(T),num2cell(x,2),'uni',0); % This is all you need...
Fs = cell2mat(Fs) % ... and the Fs values correspond exactly to the T and x values.

Fs = 1024×1

0 6.7739 21.1026 39.4898 58.2753 74.4413 86.1407 92.7790 94.9738 94.0782

function [dx, Fs] = Two_DOF_QCM_Basic_ODE(t, x, i, idot, tvec, ms, ks, cs, mu, ku, cu)

i = interp1(tvec, i, t); % this is interpolating i at t
idot = interp1(tvec, idot, t); % this is interpolating idot at t
% as Matlab is using unspecified time steps it needs a value of i for
% each t

Fs = ks.*(x(2)-x(1)) + cs.*(x(4)-x(3)); % sum of forces on sprung mass

dx = nan(4,1); % !!! preallocate column output !!!

% x(1) is sprung mass displacement
% x(2) is unsprung mass displacement
dx(1) = x(3); % sprung mass velocity. This is first column of dx
dx(2) = x(4); % unsprung mass velocity. This is second column of dx, and so on....
dx(3) = 1./ms.*(Fs).*1000; % sprung mass acceleration
dx(4) = 1./mu.*((ku.*(i(1)-x(2))) + (cu.*(idot(1) - x(4))) - Fs).*1000; % unsprung mass acceleration
end

[Simon Aldworth]

Many thanks for this Stephen, I'll study this in more detail later. For now, though, I'm a bit confused about your statement that the outputfcn approach collects all intermediate points, including when the ODE steps backwards or takes tiny steps. The help page for odeset says that "the ODE solver calls the output function after each successful time step", which to me implies that since I've called the ode with the tspan option, a successful time step will be at each element of tspan.

That has been the premise of my enquiry from the beginning - that I should only get Fs for each successful integration time step and that those steps are defined by tspan.

Regards.

[Stephen23]

"...l'm a bit confused about your statement that the outputfcn approach collects all intermediate points..."

Sorry, that was my mistake (since removed from my answer): I confused the common attempt of using persistent variables or similar to store all intermediate calculations, including those backward steps, etc., with the use of outputfcn.

I would have to check the outputFcn documentation, and probably try some examples.

In any case, it is still the more complex approach, if your goal is just to get the Fs values.