Hi, my current code is written for time responses using the lsim function but it needs to be converted to ode45 and i have no idea where to start. Here's my current code with lsim, andy help would be appreciated!

2 views (last 30 days)
% Elevator deflection code
A=[-.045 .036 0 -32.2;-.369 -2.02 176 0;.0019 -.0396 -2.948 0;0 0 1 0];
B=[0;-28.17;-11.88;0];
C=eye(4);
D=[0;0;0;0];
values=eig(A);
T = 0:1:200;
U = ones(size(T));
sys = ss(A,B,C,D);
[Y, Tsim, X] = lsim(sys,U,T);
plot(Tsim,Y)
Y1=(Y(:,1)./176.)/57.7
Y2=Y(:,2)./176
newY= [Y1 Y2 Y(:,3) Y(:,4)]
subplot(2,2,1)
plot(Tsim,Y1)
ylabel('delta(u)/V')
xlabel('Time (s)')
subplot(2,2,2)
plot(Tsim,Y2)
ylabel('delta(alpha)(deg)')
xlabel('Time (s)')
subplot(2,2,3)
plot(Tsim,Y(:,3))
ylabel('q(deg/s)')
xlabel('Time (s)')
subplot(2,2,4)
plot(Tsim,Y(:,4))
ylabel('delta(theta)(deg)')
xlabel('Time (s)')

Accepted Answer

Star Strider
Star Strider on 15 Sep 2014
You have a linear problem, so ode45 is sort of overkill. I use expm for linear problems.
Use the: Y = C*exp(A*t)*D*U representation of a linear control system:
A=[-.045 .036 0 -32.2;-.369 -2.02 176 0;.0019 -.0396 -2.948 0;0 0 1 0];
B=[0;-28.17;-11.88;0];
C=eye(4);
T = 0:1:200;
U = ones(size(T));
for k1 = 1:length(T)
Y(:,k1) = C*expm(A*T(k1))*B*U(k1);
end
figure(1)
plot(T, Y)
producing:
<<www-mathworks-com-matlabcentral-answers-uploaded_files-18154-Hi--20my-20current-20code-20is-20written-20for-20time-20responses-20using-20the-20lsim-20function-20but-20it-20needs-20to-20be-20converted-20to-20ode45-20and-20i-20have-20no-20idea-20where-20to-20start-20--202014-2009-2014.png>>
  2 Comments
Jon Miller
Jon Miller on 15 Sep 2014
Thanks a lot! I'm also pretty familiar with expm function however the class I'm currently taking insists that all of our code use the ODE45 function even if the problem is linear.
Star Strider
Star Strider on 16 Sep 2014
My pleasure!
I apologise — I didn’t realise you were required to use ode45.
In recompense, this should get you started:
A=[-.045 .036 0 -32.2;-.369 -2.02 176 0;.0019 -.0396 -2.948 0;0 0 1 0];
B=[0;-28.17;-11.88;0];
C=eye(4);
T = 0:1:200;
odelin = @(t,x) [-.045 .036 0 -32.2;-.369 -2.02 176 0;.0019 -.0396 -2.948 0;0 0 1 0]*[x(1); x(2); x(3); x(4)] + [0;-28.17;-11.88;0];
[t,y] = ode45(odelin, T, zeros(4,1));
It’s not perfect compared to the expm result (if it is correct as I believe it is). But then I can’t let you miss out on all the fun!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!