im using 'dsolve' for 2nd ode with dirac function but matlab returns 1/2 of constant

im trying to make a m.file function in order to apply to underdamped systems

here is the 2nd ode

r" + 0.9214*r' + 9.2141 * r = -4.8202*dirac(t)

the answer should be -1.6066*exp(-0.4607 * t) * sin(3.0003*t)

but matlab returns this

-(24101*900185551^(1/2)*exp(-(4607*t)/10000)*sin((900185551^(1/2)*t)/10000)*sign(t))/900185551

here, 24101*900185551^(1/2)/900185551 is 0.8033 which is 1/2 of 1.6066

what did i miss?? and how can i shorten those long numbers?

here is the code

syms r(t)

dr = diff(r);

ode = diff(r,t,2) + 0.9214*diff(r,t) + 9.2141*r == -4.8202*dirac(t);

cond1 = r(0) ==0;

cond2 = dr(0) ==0;

conds = [cond1 cond2];

ySol(t) = dsolve(ode, conds);

ySol = simplify(ySol)

ySol(t) =

-(24101*900185551^(1/2)*exp(-(4607*t)/10000)*sin((900185551^(1/2)*t)/10000)*sign(t))/900185551

Answer 1

[David Goodmanson]

Hi Song-Ha

Your statement

" the answer should be [constant] *exp(-0.4607 * t) * sin(3.0003*t) "

is missing a very important factor and should be

[constant]*exp(-0.4607*t)*sin(3.0003*t)*sign(t)

Unlike the heaviside function which has a step discontinuity of +1, the sign function has a discontinuity of +2. So you only need half the amplitude that you might think. This code, using your expression for ySol(t)

t = -.3:.0001:.3;
y = -(24101*900185551^(1/2)*exp(-(4607*t)/10000).*sin((900185551^(1/2)*t)/10000).*sign(t))/900185551;
dydt = gradient(y)./gradient(t);
plot(t,y,t,dydt)
grid on

shows that the step discontinuity in the derivative of the function is -4.82 as it should be. Half the drop is from x = (-small) to x = 0, and the other half is from x = 0 to x = (small).

As far as the long numbers, something like vpa(ySol,5) gives you something more sensible to look at, although probably not to calculate with since you are losing precision.

[Song-ha Baek]

i get that 'sign' has a discontinuity of 2, but when I do this below,

>> ySol(1); vpa(ySol(1), 5)

ans = -0.071357

>> -1.6066*exp(-0.4607)*sin(3.0003) % plugged 1 in r(t)

ans = -0.1427

So, if I want to get the second answer, I just have to multiply by 2?

Honestly, I don't get why discontinuity matters and what discontinuity does to the answer.

All I understood about 'sign' is that it returns -1 0 1 depends on the input.

[David Goodmanson]

If you want to get the second answer you would indeed multiply by two. But there is not much point in doing that since it is incorrect and would make the result (which includes the 'sign' factor) too large by that same factor of 2.

If you don't know why the discontinuity matters or what it does, you need to backtrack, go to a textbook or some other source and find out what is going on.

You have

r" + 0.9214*r' + 9.2141 * r = -4.8202*dirac(t)

and if you want to obtain a delta function on the right, there must be a step function somewhere in the definition of r. Then you can produce a delta function by differentiating the step function,. No step function, no delta function. The basic idea is that

d/dx heaviside(x) = delta(x)  
d/dx sign(x) = 2*delta(x)  

and you can multiply both sides of this by an amplitude constant and some other functions to get, in your case, -4-8202*delta(x). If sign is used instead of heaviside, only half the amplitude is required.