62 views (last 30 days)

Show older comments

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)

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

Start Hunting!