4 views (last 30 days)

Show older comments

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

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.

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.

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

Start Hunting!