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

4 views (last 30 days)
Song-ha Baek on 1 Jun 2020
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

David Goodmanson on 1 Jun 2020
Edited: David Goodmanson on 1 Jun 2020
Hi Song-Ha
" 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;
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.
##### 2 CommentsShowHide 1 older comment
David Goodmanson on 2 Jun 2020
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.