iLaplace of simple functions that include exponential and polynomial terms

Here is the laplace transfer function of a square waveform applied to a very simple LC filter.

syms s t
time=0:4e-9:8e-7;
signal=-(1.5e32*(exp(-1e-7*s)-1.0)*(5e31*s^2+5e38*s+5.0e46))/(s*(exp(-1e-7*s)+1.0)*(2.5e63*s^2+2.5e70*s+2.5e78));
figure
plot(time,subs(vpa(ilaplace(signal)),t,time))

Why this simple code cannot be calculated by Matlab?

Answer 1

[Walter Roberson]

Your floating point constants lead MATLAB to a situation in which it is internally generating derivative of the floor function applied at numeric value. MATLAB does not have derivative of floor built in.

If you rewrite so that you get pure rationals, then MATLAB is able to convert the derivative of floor into a heaviside.

syms s t
Q = @(v) str2sym(v);
time=0:Q('4*10^-9'):Q('8*10^-7');
signal=-(Q('15*10^31')*(exp(Q('-1*10^-7')*s)-Q('1'))*(Q('5*10^31')*s^2+Q('5*10^38')*s+Q('5*10^46')))/(s*(exp(Q('-1*10^-7')*s)+Q('1'))*(Q('25*10^62')*s^2+Q('25*10^69')*s+Q('25*10^77')));
iL = simplify(ilaplace(signal));
siL = subs(iL, t, time);
niL = double(siL);
plot(time, niL)

Some of the above lines can be combined. I separated operations out to make it easier to track where the problem was coming from.

Notice that sym('8e-7') creates a software floating point number, which will not work for this purpose; it is necessary to build the equivalent rational expression. That is also why 1.5e32 got converted to 15e31 and then to 15*10^31 -- using even a single decimal point in the coefficients can result in something that ilaplace cannot properly convert to heaviside.

[Star Strider]

Interesting. I kept getting a ‘D’ function and couldn’t figure out what it was or what it meant.

+1 for that elucidation!

[S H]

Great suggestion. Thank you Walter. Is there a function that automates these number conversions for automation?

[Walter Roberson]

There is an obscure way. Even knowing it exists, I always have trouble finding the documentation for it.

And I ended up having to combine it with other obscure ways... It might have been easier to rewrite it all in MuPAD.

https://www.mathworks.com/help/symbolic/sym.mapsymtype.html

rat_S = mapSymType(signal,'vpareal',@(s) feval(symengine, 'coerce', s, 'DOM_RAT'));

However, it turns out that even though rat_S will contain only rationals, that ilaplace of it will contain D(floor)(1 - 10000000*t) in at least four places, so the limitation is not really rationals: the limitation has something to do with the precision of the values.

As a general rule: mixing floating point with formula transformation or manipulation or solving, often leads to grief. Work with rationals or algebraic numbers such as sqrt(sym(2)), or sym('pi'), instead of with floating point. Pretty much leave symbolic floating point for vpaintegral() and vpasolve() to deal with and avoid it otherwise.