How do I extract an intermediate variable calculated and used inside my ode45 function?

Hello,

I'm using ode45 to produce solutions to ODE's and It works perfectly with me. However, I have a variable that depends on the current differentiated variable ,say x as in my example below. This variable is computed inside my function which then used in my differential equations. This is my ode call:

[T,X] = ode45(@(t,x) myode(t,x,P,BK,R), tspan, x0);

and here is my function:(the variable I need to extract at each instant, i.e. it should have same length as X and T, is u)

function dx=myode(t,x,P,BK,R)
x1=x(1); x2=x(2);
N=size(BK,1);
f1=x2;
f2=sin(x1)
f=[f1;f2]; %system vector field
gx=[0;cos(x1)];

phi=(ones(length(P),1)*x.').^BK(1:length(P),:);
phi=prod(phi,2);

u=-inv(R)*gx'*P'*phi;
if u>=1
    u=1;
end

dx=f+gx*u;
end

Is it possible to get u? what is the best way to do that? I know similar questions were and are being asked a lot but I searched for the solution and couldn't find a workable one to my problem!

For a simpler example (because the variables I have as my inputs above are very large matrices) and we can work on this since the main goal is the same:

x0=[1;1]; tspan=[0 1]; R=1;
[T,X] = ode45(@(t,x) myode(t,x,R), tspan, x0);
    function dx=myode(t,x,R)
    x1=x(1); x2=x(2);
    f1=x2;
    f2=sin(x1);
    f=[f1;f2]; %system vector field
    gx=[0;1];
    u=-inv(R)*gx'*x;
    if u>=1
        u=1;
    end
    dx=f+gx*u;
    end

I want to get the computed values of u at each time instant with t with x. In other words, I want to get an array with the computed values of u.

Thanks in advance. Your help is appreciated.

[madhan ravi]

upload your question(equation) as latex form and provide all the necessary details

[Stephen23]

Note that this line is probably not very efficient or accurate:

u=-inv(R)*gx'*P'*phi;

The inv and * should probably be replaced by mldivide. If you have large arrays then this will be significantly more efficient, and well as being more accurate numerically.

Answer 1

[Stephen23]

"I want to get the computed values of u at each time instant with t with x"

The simplest solution is to add a second output argument u:

function [dx,u] = myode(t,x,P,BK,R)

Then call your function normally to solve:

[T,X] = ode45(@(t,x) myode(t,x,P,BK,R), tspan, x0);

and then, now that you have the final T and X values, calculate the U values:

[~,U] = cellfun(@(t,x)  myode(t,x.',R), num2cell(T), num2cell(X,2),'uni',0);

You could try it without the 'uni' option too, it might work. Otherwise simply convert the cell array to a matrix:

U = cell2mat(U);

[Has MBK]

That would be great if it works! I've tried exactly what you said and I got this error:

Error using arrayfun
All of the input arguments must be of the same size and shape.
Previous inputs had size 1 in dimension 2. Input #3 has size 2

Error in sec_ord_invpend_sakamoto (line 60)
[~,U] = arrayfun(@(t,x) myode(t,x,P,BK), T, X);

I'm not sure what meant in the error about the inputs dimentions. I even removed R so I can have same size and shape matrices (P and BK) but I got the same error. Let's us assume my P and BK are 132x2 matrices. Any thoughts? Thanks again.

[Stephen23]

@Has MBK: I forgot about the dimensions of X. See my edited answer using num2cell.

I tried this on your simple example function:

>> [T,X] = ode45(@(t,x) myode(t,x,R), tspan, x0)
T =
   0.00000
   0.06813
   0.16813
   0.26813
   0.36813
   0.46813
   0.56813
   0.66813
   0.76813
   0.86813
   0.96813
   1.00000
X =
   1.00000   1.00000
   1.06780   0.99074
   1.16639   0.98199
   1.26435   0.97775
   1.36205   0.97671
   1.45975   0.97769
   1.55762   0.97960
   1.65567   0.98144
   1.75387   0.98230
   1.85207   0.98138
   1.95006   0.97794
   1.98120   0.97622
>> [~,U] = cellfun(@(t,x)  myode(t,x.',R), num2cell(T), num2cell(X,2),'uni',0);
>> cell2mat(U)
ans =
  -1.00000
  -0.99074
  -0.98199
  -0.97775
  -0.97671
  -0.97769
  -0.97960
  -0.98144
  -0.98230
  -0.98138
  -0.97794
  -0.97622

The only change I made to your ode function was to add u to the output arguments:

function [dx,u] = myode(t,x,R)

[Has MBK]

IT WORKS PERFECTLY! Thanks a lot Stephen. I highly appreciate your help.

Answer 2

[madhan ravi]

https://www.mathworks.com/matlabcentral/answers/92701-how-do-i-pass-out-extra-parameters-using-ode23-or-ode45-from-the-matlab-ode-suite#answer_102051

[Has MBK]

Thanks for your response. You get this error because you do not have the variables I'm feeding my ode like BK, P, etc. Take this simple example (we can work on this since the main goal is the same):

x0=[1;1]; tspan=[0 1]; R=1;
[T,X] = ode45(@(t,x) myode(t,x,R), tspan, x0);

    function dx=myode(t,x,R)
    x1=x(1); x2=x(2);

    f1=x2;
    f2=sin(x1);
    f=[f1;f2]; %system vector field
    gx=[0;1];

    u=-inv(R)*gx'*x;
    if u>=1
        u=1;
    end

    dx=f+gx*u;
    end

I want to get the computed values of u at each time instant with t and x. In other words, I want to get an array with the computed values of u.

[Has MBK]

Thanks a lot for your help but I've tried this and it didn't work well with me. Maybe I have not done it correctly! though I have followed the provided steps. If there is a clever way to save each new u in an array, we can then save it to an external file but I'm not able to figure out such a method because at each instant in the ode, it overwrites everything.

[madhan ravi]

See sir Walters answer

Answer 3

[Walter Roberson]

If the variable is to be the same length as t then you need to recalculate it afterwards. The ode* routines do not generally calculate at the returned t: they calculate at multiple times and use them to estimate the values that they return. In general the ode* routines may go forwards or backwards in time and may calculate with multiple boundary conditions at the same time point in order to be sure that they are within the error constraints.

[Has MBK]

Thanks for your response. The exact length doesn't matter but I need to get the used u's to compute dx. I've come over a method that I save the u and t by adding a save command in my function

save('test.mat', 'u', 't');

like the method mentioned here: https://www.mathworks.com/matlabcentral/answers/305698-how-do-i-extract-extra-parameters-from-the-ode45-ordinary-differential-equation-solver-function

However, this saves u and t one by one and my ode generates thousands of them so it is very inefficient and tedious work. Also, in the link the problem was solved using " OutputFcn" but I've tried that and it gives me the following errors

Error using exist The first input to exist must be a string scalar or character vector.

Error in odearguments (line 59) if (exist(ode)==2)

Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in sec_ord_invpend_sakamoto (line 60) [timeout,out] = ode45([T,X] , tspan, x0, options);

By the way, you said that I need to recalculate it afterwards. What is the best (accurate) way to do that? I've used the generated X to compute u in a for loop but I don't know why for complex or large values, it seems to not work right! It gives me an acceptable output (I don't know if it is right) when the values are small. To be clear, I used this:

for ii=1:1:length(X)
XX=X(ii,:);x1=XX(1);x2=XX(2);
gx=[0;cos(x2)];
u=-inv(R)*gx'*XX;
end
 U=[U;u];
 end

Do you think this is a reliable way? The problem is that I expect my u to be saturated because of my if condition in the myode function, but the afterwards computed u is not and I'm not sure if it is correct to have the condition also in my recalculations. I'm using this for a scientific paper so I need to be very accurate and sure about the results I provide.

Thanks a lot

Answer 4

[Raghavendra Ragipani]

Hi,

I had similar problem and I created functions to easily achieve this task. I decided to share it with the MATLAB community. Check it out at the following link.

https://www.mathworks.com/matlabcentral/fileexchange/74184-functions-to-return-additional-variables-from-ode-fun

Call save_var_in_ODE() function inside the ODE program and call get_var_in_ODE() in the main script file where you'll further manipulate it or plot it.

Example usage is given along with the package.

Cheers,

Raghav