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';
1 个评论
Jan
2016-3-14
Please use the "{} Code" button top format your code. I've done this for you this time.
采纳的回答
Star Strider
2016-3-14
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 ...
2 个评论
Mohamed Aburakhis
2016-3-14
编辑：Mohamed Aburakhis
2016-3-14
Star Strider
2016-3-14
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.
更多回答（1 个）
Jan
2016-3-14
1 个投票
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Function Creation 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
