Is the Result from dsolve() Justified for a heaviside() Input?

2 views (last 30 days)
This question is a follow up to this one that I feel warrants a separate discussion. The problem is ...
Solve this differential equation
syms y(t) u(t)
ode = diff(y(t),t,3) + 6*diff(y(t),t,2) + 11*diff(y(t),t) + 6*y(t) == diff(u(t),t,2) + 2*diff(u(t),t) + 3*u(t)
ode = 
with u(t) = heaviside(t) and initial condtions zero for t = -1. Using dsolve we get
ydsol(t) = simplify(dsolve(subs(ode,u(t),heaviside(t)),[y(-1)==0, subs(diff(y(t),t),t,-1)==0 subs(diff(y(t),t,2),t,-1)==0]))
ydsol(t) = 
If ydsol(t) is the solution, it should satisfy the IC's
subs([ydsol(t), diff(ydsol(t),t) diff(ydsol(t),t,2)],t,-1)
ans = 
which it does. But it should also satisfy the differential equation
simplify(subs(lhs(ode),y(t),ydsol(t))) == simplify(subs(rhs(ode),u(t),heaviside(t)))
ans = 
which it does not.
The result ydsol(t) doesn't look correct. Is this user error on my part?
The correct answer can be obtained by first finding the impulse response
yimp(t) = simplify(dsolve(subs(ode,u(t),dirac(t)),[y(-1)==0, subs(diff(y(t),t),t,-1)==0 subs(diff(y(t),t,2),t,-1)==0]))
yimp(t) = 
and then integrating the impulse response to get the step response
syms tau
ystep(t) = simplify(int(yimp(tau),tau,0,t)*heaviside(t))
ystep(t) = 
Verify the initial conditions
subs([ystep(t), diff(ystep(t),t) diff(ystep(t),t,2)],t,-1)
ans = 
And verify the solution satisfies the differential equation
simplify(subs(lhs(ode),y(t),ystep(t))) == simplify(subs(rhs(ode),u(t),heaviside(t)))
ans = 
ystep(t) is clearly (IMO) correct.
We can also obtain the solution via the Laplace transform, pretty much by inspection
syms s
ylap(t) = ilaplace((s^2 + 2*s + 3)/(s^3 + 6*s^2 + 11*s + 6)/s)*heaviside(t)
ylap(t) = 
which is identical to ystep(t).
The appoximate solution from the Control System Toolbox is
tsim = 0:1e-3:10;
ycst = step(tf([1 2 3],[1 6 11 6]),tsim);
FWIW, here's how the solutions compare graphically
fplot([ydsol(t) ystep(t)],[-1 10])
hold on
Can anyone explain how the result ydsol(t) might be justified?

Answers (0)




Community Treasure Hunt

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

Start Hunting!