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).

[marcel hendrix]

Your post covers my top-level problem pretty well. My concern was that I wasn't using the toolbox right, and that that was the reason of the failure to plot. The answer that I will go with is that the toolbox simply does not (yet) have the 'function vocabulary' needed for my problem. I have since that time found numerical ilaplace functions on MATLAB Central that fit my problem very well with quite simple code and no need for the symbolic toolbox.

Sorry that I got carried away and obfuscated the problem by mentioning and using the rectangularPulse function with only five terms (because I needed to plot only between 0 and 5*T). Indeed using only 5 terms makes the infinite sum different from a function with only 5 terms, and some 'endeffects' could be expected. However, it seems to work without such artifacts.

[Paul]

Here's an example of the second option I mentioned in this comment regarding a closed form solution expressed as a perodic summation, demonstrated with the example problem of @David Goodmanson so we can compare solutions.

Laplace transform of the current across the resistor

syms R L C U(s) V(s)
eqn = (s*L + R + 1/(s*C))* U(s) == V(s)

eqn = 

eqn = isolate(eqn,U)

eqn = 

Laplace transform of the voltage across the resistor

syms V_r(s)
eqn = V_r(s) == R*rhs(eqn)

eqn = 

[num,den]=numden(rhs(eqn));
eqn = lhs(eqn) == num/den

eqn = 

Define the input voltage, square wave of amplitude V_sq, period T = 2*a

syms a positive 
syms V_sq real
eqn = lhs(eqn) == subs(rhs(eqn),V,V_sq*tanh(a*s/2)/s)

eqn = 

Sub in the values of the constants from David's example

eqn = subs(eqn,[V_sq,L,R,C,a],[10,1,20,6e-5,1/2])

eqn = 

V_r is the of the form tanh(b*s)*Z(s)

Stating w/o proof that v_r(t) is then

syms t b real
syms z(t) v_r(t)
syms n integer
try
    v_r(t) = z(t) - 2*symsum((-1)^n*heaviside(t - 2*b - 2*b*n)*z(t - 2*b - 2*b*n), n, 0, inf)
catch ME
    ME.message
end

ans = 'Recursive definition: Reached maximal depth for nested procedure calls.'

Not sure why there is an error here, I was expecting to see an unevaluated symbolic expression.

Proceeding ...

b = 1/sym(4);

Inverse transform of Z(s)

z(t) = ilaplace(rhs(eqn)/tanh(s/4),s,t)

z(t) = 

We only need a finite upper bound of the periodic summation to represent v_r(t) for a finite time. For example, choosing an upper bound of n = 20

v_r(t) = z(t) - 2*symsum((-1)^n*heaviside(t - 2*b - 2*b*n)*z(t - 2*b - 2*b*n), n, 0, 20);

And evaluating to t = 4

figure
fplot(v_r(t),[0,4]);
ylim([-3,3]),grid

Answer 3

[David Goodmanson]

Hi marcel,

There seems to be some uncertainty about the question, but this answer is for the voltage across the resistor, for a series RLC circuit driven by a square wave voltage about zero, amplitude V. In what follows the square wave period is 2a rather than T, so each half of the wave has length a, which I think is a bit easier for derivation. To get back to period T, it's easy enough to replace 'a' by T/2 in the result. Here are a couple of examples from the code below:

T = 1;  V = 10;  L = 1;  R = 80;  C = .2

Here the values are such that a not speedy rise time is evident and there is some droop across the top of the wave. If you were to use, say L = 1 R = 1000 C = 1e6 you get basically a perfect square wave since the resistor has all the impdance.

For

T = 1;  V = 10;  L = 1;  R = 20;  C = 6e-5

you get an underdamped circuit with a lot of ringdown.

Details:

As has been noted, ilaplace does not work for the function tanh(sT/4)/s = tanh(sa/2)/s, but when tanh is written out, you get

syms s t a b
f =  (exp(a*s/2)-exp(-a*s/2))/(exp(a*s/2)+exp(-a*s/2));     % = tanh(a*s/2)
g0 = ilaplace(f/s,s,t)
     g0 = (-1)^floor(t/a)/2 - (-1)^floor(-t/a)/2

which works. Using

     floor(x) + floor(-x) = -1    % x not an exact integer

the two terms are identical, so you get

     g0 = (-1)^floor(t/a)     % square wave, period 2a

Another ilaplace that will be needed is the same thing with the denominator = (s+b) instead of s:

gb = ilaplace(f/(s+b),s,t)
          gb  = (exp(-b*t)*((-exp(a*b))^floor(t/a) - 1))/(exp(-a*b) + 1) ...
              - (exp(-b*t)*((-exp(-a*b))^floor(-t/a) - 1))/(exp(a*b) + 1)

I was suprised that syms was able to come up with this, but it did. In this case the two terms are not identical but you can define

          n = floor(t/a)

and combine them to get

gb = ((1-exp(a*b))/(1+exp(a*b)))*exp(-b*t) + (2/(1+exp(a*b)))*(-1)^n*exp(a*b*(n+1))*exp(-b*t)

(This can be derived using a delta function pulse train).

The second term has a problem that exp(a*b*(n+1)) is rapidly increasing with t so the calc will bomb out for large enough t before the exp(-b*t) factor can bring it back down. The way out is to define

t = a*floor(t/a)+tbar = a*n+tbar  with  0<= tbar <a   (tbar is local to each interval)

and then

gb = ((1-exp(a*b))/(1+exp(a*b)))*exp(-b*t) + (2*exp(a*b)/(1+exp(a*b)))*(-1)^n*exp(-b*tbar)    (1)

For the circuit aspect, (current = u since cap I shows up badly in this font)

(sL + R + 1/(sC)) u(s) = V(s) =  Vf(s)/s         % f = tanh(a*s/2)

and

Vres(s)  = (VR/L) f(s) / (s^2 + (R/L) s + 1/(LC) )

Now factor the polynomial into its roots, (change their sign since we need s+b1 not s-b1 etc.) and use partial fractions

1/(s+b1)(s+b2) = (1/(s+b1) -1/(s+b2)) / (b2-b1)
so
Vres(s) = (VR/L)/(b2-b1) (f(s)/(s+b1) - f(s)/(s+b2))

the ilaplace of this is known, so going back to the time domain, the expression (1) for gb will be evaluated for b=b1 and b=b2:

V = 10;
L = 1;
% R = 80;  C = .2;
R = 20;   C = 6e-5;
T = 1;
a = T/2;

t = 0:.001:4;
n = floor(t/a);
tbar = t-n*a;
gb = @(a,b,t,n,bar) ((1-exp(a*b))/(1+exp(a*b)))*exp(-b*t) ...
    + (2*exp(a*b)/(1+exp(a*b)))*(-1).^n.*exp(-b*tbar);
b =-roots([L R 1/C])
g = (V*R/(L*diff(b)))*(gb(a,b(1),t,n,tbar)-gb(a,b(2),t,n,tbar));

figure(1)
plot(t,g)
grid on

To this solution can be added any linear combination of the homogeneous solutions exp(-b1*t) and exp(-b2*t). This will change the boundary conditions at t=0. However, the solution as is gives the right boundary conditions for a quiescent circuit before the square wave is applied.

[marcel hendrix]

Thank you for this thoughtful post. I will remark that using laplace and ilaplace in pairs for solving ODE's automatically handles initial conditions, inconsistent state variables, and the need to laboriously piece zero-state and zero-input solutions together. Those features alone make it very useful for my work. Being able to do this symbolically is in most cases only a party-trick. It is nice to impress colleagues and infrequently it adds unique features, but in general it is not very practical to carry along page-long equations and numerical inversion is fine.

[David Goodmanson]

Your question seemed oriented toward symbolics so I included some, but in general my opinion of symbolics is similar to yours. In some cases 'party trick' is apt, because symbolics will spew out some crazy long expression that's almost impossible to interpret. Sometimes trying to wrangle symbolics into doing what you want just isn't worth it. Giving credit where it's due, in this case symbolics was a useful check, but I had done a derivation already using the lemma

invL(f(s)) = f(t)
then
invL(f(s)/(s+b)) = e^(-bt) Int{0,t} f(t') e^(bt') dt'

on a delta function pulse train. Thanks for mentioning zero state and zero input. Are there good numerical inverse-laplace programs/code that you could recommend?

[marcel hendrix]

The MATLAB file-exchange has 3 implementations. I think that "A CME-based numerical inverse Laplace transformation method" is the best one (based on running the examples) but difficult to read because of its optimizations. The invlab code is straightforward and gives good results. I would avoid the gavsteh procedure based on the examples and my simple initial tests.

[Paul]

If only a numerical solution to the ODE is needed, then why bother using the Laplace transform in the first place?

[marcel hendrix]

The features of symbolic/numerical laplace and ilaplace mentioned here, benefit numerical code that uses them. This is a significant simplification in programming that code and I intend to make use of that as mentioned here.

[David Goodmanson]

Hi marcel, the links sound interesting but appear to be self-referential.

[marcel hendrix]

Here is the IEEE Xplore reference for a good paper if you have an engineering background: Ajoy Opal, 'The Transition Matrix for Linear Circuits', IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, Vol. 16, No 5, May 1997. But maybe The_Matrix_Exponential_via_the_Laplace_Transform is enough?