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

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?

[Jan]

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

Answer 1

[Jan]

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

Answer 2

[prabhjeet singh]

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.

[Jan]

Please explain what "still not able to plot (t,u)" means. I've showed you a method to obtain u. So why ist plut(t,u) not sufficient?

The output will draw exactly what way calculated. If you expect something else, maybe your expectation is not matching the code. But as long as I see the code only, I cannot guess, why you expect something else.

[prabhjeet singh]

Thankyou Jan for showing the method, you have certainly made my day. This problem was bugging me since days. One last thing though, can you please plot t and u and display the image here.

Also why there is a need to use element wise operators as you have suggested before.

Once again, thankyou very much.

[Jan]

When you call fun() with matrices, e.g. x1 will be a vector. Then x1^2 is the multiplication of a [1 x n] with a [1 x n] vector. This is not defined mathematically. Using the elementwise .^ instead squares each element of the vector. Try it:

x = 1:3
x * x   % Failing
x^2     % Failing also
x .* x  % [1, 4, 9]

"can you please plot t and u and display the image here" - so you ask me to apply all needed changes of the code for the elementwise processing by my own to produce the needed diagram? Is this efficient to let me do the work? If you have done this, you can plot u by your own.

[prabhjeet singh]

Thankyou very much again Jan.

Now I have another doubt. If I extract values of 'yr' from function file as per your suggested code, I can easily plot (t,yr). Now instead of using 'yr' from function file and I simpy write formula for 'yr' in the main file and then I plot (t,yr), graph is not same. I get different values of 'yr' (typical sine curve graph which is correct).

I am perplexed why the output 'yr' of function file and 'yr' of main .m file is not same. Hope I have made myself clear.

Thankyou for helping me out. It really means alot :)

Thanks

[Jan]

The yr function insice the function to be integrated is simply:

yr = sin(t)

So please post your code, because I do not undestand, what differs from what.

[prabhjeet singh]

% 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]

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