MATLAB Answers

Has MBK
0

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

Asked by Has MBK
on 11 Nov 2018 at 4:41
Latest activity Commented on by Has MBK
on 11 Nov 2018 at 17:18
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.

  2 Comments

upload your question(equation) as latex form and provide all the necessary details
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.

Sign in to comment.

Tags

Products


Release

R2017b

3 Answers

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);

  3 Comments

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.
@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)
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

  3 Comments

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.
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.

Sign in to comment.


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

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.

  1 Comment

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');
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.