How to get the time vector from the function

Dear all

I have a program the uses two function but I couldn't understand how to get the time to plot the system response. I saved t inside the function but it gives me a different length. So please any one have an idea The code is attached

Many thanks

close all;clear all; clc
global u_save t_save y_save  a_u eta_u segma_r eta_z eta_w a_w vhat;
tf=10;
x1_0=1;
x2_0=3;
uhat_0=2;
a_u=0.232;
eta_u=0.0866;
p11_0=3;
p21_0=0;
p12_0=0;
p22_0=3;
segma_r=0.000001;
eta_z=0.822;
zhat1_0=1;
zhat2_0=1;
eta_w=0.8;
a_w=0.2;
vhat=0;

time=[0:1:tf];       
x0=[x1_0;x2_0;uhat_0;p11_0;p21_0;p12_0;p22_0;zhat1_0;zhat2_0];

[t,x]=ode45('eqp1',time,x0);

figure(1)
plot(t_save,y_save);
xlabel('time (sec)');
ylabel('y');
%legend('y','')
grid on
figure(2)
plot(t_save,u_save);
xlabel('time (sec)');
ylabel('u_hat');
grid on
%%%%%%%%%%%%%%%this is the function
  function [x_dot] = eqp1(t, x)

t        

global u_save t_save y_save  a_u eta_u segma_r eta_z eta_w a_w vhat;

x1=x(1);
x2=x(2);
uhat=x(3);
p11=x(4);
p12=x(5);
p21=x(6);
p22=x(7);
zhat1=x(8);
zhat2=x(9);
u=uhat+a_w*sin(eta_w*t);

y_save=[y_save x2];
u_save=[u_save u];
t_save=[t_save t];

x_dot(1)=-x1+u;
x_dot(2)=-x2+0.5*x1^2;
x_dot(3)=-(a_u*eta_u*zhat2)/(eta_u+a_u*abs(zhat2));
x_dot(4)=eta_z*p11+1*(p12+p21)-eta_z*(p11^2+segma_r*p12*p21);
x_dot(5)=eta_z*p12+1*p22-eta_z*(p11*p12+segma_r*p12*p22);
x_dot(6)=eta_z*p21+1*p22-eta_z*(p11*p21+segma_r*p21*p22);
x_dot(7)=eta_z*p22-eta_z*(p12*p21^2+segma_r*p22^2);
x_dot(8)=(0-eta_z*segma_r*p12)*zhat2+eta_z*p11*(x2-zhat1-vhat);
x_dot(9)=(-eta_z*segma_r*p22)*zhat2+eta_z*p21*(x2-zhat1-vhat);
% x_dot(10)
x_dot=x_dot';

[Jan]

Please use the "{} Code" button top format your code. I've done this for you this time.

Answer 1

[Star Strider]

I can’t understand what you want to do. It’s difficult to follow your code.

Also, I would not use global variables. They cause more problems than they solve. They used to be necessary to pass extra parameters to ODE functions, but that has not been true for at least 15 years (at least as I remember). You don’t say what version of MATLAB you’re using.

Looking through your code, if you want to plot ‘x(2)’ as a function of time (since that seems to be what ‘y_save’ is), I would just do:

... PREVIOUS CODE ...
[t,x]=ode45(@eqp1,time,x0);

figure(1)
plot(t, x(:,2));
... REMAINING CODE ...

[Mohamed Aburakhis]

thanks very much, your command works but the length of "t" is not right, how can i make it as the tf of my codes? I use matlab 15a

[Star Strider]

That value of ‘t’ is created by ode45. You gave ode45 a vector of times in your ‘time’ array, not a range, so if you plot all your variables as functions of the same ‘time’ vector, everything should work correctly. (The ‘t’ vector returned by ode45 will be identical to your ‘time’ vector.)

I cannot follow what you are doing in your code, so beyond that and what I wrote previously, I cannot help. It is extremely confusing with all the global statements — none of which are necessary. If you want to plot more of the ‘x’ array variables that ode45 solves, address them as I addressed ‘x(:,2)’ in the plot call, and plot them the same way.

Also, delete all these from your ‘eqp1’ ODE function:

y_save=[y_save x2]; 
u_save=[u_save u]; 
t_save=[t_save t];

If you want to create your ‘u’ array in your script workspace, you can do that simply by calculating it there after your ode45 call:

u = x(:,3) + a_w*sin(eta_w*t);

There should be no problems because both ‘x(:,3)’ and ‘t’ are column vectors, so ‘u’ will also be a column vector.

Answer 2

[Jan]

Storing the time inside the function to be integrated is not meaningful: The integrator can reject steps and calls the function several time to compute one step. So rely on the output "t" from ODE45 only.