inverse laplace and more and more and more

Hi, I'm trying to do the following: 1.inverse laplace the transfer function, Hs 2.this inverse laplace, ht, is then converted into numeric (from symbolic) 3.and finally convolve this function, htn, with function x1 which is already defined.

>>syms s t1; >>Hs = 100000/(s^5+32.36*s^4+523.59*s^3+5235.924*s^2+32360*s+100000);

>>ht = ilaplace(Hs,t1); %change it to impulse function, and instead of using %default 't', use 't1' to prevent the overlap of variables (t-variable already defined)

>>htn = subs(ht,t); %make it numeric by substituing 't' into 't1' %4. Compute the impulse response of filter, use 'conv' function to %apply the filter to each signal %use 'conv' function: htn*x1,htn*x2,htn*x3

>>YT1=conv(htn,x1);

ISSUE: when I ilaplace 'Hs' I get a really complex result with r3.. (I think that's something to do with quartic equations..right?) Therefore, when I try to change 'ht' to numeric using 'subs' command, the result is another symbolic function, instead of numeric.... and 'conv' command also doesnt work because it only works on numeric functions... T-T

I tried expanding the 'ht' using vpa() command AND 'subs' but that did not work either. What should I do in order to get 'subs' command to work??

Respectfully, Dan

[RahulTandon]

Tried Simulink to solve the problem? Could you please state the final result that yo want?

Answer 1

[Teja Muppirala]

Do you have the Control Systems Toolbox? I think this entire problem would be much easier if you used that instead of symbolic variables and convolutions:

s = tf('s')
Hs = 100000/(s^5+32.36*s^4+523.59*s^3+5235.924*s^2+32360*s+100000)

% Here is a plot of the impulse response of Hs:
impulse(Hs);

% Or if you really wanted the data, and not just a plot:
[ht,tdata] = impulse(Hs);

% Let's apply Hs to a vector x
t = linspace(0,10,1001);
x = sin(t.^2); % A test vector x
figure;
lsim(Hs,x,t); % Apply Hs to x (Same as convolution with ht)

[Daniel]

okay,...

so basically,

instead of using ilaplace, you used "impulse(Hs)"

>>ht=impulse(Hs)

and instead of [conv] command, you used lsim(Hs,x,t)

A quick question. In my original code, I had >>ht=ilaplace(Hs) >>htn=subs(ht,t1)....... Does [impulse()] command do both ilplace and subs?

[Teja Muppirala]

Does impulse(...) do both ilaplace and subs?

Yes, that's one way to think of it. Interally, it's doing it using a different method, but the result is the same.

[Daniel]

kk. thanks alot!

however, if you dont mind, could you explain what lsim(..) exactly does? I simply dont see how lsim(...) works like a conv(...) command.

ALSO!

when using lsim(a,b,c), if the number of variables doesnt match, lsim doesnt work, right? i.e. if a=<449x1 double>,b=<1x15000 double>, and c=<1x15000 double>. How can I remedy this issue?

[Teja Muppirala]

I answered below.

Answer 2

[Teja Muppirala]

The length of the variables only needs to match on the 2nd and 3rd inputs because that is where you are defining the input x(t) to system. And when you define "x(t)" you need to specify a vector for both "x" and for "t", and if they're not the same size, then that wouldn't make any sense right?

How does lsim work? A very good question! Basically it converts the transfer function back into a differential equation and integrates the ODE. You can find out a little bit more in the documentation.

doc lsim

You stumble upon a very enlightening point though. Say I had a differential equation to solve:

y'' + 0.1y' + y = sin(t)

y(0) = 0 and y'(0) = 0

And I want to find out the solution from t = 0 to 10. It turns out there are lots of ways to solve it.

Here I solve it 5 seemingly different ways, but they all give me the same answer to within very minor numerical errors (copy and paste this into the editor and run it). The algorithm in "lsim" basically does something like method 2.

%%1. Solve the ODE symbolically
tspan = 0:.01:10;
y1 = dsolve('D2y + 0.1*Dy + y = sin(t)','Dy(0)=0', 'D2y(0)=0');
y1 = subs(y1,tspan);
plot(tspan,y1);

%%2. Solve the ODE using a numerical ODE solver
hold all;
F = @(t,y) [y(2); -0.1*y(2)-y(1)+sin(t)];
[t,y2] = ode45(F,tspan,[0; 0]);
plot(tspan,y2(:,1));

%%3. Call "lsim" on the transfer function G(s)
s = tf('s');
G = 1/(s^2 + 0.1*s + 1);
u = sin(tspan);
y3 = lsim(G,u,tspan,0);
plot(tspan,y3);

%%4. Get the impulse response and use convolution "conv"
g = impulse(G,tspan);
dt = tspan(2) - tspan(1);
y4 = conv(g,u) * dt;
y4 = y4(1:1001);
plot(tspan,y4);

%%5. Multiply in the Laplace domain, and then take the impulse
U = eval(char(laplace(sym('sin(t)'))));
y5 = impulse(G * U,tspan);
plot(tspan,y5,':');

%%Add a legend
legend({'Method 1','Method 2','Method 3','Method 4','Method 5'})