Why is dsolve solution halved from real answer?

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?

Answer 1

[David Goodmanson]

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.

[Walter Roberson]

Nice analysis!

[Miles Robertson]

Thanks very much. I realized after posting this that the impulse occuring at the same time as the initial conditions was the root of the problem. I believe the reason the answer was halved from what I expected is because according to Matlab the integral of dirac from negative infinity to zero is 1/2 (not 0 or 1, which makes sense with dirac being the limit of the normal curve as the SD approaches zero), and the initial conditions are being applied there, making an unexpected result.

Since I intended for the initial condition to be at t=0-, and since I only needed this for numerical answers, I shifted the dirac function to the right by a ten-thousandth of a unit, and the numeric answer worked fine and matched the other methods up to 4 decimal places.

However, like you said David, an even better approach would be by hand. I just wanted to get to the bottom of why Matlab wasn't agreeing with my solutions by hand!

Thanks for your time.

[Paul]

@Miles Robertson

I tried shifting the dirac forcing function to the right by 1/10000 and saw the real part of the result accurate to to about 4 decimal places (as you said) with a very small imaginary part.

Edit: the small imaginary part showed up in fplot(imag(h1(t))), where h1(t) was the solution from dsolve. However, simplifying h1(t) showed that its imaginary part is really 0.

An alternative approach is to specify the intial conditions (at least for the derivative) slightly the left of 0, which is a better representation of the problem and yields the same result as the Laplace transform approach:

cond1 = [y(0)==0, dy(-eps)==0];
h1(t) = dsolve(eqn,cond1);
>> h1(t)
 
ans =
 
sin(2*t)/4 + (sin(2*t)*sign(t))/4

h1(0) = 0 and h1(t >0) = sin(2*t)/2, or h1(t) = sin(2*t)/2, for t >= 0.

It gets a bit trickier when the initial conditions are not zero. For example:

cond2 = [y(0)==1, dy(0)==0];
laplace_eqn = subs(laplace(eqn,t,s), laplace(y(t),t,s), Y(s));
laplace_eqn = isolate(subs(laplace_eqn, lhs(cond2), rhs(cond2)), Y(s));
laplace_y2(t) = rhs(ilaplace(laplace_eqn));

cond3 = [y(0)==1, dy(-eps)==0];
h3(t) = simplify(dsolve(eqn,cond3));

Plotting h3(t) differs from laplace_y2 by up to 1e-16. One can do a bit better with:

cond4 = [y(0)==1, dy(-sym(eps)/10)==0];
h4(t) = simplify(dsolve(eqn,cond4));

[Miles Robertson]

Interesting. Thanks for taking a look at that!

[Paul]

@David Goodmanson

I would expect that the Laplace solution has an implied B = 0, i.e,. y(t) = 0 for t<0, but I'm a bit surprised that it doesn't expicitly show that condition in the solution..

Can you say what the additional requirement is that needs to be imposed on dsolve for it to yield the same solution as the others, and how to actually impose that requirement on dsolve? Would that requirement be causality?

[David Goodmanson]

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]

@Miles Robertson

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)