dsolve and heaviside (unit step response)

9 views (last 30 days)
Hi, I'm trying to understand the behaviour of dsolve and heaviside for solving a simple ODE.
syms y(t)
Y = dsolve(diff(y)+y == 1, y(0) == 0)
ezplot(Y,[0,5])
Gives the desired response of a first order LTI ODE to a step, i.e. Y = 1 - exp(-t) (not too worried about negative time). However,
sympref('HeavisideAtOrigin', 1)
Y = dsolve(diff(y)+y == heaviside(t), y(0) == 0)
ezplot(Y,[0,5])
Produces Y = sign(t)/2 - exp(-t)/2 + 1/2 so the value at Y(0) is 0.5, rather than starting from Y(0)=0.
Given the same initial value and forcing function 1 for t>=0, they should produce the same result? Obviously, the first result is the one I was expecting.

Accepted Answer

Martin Brown
Martin Brown on 27 Feb 2018
Edited: Martin Brown on 27 Feb 2018
Got in touch with Mathworks "bugs" and this is only correct when you include the 'IgnoreAnalyticConstraints', false flag. Without this, the step response is wrong for 1st order, 2nd order .... Apparently "By default , the solver applies some simplifications while solving differential equations which may lead to unexpected results in rare cases". Not sure how rare a step response is though.
syms y(t) Dy = diff(y); D2y = diff(y,2);
% 1st order step Y1 = dsolve(Dy+y==heaviside(t), y(-1)==0); Y2 = dsolve(Dy+y==heaviside(t), y(-1)==0, 'IgnoreAnalyticConstraints', false); figure(1); ezplot(Y2-Y1);
% 2nd order step Y1 = dsolve(D2y+Dy+y==heaviside(t), y(-1)==0, Dy(-1)==0); Y2 = dsolve(D2y+Dy+y==heaviside(t), y(-1)==0, Dy(-1)==0, 'IgnoreAnalyticConstraints', false); figure(2); ezplot(Y2-Y1);

More Answers (1)

Torsten
Torsten on 21 Feb 2018
Why Y(0)=0.5 ?
Y(0) = sign(0)/2-exp(-0)/2+1/2 = 0-1/2+1/2 = 0.
Best wishes
Torsten.
  1 Comment
Martin Brown
Martin Brown on 21 Feb 2018
Edited: Martin Brown on 21 Feb 2018
Sorry, I was being a bit lax with the question, I should have put Y(0+) = 0.5 You're correct at the point Y(0) = 0. However, this jump isn't really what would be expected. I've tried it on an old Matlab version at work (2014) and the heaviside / dsolve combination gives the "correct" answer Y = exp(-t)*heaviside(t)*(exp(t) - 1) where Y(0) = Y(0+) = 0 I don't understand why the later versions 2016 and 2017 give different (wrong?) answers, as described above.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!