Third Order ODE with unit step input

I have been trying to solve this differential equation for two days now. I do not know what to do with the right hand side of the ODE. The only way I have seen to solve it does not include the derivative of the input as well. Would really appreaciate some help atleast to know how to start it up.

y^''' (t)+6y^'' (t)+11y^'(t) +6y(t)=u^'' (t)+2u^' (t)+3u(t)

y’’(0) = 1 ; y’(0) = -1; y(0) = 1

where u=Unit step Us(t).

Ive tried to do it in simulink but the answers there havent been coming out right either.

[Paul]

Have you tried using the Laplace transform technique?

[Parveen Ayoubi]

I have tried it but I dont know how to do it without giving it the entire transfer function. The goal is to not solve anything by hand before inputting the ode into matlab. Matlab should be able to solve the entire ode without any analytical help from the user.

[Paul]

Solve the ode numerically or exactly? Can you post some code of what you've already tried and show where the sticking point is?

The third example on this doc page may be of interest.

[Parveen Ayoubi]

I checked out that doc page and I will try some of those things, although I was using ode45 and dsolve both to try to figure it out. Ill write out the last code I tried and send the frequency domain one.

-One note, I am a beginner at matlab

syms y(t);
syms u(t);
Dy = diff(y,t);
D2y = diff(y,t,2);
u = heaviside(t); // I saw online that this is unit step
Du = diff(u,t);
ode = diff(y,t,3) + 6*diff(y,t,2)+11*diff(y,t)+6*y == diff(u,t,2)+2*diff(u,t)+3*u;
cond1 = Dy(0) == -1;
cond2 = D2y(0) == 1;
cond3 = y(0) == 1;
conds = [cond1 cond2 cond3];
ySolved(t) = dsolve(ode,conds); 
uSolved(t) = dsolve (ode);

I dont really know where to start with it so ive been trying all kinds of different codes. I dont know how to add the input to the problem properly. I do know what the final answer is because I solved it analytically and then applied the transfer function and laplace inverse on matlab.

X = ilaplace((s^2+5*s+6+3/s)/(s^3+6*s^2+11*s+6));

=> 3*exp(-2*t))/2 - exp(-t)/2 - exp(-3*t)/2 + 1/2

This would be the transfer function for the given u=Us(t) in the ode.

Is there a way for me to give the full ode (y^''' (t)+6y^'' (t)+11y^'(t) +6y(t)=u^'' (t)+2u^' (t)+3u(t)) to matlab and it give me the laplace of it. so that I can find the transfer function and ilaplace to find the answer without me having to solve it analytically.

[Paul]

Can you show how you arrived at the equation X = ilaplace(....) ?

[Walter Roberson]

X = ilaplace((s^2+5*s+6+3/s)/(s^3+6*s^2+11*s+6));

Perhaps you meant

X = ilaplace(((s^2+5*s+6+3)/s)/(s^3+6*s^2+11*s+6));

[Parveen Ayoubi]

I meant it as X = ilaplace((s^2+5*s+6+3/s)/(s^3+6*s^2+11*s+6)); correctly, I got this by solving it by hand. I was getting what Y(s)/U(s) is equal to, the way I did this was by dividing the whole equation by U(s)=1/s which ended up giving me that answer.

[Paul]

Try as I might, I can't recreate your result. Here's what I get

syms y(t) u(t)
ode = diff(y,t,3) + 6*diff(y,t,2) + 11*diff(y,t) + 6*y(t) == diff(u,t,2) + 2*diff(u,t) + 3*u(t);
Leqn = laplace(ode);
syms Y U s y0 Dy0 D2y0
Leqn = subs(Leqn,[laplace(y(t),t,s) subs(diff(u(t), t), t, 0) u(0) laplace(u,t,s)],[Y 0 0 U]);
Leqn = subs(Leqn,[laplace(y(t),t,s) y(0) subs(diff(y(t), t), t, 0) subs(diff(y(t), t, t), t, 0) subs(diff(u(t), t), t, 0) u(0) laplace(u,t,s)],[Y y0 Dy0 D2y0 0 0 U])

Leqn = 

Y = solve(Leqn,Y);
Y = subs(Y,U,1/s);
Y = subs(Y,[y0 Dy0 D2y0],[1 -1 1])

Y = 

[Parveen Ayoubi]

Wouldnt you be missing the initial conditions of U when changing from time domain to frequency? I think that is the difference, Because if we have unit step for the input then u(0)=1 u'(0) is equal to one as well. If you can show me how i would add those conditions to the code you did it would be great. Or atleast a hint

[Parveen Ayoubi]

Nevermind im incorrect, which makes all my work wrong and im at another wall. Because if im assuming that t>0 like walter robinson had in his code it wouldnt be the correct answer

Answer 1

[Walter Roberson]

sympref('HeavisideAtOrigin', 0);
syms y(t)
assume(t > 0)
Dy = diff(y,t);
D2y = diff(y,t,2);
u(t) = heaviside(t); %// I saw online that this is unit step
Du = diff(u,t);
ode = diff(y,t,3) + 6*diff(y,t,2)+11*diff(y,t)+6*y == diff(u,t,2)+2*diff(u,t)+3*u;
cond1 = Dy(0) == -1;
cond2 = D2y(0) == 1;
cond3 = y(0) == 1;
conds = [cond1 cond2 cond3];
sol = dsolve(ode,conds)

sol = 

The assumption that t > 0 is not strictly valid, but if you use t >= 0 you get sign(t) and heaviside() in the equations, and you would have to break it into two cases. You could continue on with

assume(t == 0)
sol0 = dsolve(ode,conds)

sol0 = 

... which happens to give the same result.

[Walter Roberson]

sympref('HeavisideAtOrigin', 0);
syms y(t)
assume(t,'real');

Dy = diff(y,t);
D2y = diff(y,t,2);
u(t) = heaviside(t); %// I saw online that this is unit step
Du = diff(u,t);
ode = diff(y,t,3) + 6*diff(y,t,2)+11*diff(y,t)+6*y == diff(u,t,2)+2*diff(u,t)+3*u;
cond1 = Dy(0) == -1;
cond2 = D2y(0) == 1;
cond3 = y(0) == 1;
conds = [cond1 cond2 cond3];

sol = simplify(dsolve(ode,conds))

sol = 

sol0 = limit(sol, t, 0)

sol0 = 

1

syms tpos tneg
assume(tpos > 0)
assume(tneg < 0)
solpos = subs(simplify(subs(sol, t, tpos)),tpos,t)

solpos = 

solneg = subs(simplify(subs(sol, t, tneg)),tneg,t)

solneg = 

[Walter Roberson]

However, if you take the laplace transform approach, then those are typically only valid for non-negative time.