How to extract variable from a function file while using ode45?

Question

prabhjeet singh 2021-3-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/778452-how-to-extract-variable-from-a-function-file-while-using-ode45

评论： Jan 2021-3-22

x = [0.8 -0.2 1 0 0 0.8 -0.2]';
tspan = [0 20];
[t,x]= ode45(@(t,x) fun(t,x),tspan,x);
function [dot] = fun(t,x)
k1 = 3;
k2 = 3;
k3 = 3;
x1 = x(1);
x2 = x(2);
x3 = x(3);
z2 = x(4);
z3 = x(5);
S1 = x(6);
S2 = x(7);
yr = sin(t);
s1 = x1 - yr;
s2 = x2 - z2;
s3 = x3 - z3;
tau2 = exp(-t) + 0.05;
tau3 = exp(-t) + 0.05;
alpha2 = -k1*s1 - (x1^2 + x2^3 +x3 -x2) + cos(t);
z2dot = (alpha2 - z2)/tau2;
alpha3 = -k2*s2 - (x1^2*x2 + x3^5 - x3) + z2dot;
z3dot = (alpha3 - z3)/tau3;
u = -k3*s3 - x1*x2*x3^2 + z3dot; 
x1dot = x1^2 + x2^3 + x3;
x2dot = x1^2*x2 + x3^5;
x3dot = u + x1*x2*x3^2;
S1dot = x1dot - cos(t);
S2dot = x2dot - z2dot;
dot = [x1dot x2dot x3dot z2dot z3dot S1dot S2dot]';

For example, I want to extract variable u from this code, is it possible to do so?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan 2021-3-20

Please use the tools for formatting code in the forum. I've done this for you this time.

请先登录，再进行评论。

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

Answer 1

Jan 2021-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/778452-how-to-extract-variable-from-a-function-file-while-using-ode45#answer_653342

编辑：Jan 2021-3-20

在 MATLAB Online 中打开

This question is asked several times per day. I ask MathWorks to include an explanation in the documentation.

Yes, of course it is possible. Simply modify the function to be integrated such, that it accepts vectors as input and use the output of ODE45 as input. There is a tricky need to transpose the inputs:

[t, x] = ode45(@fun, tspan, x);  % @fun is faster than @(t,x) fun(t,x)
[~, u] = fun(t.', x.');
function [dot, u] = fun(t, x)
k1 = 3;
k2 = 3;
k3 = 3;
x1 = x(1, :);
x2 = x(2, :);
x3 = x(3, :);
... etc.
% Use elementwise operators: * ==> .* ,  ^ ==> .^ ,  / ==> ./
end

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

prabhjeet singh 2021-3-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/778452-how-to-extract-variable-from-a-function-file-while-using-ode45#answer_653462

Thankyou very much for the instant respone.

I am still not able to plot (t,u). I am expecting response similar to the attachment. Can you please help regarding this. Thanks a ton in advance.

7 个评论
显示 5更早的评论隐藏 5更早的评论

prabhjeet singh 2021-3-22

在 MATLAB Online 中打开

% Main file
x = [0.8 -0.2 1 0 0 0.8 -0.2]';
tspan = [0 20];
[t,x]= ode45(@fun,tspan,x);
[~ , u] = fun(t.',x.');
yr = sin(t);
figure(1)
plot(t,yr,'--r','Linewidth',1.5);
% Function file
function [dot,u] = fun(t,x)
k1 = 3;
k2 = 3;
k3 = 3;
x1 = x(1);
x2 = x(2);
x3 = x(3);
z2 = x(4);
z3 = x(5);
S1 = x(6);
S2 = x(7);
yr = sin(t);
s1 = x1 - yr;
s2 = x2 - z2;
s3 = x3 - z3;
tau2 = exp(-t) + 0.05;
tau3 = exp(-t) + 0.05;
alpha2 = -k1.*s1 - (x1.^2 + x2.^3 + x3 - x2) + cos(t);
z2dot = (alpha2 - z2)./tau2;
alpha3 = -k2.*s2 - (x1.^2.*x2 + x3.^5 - x3) + z2dot;
z3dot = (alpha3 - z3)./tau3;
u = -k3.*s3 - x1.*x2.*x3.^2 + z3dot; 
x1dot = x1.^2 + x2.^3 + x3;
x2dot = x1.^2.*x2 + x3.^5;
x3dot = u + x1.*x2.*x3.^2;
S1dot = x1dot - cos(t);
S2dot = x2dot - z2dot;
dot = [x1dot x2dot x3dot z2dot z3dot S1dot S2dot]';
end

Now if you execute this, I get plot of yr and t. But if I edit code as per your suggestion, I get different plot. Hope it is clear now.

Sorry for bothering you alot for this question.

Jan 2021-3-22

在 MATLAB Online 中打开

"if I edit code as per your suggestion" - this still does not allow to understand, what you are doing to obtain different values for yr. You just show one of the methods, but what is the other?

[~ , u, yr] = fun(t.', x.');
figure;
plot(t, yr);
hold on
plot(t, sin(t), 'ro');
function [dot,u,yr] = fun(t,x)
...
yr = sin(t);
...
end

Except for transposition of the vector both outputs are equal.

You questions about Matlab are welcome in this forum. You do not "bother", but try to solve a problem. This is the purpose of this forum and if I do not find the time to assist you, somebody else will do this.

请先登录，再进行评论。