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!

% 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)')

Answer 1

[Star Strider]

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

[Jon Miller]

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]

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!