Passing extra output arguments from ode

I've the following code to obtain [t,x] after simulating odes . In addition to the solution, I also want to obtain the functiion values F(x,t) in the below equation.

% obtains the dynamics from a system of odes

x0      =  initial values;
tspan   = 0:0.0001:0.5;
options = odeset('abstol', 1e-10, 'reltol', 1e-9);

% generate  solution
[t, x]  = ode15s(@(t,x) fun(t,x), tspan , x0 ,options);

% obtain function values ,f, from ode
[dy f]  = fun([ ], [ ])

function [dy f] = fun(t,x)
    persistent f
     f  = []
     A  = % Matrix
     dy = A*x
     if nargout >1
         f = horzcat(f, A*x)
      end
end

In the code, I've defined a persistent variable, f, that stores F for each time step that I solver takes. And I receive f using `[dy f] = fun([ ], [ ])``

But the size of f and x are not equl because from what I understand the output of x is given (after interpolation I guess) only at the time points defined in tspan. Whereas, f is saved for very time instant t that the solver takes internally.

In summary my question is, how do I obtain f for the same time instants at which x is obtained.

Answer 1

[Stephen23]

"In summary my question is, how do I obtain f for the same time instants at which x is obtained."

Using a persistent variable or OutputFcn is a red-herring. The simple solution is to call the function after calling the ODE solver, using exactly the values that the ODE solver returns. Here is a complete working example simply based on the first example from the ode15s documentation:

tspan = [0,2];
x0    = 1;
fun   = @(t,x) -10*t;
[t,x] = ode15s(fun, tspan, x0)

Then you just need this:

dy = arrayfun(fun,t,x)

[Deepa Maheshvare]

Thanks a lot for the response. Could you please explain how to use arrayfun when x is a vector (i.e when x = [x1 x2 x3..] like the e.g. provided in my original post)?

[Stephen23]

Either use a loop or cellfun, e.g. using the second example from the ode15s documentation:

>> [t,y] = ode15s(@vdp1000,[0 3000],[2 0]);
>> dydt = cellfun(@vdp1000,num2cell(t),num2cell(y,2),'UniformOutput',false);
>> dydt = [dydt{:}].' % optional
dydt =
            0           -2
   -2.817e-05      -1.9155
  -5.5179e-05      -1.8345
   -8.105e-05      -1.7569
   -0.0001866      -1.4402
  -0.00027362      -1.1791
  -0.00034496     -0.96512
  -0.00040334     -0.78997
... lots of lines here
   0.00068031   1.7023e-05
   0.00068638   2.6636e-05
   0.00072007   1.1815e-05
   0.00075974   0.00012112
   0.00080746   0.00037099
   0.00086676   0.00031263
   0.00095752   0.00041708
     0.001089   0.00047778
    0.0011951   0.00022839

No persistent variables, no changing the function code, no complicated interpolation, no mismatched sample times, no messing around with OutputFcn... just two lines of simple code.

[Deepa Maheshvare]

For instance , if my y has size 500*10

num2cell(y,2) returns

500×1 cell array

    {1×10 double}
    {1×10 double}
    :
    : 
    {1×10 double}

Now the problem is , I want to pass y to function fun as 10 x 1 double for matrix multiplication to be consistent.

[Stephen23]

"Now the problem is , I want to pass y to function fun as 10 x 1 double for matrix multiplication to be consistent."

You could:

• make a small change to fun so that it reshapes the input into a column vector.
• call fun from an anonymous function, which reshapes the input into a column vector.
• use a loop or cellfun to reshape the cell contents, before calling fun via cellfun.

If you showed your code and gave some sample data then I could test these for you too.

[Deepa Maheshvare]

I could implement your first suggestion.

Here I am presenting a complete code

clc;
clear all;
global ngrid

ngrid = 10;
phi0    = 300*ones(ngrid,1);
tspan   = 0:0.01:0.5;
options = odeset('abstol', 1e-10, 'reltol', 1e-9);
D  = 9*ones(ngrid-1,1);

%% generate exact solution
[t, y]  = ode15s(@(t,phi) fun(t,phi,D), tspan , phi0 ,options);

dydt = cellfun(@fun,num2cell(t),num2cell(y,2),num2cell(D),'UniformOutput',false);
dydt = [dydt{:}].' % optional

function dphi = fun(t,phi,D)
 phi
[m,n] = size(phi);
if n ~= 1
    phi = transpose(phi);  % reshapes row vector to a column vector
end

 MT  = [-1  0  0  0   0  0  0  0  0;...
        1  -1  0  0   0  0  0  0  0;...
        0   1  -1 0   0  0  0  0  0;...
        0   0  1  -1  0  0  0  0  0;...
        0   0  0   1  -1 0  0  0  0;...
        0   0  0   0   1 -1 0  0  0;...
        0   0  0   0   0  1 -1 0  0;...
        0   0  0   0   0  0  1 -1 0;...
        0   0  0   0   0  0  0  1 -1;...
        0   0  0   0   0  0  0  0  1];

 M  = transpose(MT);
 A  = MT*diag(D)*M
 
 dphi =  A*phi;

 end

I have a new addition here,

% generate  solution
[t, x]  = ode15s(@(t,x) fun(t,x), D tspan , x0 ,options);

In addition to x and t , fun takes an additional argument D as an input.

In such a case, I am not sure how to use cellfun

I understand

dydt = cellfun(@fun,num2cell(t),num2cell(y,2),num2cell(D),'UniformOutput',false);

will not work

But I am not sure how to pass D in cellfun.. Should I convert D as a 51 x 9 matrix, convert to cell array and then pass? Is there an alternate way?

Please note:

I can't define D inside fun, D comes from an upstream function.

[Stephen23]

For cellfun use exactly the same anonymous function like you did with ode15s:

@(t,phi) fun(t,phi,D)

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

[Deepa Maheshvare]

Yes, I tried the same way. The problem is, D is not a function of time. So the size of D is not equal to size of t or size(x,1). I have to pass D for every t, time step that the cellfun calls the anonymous function fun

[Stephen23]

"The problem is, D is not a function of time.

Sure, that is clear.

So the size of D is not equal to size of t or size(x,1). I have to pass D for every t, time step that the cellfun calls the anonymous function fun"

That is exactly what the anonymous function does.

When fun is parameterized exactly as you did for ode15s then array D is passed complete as an input to fun every time ode15s calls the anonymous function, or every time cellfun calls the anonymous function. The code I showed you does NOT require/assume that D has the same number of elements as t.

You seem a bit confused about parameterizing the function. I recommend that you reassure yourself by using the anonymous functions exactly as I wrote and printing D to the command window (from inside fun) and you will see that in fact both ode15 and cellfun get the complete D array (and so its size is completely unrelated to t or anything else, D can be literally any size or class that you want).

[Deepa Maheshvare]

"When fun is parameterized exactly as you did for ode15s then array D is passed complete as an input to fun every time ode15s calls the anonymous function, or every time cellfun calls the anonymous function"

Yes, I tried this

[t, y]  = ode15s(@(t,phi) fun(t,phi,D), tspan , phi0 ,options);
dydt = cellfun(@fun,num2cell(t),num2cell(y,2),D,'UniformOutput',false);

Error:

Error in ----
f = cellfun(@predicted,num2cell(t),num2cell(phi_tilde,2),Dhat,'UniformOutput',false);

Error using cellfun
Input #4 expected to be a cell array, was double instead.

[Stephen23]

What you did:

cellfun(@fun, ...

What I told you to do:

cellfun(@(t,phi) fun(t,phi,D), ...

Just like I wrote, you need to use exactly the same anonymous function as you did with ode15s. For whatever reason you continue to call fun directly, even though I told you to use that anonymous function.

You could even just define it once (which is what I would do):

D = [...];
afun = @(t,phi) fun(t,phi,D); % define the anonymous function once!
[t,y] = ode15s(afun, ...);
dydt = cellfun(afun, ...);

[Deepa Maheshvare]

Thank you

[Deepa Maheshvare]

I tried the second suggestion

afun = @(t,phi) fun(t,phi,D); % define the anonymous function once!
[t,y] = ode15s(afun, ...);
dydt = cellfun(afun, ...);

Code:

clc;
clear all;
global ngrid

ngrid = 10;
phi0    = 300*ones(ngrid,1);
tspan   = 0:0.01:0.5;
options = odeset('abstol', 1e-10, 'reltol', 1e-9);
D  = 9*ones(ngrid-1,1);

%% generate exact solution

% [t, y]  = ode15s(@(t,phi) fun(t,phi,D), tspan , phi0 ,options);
% f = cellfun(@(t,phi) fun(t,phi,D),'UniformOutput',false)

afun = @(t,phi) fun(t,phi,D); % define the anonymous function once!
[t,y] = ode15s(afun,tspan , phi0 ,options);
f = cellfun(afun,'UniformOutput',false);

f = [f{:}].' % optional

function dphi = fun(t,phi,D)
 phi
[m,n] = size(phi);
if n ~= 1
    phi = transpose(phi);
end

 MT  = [-1  0  0  0   0  0  0  0  0;...
        1  -1  0  0   0  0  0  0  0;...
        0   1  -1 0   0  0  0  0  0;...
        0   0  1  -1  0  0  0  0  0;...
        0   0  0   1  -1 0  0  0  0;...
        0   0  0   0   1 -1 0  0  0;...
        0   0  0   0   0  1 -1 0  0;...
        0   0  0   0   0  0  1 -1 0;...
        0   0  0   0   0  0  0  1 -1;...
        0   0  0   0   0  0  0  0  1];

 M  = transpose(MT);
 A  = MT*diag(D)*M
 
 dphi =  A*phi;

 end

Error:

Error using cellfun
Input #2 expected to be a cell array, was char instead.

Error in simulate (line 18)
f = cellfun(afun,'UniformOutput',false);

[Stephen23]

For some reason you decided to call cellfun without any data input arguments. In your previous comments you showed that you called cellfun correctly with fun and two cell array input arguments (corresponding to t and y). Now you need to do the same with afun: you still need to provide input arguments to cellfun otherwise it cannot provide the data (i.e. t and y) when it calls the anonymous function.

What you did:

dydt = cellfun(afun,'UniformOutput',false); % where are t and y ?

What you need to do:

dydt = cellfun(afun,num2cell(t),num2cell(y,2),'UniformOutput',false);

Earlier I told you to use the anonymous function. I did NOT tell you to remove all of the other inputs!

[Deepa Maheshvare]

Thank you. After correcting the stupid mistake from my side,

line 18 is now:

f = cellfun(afun,num2cell(t),num2cell(y,2),num2cell(D),'UniformOutput',false);

Error:

Error using cellfun
All of the input arguments must be of the same size and shape.
Previous inputs had size 51 in dimension 1. Input #4 has size 9

Error in simulate (line 18)
f = cellfun(afun,num2cell(t),num2cell(y,2),num2cell(D),'UniformOutput',false);
 

[Stephen23]

D is not an input to cellfun or ode15s.

The array D is passed to fun solely via the anonymous function definition, which access D directly from the workspace where the anonymous function is defined. You do NOT need pass D as an input to cellfun.

What you did:

f = cellfun(afun,num2cell(t),num2cell(y,2),num2cell(D),'UniformOutput',false);
%                                          ^^^^^^^^^^^^ does not belong here

What you need to do:

dydt = cellfun(afun,num2cell(t),num2cell(y,2),'UniformOutput',false);
%              ^^^^ D is already in this function definition!!!!!!!

I suggest that you read the link I gave you earlier on how to parameterize functions. Then you would get it working in 20 seconds. Making code changes randomly, without understanding how the anonymous function is used, is not an efficient use of your time:

https://en.wikipedia.org/wiki/Programming_by_permutation

Answer 2

[Ameer Hamza]

You need to use OutputFcn to get the value of a variable at the same time instant as the value in the output vector: https://www.mathworks.com/help/matlab/ref/odeset.html#d120e857290. Following code shows an example

clear global
global x_vector;

opts = odeset('OutputFcn', @outfcn);

[t,x] = ode45(@(t,x) [x(2); -x(1)], [0, 2*pi], [0; 1], opts);

plot(t, x_vector)

function s = outfcn(t,x,f)
    if strcmp(f, 'done')
        return
    end
    global x_vector;
    x_vector = [x_vector x(1, :).^2];
    s = 0;
end

Note that in this code, I used the global variable to get value out of outfcn, to show a general idea. In practice, global variables are not recommended. I suggest you look at John D'Errico's answer here: https://www.mathworks.com/matlabcentral/answers/510713-is-it-possible-to-store-the-intermediate-values-of-fmincon to see a more efficient way to get value out of the outfcn.