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

39 次查看（过去 30 天）
Has MBK
Has MBK 2018-11-11
回答： Raghavendra Ragipani 2020-2-6
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 个评论
Stephen23
Stephen23 2018-11-11
编辑：Stephen23 2018-11-11
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.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Stephen23
Stephen23 2018-11-11
编辑：Stephen23 2018-11-11
"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 个评论
显示 1更早的评论
Stephen23
Stephen23 2018-11-11
编辑：Stephen23 2018-11-11
@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)

请先登录，再进行评论。

更多回答（3 个）

madhan ravi
madhan ravi 2018-11-11
编辑：madhan ravi 2018-11-11
  3 个评论
显示 1更早的评论
Has MBK
Has MBK 2018-11-11
编辑：Has MBK 2018-11-11
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.

请先登录，再进行评论。

Walter Roberson
Walter Roberson 2018-11-11
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 个评论
Has MBK
Has MBK 2018-11-11
编辑：Has MBK 2018-11-11
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

请先登录，再进行评论。

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

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2017b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by