# Why is dsolve solution halved from real answer?

62 views (last 30 days)
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?

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.
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)