Asked by Has MBK
on 11 Nov 2018 at 4:41

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.

Answer by Stephen Cobeldick
on 11 Nov 2018 at 7:09

Edited by Stephen Cobeldick
on 11 Nov 2018 at 10:36

Accepted Answer

"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
on 11 Nov 2018 at 7:55

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.

Stephen Cobeldick
on 11 Nov 2018 at 10:31

@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
on 11 Nov 2018 at 17:18

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

Sign in to comment.

Answer by madhan ravi
on 11 Nov 2018 at 5:22

Edited by madhan ravi
on 11 Nov 2018 at 6:21

Has MBK
on 11 Nov 2018 at 5:48

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
on 11 Nov 2018 at 6:21

madhan ravi
on 11 Nov 2018 at 6:22

See sir Walters answer

Sign in to comment.

Answer by Walter Roberson
on 11 Nov 2018 at 6:21

Has MBK
on 11 Nov 2018 at 6:35

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

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 2 Comments

## madhan ravi (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/429179-how-do-i-extract-an-intermediate-variable-calculated-and-used-inside-my-ode45-function#comment_635553

## Stephen Cobeldick (view profile)

Direct link to this comment:https://www.mathworks.com/matlabcentral/answers/429179-how-do-i-extract-an-intermediate-variable-calculated-and-used-inside-my-ode45-function#comment_635590

Sign in to comment.