ilaplace doesn't handle tanh()

While doing convolutions, I wanted to convolve a function f() with a square wave. There are two Laplace transformations of a square wave that I tried, as shown below. The problem occurs when I do the inverse Laplace on the Laplace of the original square wave, or more precisely when I attempt to plot that result. (In the code below I only show the inverse laplace of the laplace transformed squarewave, not the convolution.)

% P1 = ilaplace( 1/s*(1-2/(exp(T/2*s)+1)) );
  P2 = ilaplace( 1/s*tanh(s*T/4) );
% Unfortunately, both don't work with the symbolic toolbox.
assume(P2,'real')
Tsteps = 0.1:0.5:20;
num = vpa(subs(P2,t,Tsteps),6)'; 
stairs(Tsteps(:),num(:,1))

When I inspect the contents of num, the (voluminous!) contents are largely symbolic:

[conj(subs(diff(floor(s), s), s, -400))*(35.334851093590259552001953125 + 0.00000000023283064365386962890625i) + conj(s ... 

My program works with 'simpler' functions, like sinewaves, and then plots OK.

The stairs function faults in the (undocumented) getRealData().

Do I have to assume() something on the arguments of laplace(), ilaplace(), subs() or vpa()?

% getRealData returns the real components of ARGS. If ARGS contains any
% complex data, a warning is displayed. If any of ARGS are not numeric
% and cannot be converted to a double, an error is thrown.

[marcel hendrix]

My apologies. T is a simple constant to scale the period of the square wave, e.g., 10 (seconds) is a good value.

Sorry for not mentioning this in the original question.

Answer 1

[Torsten]

You mean

syms s T
P2 = ilaplace( 1/s*(exp(s*T/4)-exp(-s*T/4))/(exp(s*T/4)+exp(s*T/4)) ,s);
Tsteps = 0.1:0.5:20;
num = vpa(subs(P2,T,Tsteps),6)'

num = 

?

[Paul]

Hi Torsten,

With only two arguments, ilaplace uses the second as the transformation variable. So it should be

syms s T t
%P2 = ilaplace( 1/s*(exp(s*T/4)-exp(-s*T/4))/(exp(s*T/4)+exp(s*T/4)) ,s);
P2 = ilaplace( 1/s*(exp(s*T/4)-exp(-s*T/4))/(exp(s*T/4)+exp(s*T/4)) ,t)

P2 = 

assume(T,'positive')
P2 = simplify(P2,10)

P2 = 

But that's not a square wave. Why is P2 defined that way?

[Torsten]

I forgot a minus sign:

syms s T
P2 = ilaplace( 1/s*(exp(s*T/4)-exp(-s*T/4))/(exp(s*T/4)+exp(-s*T/4)) ,s)

P2 = 

Answer 2

[Paul]

Hi Marcel,

I attempted to recreate the example with the information provided, but got a different result.

syms s
syms T positive
syms t real
P2 = ilaplace( 1/s*tanh(s*T/4),s,t)

P2 = 

% Unfortunately, both don't work with the symbolic toolbox.
assume(P2,'real')
Tsteps = 0.1:0.5:20;

Assume T = 1, for example

num = vpa(subs(subs(P2,T,1),t,Tsteps),6)'; 
num

num = 

stairs(Tsteps(:),num(:,1))

Error using stairs (line 106)
Unable to convert symbolic expression to double array because it contains symbolic function that does not evaluate to number. Input expression must evaluate to number.

I think the fundamental problem is that ilaplace doesn't yield a closed form expression. The only way it could, I believe, is if the symbolic engine actually had a Laplace transform pair defined for a square wave. But the toolbox doesn't have a "square" function like, for example, rectangularPulse, so it's not surprising the toolbox doesn't have a Laplace transform for a function that can't be expressed by the user.

Also, the toolbox doesn't compute the inverse Laplace transform numerically, so I wouldn't expect vpa to work on an unresolved ilaplace expression.

[Paul]

However, we can do this:

syms s
syms T positive
syms t real
S(s) = 1/s*tanh(s*T/4)

S(s) = 

S(s) = rewrite(S(s),'exp')

S(s) = 

s(t) = ilaplace(S(s))

s(t) = 

Square wave with T = 2.

figure
fplot(subs(s(t),T,2),[0 10])
axis padded

[Paul]

We can get to a simpler form of the square wave

syms s
syms T positive
syms t real
F(s) = 1/s*tanh(s*T/4);
F(s) = rewrite(F(s),'exp');
F(s) = simplify(F(s),10);
f(t) = ilaplace(F(s),s,t)

f(t) = 

figure
fplot(subs(f(t),T,2),[0,10])
axis padded

●

Interesting that ilaplace directly yielded f(t), but going the other way doesn't yield F(s)

laplace(f(t),t,s)

ans = 

[marcel hendrix]

This gives (-1)^floor(t)/2 - (-1)^floor(-t)/2 as another symbolic toolbox compatible substitute for the non-existing square wave function. However, in my usage I need the laplace transformed version that you give as F(s), which expands to (exp(s) - 1)/(s*(exp(s) + 1). Unfortunately, this F(s) has the same problems as the previous alternatives: the toolbox does not know how to transform it to something that vpa can handle.

You mentioned rectangularPulse(). I didn't know that one yet. It appears that it is perfectly suited to my work (it supports incredibly weird^H^H^H^H^H flexible options), and allows "easy" dutycycle, frequency, and phase modulation.

[Paul]

Do you mean that you have to use F(s) written in terms of tanh(s) and that you're not allowed to have your program express it in different terms? If that's the case, I don't know what else can be done.

rectangularPulse is only a single pulse, so I don't know how duty cycle, frequency, and phase modulation are applicable.

[marcel hendrix]

Sorry, I don't understand what you say about tanh(s). That was the function I started with and was not understood by the toolbox? Some rewrite of tanh might work with the symbolic toolbox, but I have not seen that in this thread yet.

The function rectangularPulse allows to use functions to define its edges. (I really like how well MATLAB toolboxes are integrated with each other). Below I simply use literals.

%--- Local functions ---%
function res = excitation(t)
res =   (rectangularPulse( 0, 2,t)) + ...
        (rectangularPulse( 4, 6,t)) + ...
        (rectangularPulse( 8,10,t)) + ...
        (rectangularPulse(12,14,t)) + ...
        (rectangularPulse(16,18,t));

P    = ilaplace( M * laplace( E*excitation(t) ) ); % M is an arbitrary laplace function
num  = vpa(subs(P,t,Tsteps),6)';                   % numerical result of P to 6 decimals

For the example, a squarewave source driving a series RLC filter, M is defined as follows:

E = 1; R = 1; L = 1; C=1; T = 10;
A = [ 0  1/C;  -1/L -R/L ];
B = [ 0  R/L ]; 
M = (s*I-A)^-1 * B;

[Paul]

"Some rewrite of tanh might work with the symbolic toolbox, but I have not seen that in this thread yet."

rewrite of F(s) originally given in terms of tanh to instead express F(s) in terms of exp was explicitly shown here and here, not to mention what was shown here.

The example shown above is not a square wave as posed in the original question. It's a different signal, which does not have Laplace transform as originally posed, and is more easily handled by the software.

syms t real
syms s
F(s) = laplace(excitation(t))  % could be expressed more compactly

F(s) = 

However, if the only goal is to get a plot of the output up to a specified time, then computing the output with the square wave (or any other) input cut off at that time is fine (assuming the system under excitation is causal).

Finally, I don't understand the plots. Sending a square wave through that system should not result in staircase like behavior. The output should be continuous.

After defining I and correcting the defintion of B

E = 1; R = 1; L = 1; C=1; T = 10; I = sym(eye(2));
A = [ 0  1/C;  -1/L -R/L ];
%B = [ 0  R/L ];
B = [ 0 ; R/L ];
M = (s*I-A)^-1 * B;
P    = ilaplace( M * laplace( E*excitation(t) ) );

fplot shows a continuous output

figure
subplot(211);
fplot(P(1),[0,20]);grid
subplot(212);
fplot(P(2),[0,20]);grid

●

As does vpa

figure
tval=0:.001:20;
plot(tval,vpa(subs(P(1),t,tval))),grid

function res = excitation(t)
res =   (rectangularPulse( 0, 2,t)) + ...
        (rectangularPulse( 4, 6,t)) + ...
        (rectangularPulse( 8,10,t)) + ...
        (rectangularPulse(12,14,t)) + ...
        (rectangularPulse(16,18,t));
end

[marcel hendrix]

Thank you for the additions!

"rewrite of F(s) originally given in terms of tanh to instead express F(s) in terms of exp was explicitly shown here and here, not to mention what was shown here.

We are talking past each other? Any of the approaches shown give results (a laplace function, e.g. tanh(s)/s) that the symbolic toolbox can't use in symbolic convolution and that can't be plotted.

My intended use is discrete-time (sampled-data) processing for piece-wise linear (power) circuits. None of the details are worked out yet, but the stair-case plot is deliberate.

In a PWL circuit there are complicated consistency (initial condition) problems when going from one "piece" to the next. Laplace and inverse Laplace integration automatically handle this. Standard PWL methods use complementarity, MQCP, and the Katzenelson algorithm. I'd like to see how far I get with Laplace / z-transforms.

[Paul]

"We are talking past each other?"

Apparently. Let's review.

You started this thread by defining the Laplace transform, F(s), of an infinite duration square wave in terms of tanh like so:

syms s
syms t real
syms T positive
F(s) = 1/s*tanh(s*T/4)

F(s) = 

Then we see that applying ilaplace to F(s) fails to return a time domain expression for a square wave

f(t) = ilaplace(F(s),s,t)

f(t) = 

So we are in agreement that leaving F(s) in terms of tanh is not that useful.

But, if we rewrite F(s) by rewriting tanh in terms of exp we get

F(s) = simplify(rewrite(F(s),'exp'),10)

F(s) = 

from which a closed-form expression from the square wave is obtained

f(t) = ilaplace(F(s),s,t)

f(t) = 

At this point, you introduced a new function that is a finite duration sequence of square pulses (that also includes a positive bias), which is a different problem than you asked about from the outset. Such an input may be useful for whatever it is you're trying to do, but it's different.

What we haven't shown yet is if the "rewrite tanh in terms of exp" trick will work when F(s) is multiplied by M(s) to represent convolution of f(t) and m(t) in the time domain. Let's try with the M(s) that you provided

F(s) = 1/s*tanh(s*T/4);
E = 1; R = 1; L = 1; C=1; I = sym(eye(2));
A = [ 0  1/C;  -1/L -R/L ];
%B = [ 0  R/L ];
B = [ 0 ; R/L ];
M(s) = (s*I-A)^-1 * B;
Y(s) = M(s)*F(s);
y(t) = ilaplace(Y(s),s,t)

y(t) = 

Unsurprisingly, ilaplace fails.

Now try the rewrite

Y(s) = rewrite(Y(s),'exp');
y = ilaplace(Y(s),s,t)

y = 

●

If we inspect that expression, we see it contains some strange terms involving the partial derivative of floor(s). I have no idea what that means. I tried a few other approaches to get a usable, closed-form expression, but was not successful. I suppose that this example illustrates your fundamental concern.

At this point, we have two options that I can think of for getting the output of the LTI system in response to the square-wave-like input.

The first would be to window f(t) and proceed from there to get y(t) = ilaplace(M(s)*laplace(f(t)*w(t)),s,t), exactly as you did with your excitation function. The end result of will be valid only up to the time the window cuts off. The disadvantage of this approach is that you need to decide the cutoff time a priori.

The second option would be to develop a symbolic expression for y(t) as a periodic summation. That expression can, in principle, be evaluated out to any finite time. At least I think this option is feasible.

Both options assume that ilaplace(M(s)*P(s)), where p(t) is one period of the input, returns an expression that can be evaluated, which might not be the case for moderately complicated M(s) and/or P(s).