MATLAB Answers

Why is dsolve solution halved from real answer?

62 views (last 30 days)
Miles Robertson
Miles Robertson on 14 Feb 2021
Commented: Paul on 18 Feb 2021
I'm trying to figure out why dsolve is being a pill... when I try to solve second order linear differential equations with delta forcing functions, the resulting solution is always half of what it should be.
syms y(t) Y(s) g(t)
d2y = diff(y,t,2);
dy = diff(y,t);
eqn = d2y + 4*y == dirac(t);
cond = [y(0)==0, dy(0)==0];
% dsolve method
dsolve_y(t) = dsolve(eqn,cond);
% laplace transform
laplace_eqn = subs(laplace(eqn,t,s), laplace(y(t),t,s), Y(s));
laplace_eqn = isolate(subs(laplace_eqn, lhs(cond), rhs(cond)), Y(s));
laplace_y(t) = rhs(ilaplace(laplace_eqn));
% green's equation (only solves for t>0)
roots = solve(subs(lhs(eqn),[d2y dy y],[s^2 s 1])==0);
coef = formula(coeffs(lhs(eqn)));
g(t) = (exp(roots(1)*t)-exp(roots(2)*t))/(coef(1)*(roots(1)-roots(2)));
greens_y(t) = simplify(g(t));
% delta change in initial conditions for homogeneous eqn (solves for t>0)
homogeneous_eqn = lhs(eqn) == 0;
new_cond = cond;
new_cond(2) = dy(0)==1/coef(1);
new_cond_y(t) = dsolve(homogeneous_eqn, new_cond);
disp("Functions calculated (last two restricted to t>0)")
disp([dsolve_y(t), laplace_y(t), greens_y(t), new_cond_y(t)])
Why is it that dsolve puts the correct solution at half of what it is supposed to be?

Answers (1)

David Goodmanson
David Goodmanson on 15 Feb 2021
Edited: David Goodmanson on 15 Feb 2021
Hi Miles,
One of your conditions is dy/dt = 0 which is iffy, since there is a dirac delta function at the origin and you are not specifying exactly where that conditions is applied, t = 0+ or 0-.
The general solution to the problem is
y = Acos(2t) + Bsin(2t) t<0
y = Dcos(2t) + Esin(2t) t>0
Interpreting the condition y(0) = 0 to apply to both t=0+ and t=0- forces A=D=0, leaving just the sine terms. If you integrate both sides of the equation
y'' + 4*y = dirac(t)
across the origin from -eps to eps, the result is
y'(0+) - y'(0-) = 1 so
2E-2B = 1
E-B = 1/2
Three of the solutions have E = 1/2 and an implied B = 0, so the function is 0 for x<0. The dsolve solution has E = 1/4, B = -1/4, so it's an even function of t and nonzero for negative t. But the discontinuity in the derivative is correct, so it's a perfectly good solution unless you wish to impose some additional requirement.
  7 Comments
Paul
Paul on 18 Feb 2021
Another option to get the the impulse response from the SMT is to compute the step response and then differentiate the result:
>> g(t) = simplify(dsolve(lhs(eqn) == heaviside(t),cond,'IgnoreAnalyticConstraints',false))
g(t) =
-((sign(t) + 1)*(cos(2*t) - 1))/8
>> h(t) = simplify(diff(g(t),t,1))
h(t) =
(sin(2*t)*(sign(t) + 1))/4
>> simplify(h(t)-sin(2*t)/2*heaviside(t)) % verify against expected result
ans =
0
Those additonal inputs to dsolve are important in this case, as you would see if you run the code above without them, or from this even simpler example:
>> g(t) = dsolve(diff(y(t),t,1)+y(t) == heaviside(t),y(0) == 0) % weird answer
g(t) =
sign(t)/2 - exp(-t)/2 + 1/2
>> g(t) = dsolve(diff(y(t),t,1)+y(t) == heaviside(t),y(0) == 0,'IgnoreAnalyticConstraints',false) % correct answer
g(t) =
exp(-t)*heaviside(t)*(exp(t) - 1)

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!