Why is dsolve solution halved from real answer?

7 views (last 30 days)
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
David Goodmanson
David Goodmanson on 17 Feb 2021
Hi Paul,
Since the usual Laplace transform assumes y=0 for t<0, maybe they don't see the need to state that result explicitly.
Causality can be seen as the condition that imposes B = 0. I am not that familiar with dsolve, and not much motivated to try and figure out what scheme of inputs might get it done. I suppose that's being lazy, but it's pretty hard to get a general function such as dsolve to behave itself for all possible limiting-case type boundary conditions. And your own judgment is going to be the final arbiter anyway, so I am not sure how much effort is justified in trying to coddle dsolve into giving a particular result.
For causality it helps to look at the fourier transform. And since the inverse laplace transform is basically the fourier transform turned by 90 degrees, they are not all that different.
The diffyq
y'' + w0^2 y = delta(t)
has the solution
y(t) = -1/(2pi) Integral{-inf,inf} 1/((w-w0)(w+w0)) exp(iwt) dw
and has two poles on the axis of integration. Causality requires the path to go just under each pole, so that for t<0 the path of integration can be closed around the lower half plane, enclosing no poles, resulting in y=0. For t>0 the path is closed in the upper half plane, enclosing both poles and resulting in
y = sin(w0t)/w0 t>0.
If the path goes over either or both poles, there is nonzero y(t) for t<0.
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.

Categories

Find more on Mathematics in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!