how to solve laplace transform problem in matlab

[David Goodmanson]

Why didn't anyone ever write down the correct answer? If someone wants to use the laplace transform, it helps to at least know what to aim for.

y = (1/3) - (1/2) e^(-t) + (1/6) e^(-3t)

[Paul]

Hi David,

When I seen an initial value problem like this with direction to solve with the Laplace transoform, I automatically assume

a) we're talking about the unilateral Laplace transform (because the IC's are given),

b) the initial conditions apply at t = 0- (typical when using the unilateral LT, and is also implicity assumed in Matlab implementation of laplace),

c) the forcing function on the right hand side is actually supposed to be a unit step function (because the unilateral LT is used for causal signals).

So my solution would be

syms t
y(t) = ((1/3) - (1/2)*exp(-t) + (1/6)*exp(-3*t))*heaviside(t);

which satisfies the IC at t = 0- and satisifes the differential equation (assuming the RHS is supposed to be a step)

simplify(diff(y(t),t,2) + 4*diff(y(t),t) + 3*y(t))

ans = 

[David Goodmanson]

Hi Paul,

That's all true. I assumed heaviside without even really thinking about it. On the other hand, the solution without the multiplicative heaviside function satisfies the diffeqn and meets the boundary conditions for all values of t, negative and positive. And, for example, in a standard reference like Schaum's Outlines (Murray Spiegel) in the Ch. 3 worked problems for situations exactly like this one, he gives the solutions without the multiplicative heaviside function. But your clarification is worthwhile and it would have been better had I thought enough to mention something to begin with.

[Paul]

David, I'm sure you know this but in case anyone else is interested ....

Under my assumption that the input is really a step function, including the heavisde(t) on the output will make a difference when differentating the expression for y(t) to get its derivatives. Taking the derivative enough times, in this problem it's three times, will start to yield Dirac deltas and Dirac derivatives.

Side note: I was pleasantly surprised to see that the definition of the unilateral Laplace transform in 2023a doc laplace shows the lower limit of the defining integral at t = 0-, which changed somewhere along the way from when it was shown as just t=0, e.g., in laplace 2018a

However, the doc page also states: " ... assumes that the signal f(t) is only defined for all real numbers t ≥ 0, or f(t) = 0 for t < 0." I don't agree with the "only defined" clause. I think it should just say " ... assumes that the signal f(t) satisfies f(t) = 0 for t < 0-."

[David Goodmanson]

Hi Paul,

I agree with all of this and the importance of t = 0-, except that on the last line, ' f(t) = 0 for t<0 ' is true most of the time but not in general. The classic case is the voltage across a capacitor that is initially charged to voltage V0 and discharged at t=0 with a switch through a resistor, showing exponential decay down to zero voltage, in which case f(t) = V0 for t<0.

[Paul]

Hi David,

The few references that I have handy discuss the class of functions operated on by the unilateral Laplace transform as either:

a) functions defined over the entire real line, but limited to be non-zero only for t>=0 (one reference uses t > 0),

or

b) functions defined over the (limited) domain t >= 0 (basically, negative time doesn't exist).

I've always wondered why no reference (that I've found) says something to the effect of the unilateral Laplace transform can be applied to the same class of functions that apply for the bilateral Laplace transform, but the unilateral Laplace transform just ignores the part of the function for t < 0. I've always assumed, but don't know, that there's a good reason for not stating it this way.

One thing that appeals to me using either a) or b) is that either condition maintains a 1-1 mapping from t to s and from s to t (at least I think it does). Whether or not that appeal is important mathematically is another question.

For the example you posed, would you write the solution for voltage as

v(t) = V0*exp(-t) , or

v(t) = V0*exp(-t)*heaviside(t) , or

v(t) = V0*exp(-t)*heaviside(t) + V0*heaviside(-t) , or

something else?

[David Goodmanson]

Hi Paul,

Possibly the reason that you have not seen a statement about unilateral being applied to bilateral situations is that the number of functions with a bilateral transform is very limited compared to functions with a unilateral transform.

For the capacitor I tend to think of this from physics point of view where all times exist, so -inf <= t <= inf. Before the switch is closed at t = 0, the capacitor is sitting at a constant V0 for all t<0. For a time constant RC = tau,

V = V0                t<0                   (a)
V = V0*exp(-t/tau)    t>=0

which in terms of heaviside is your last expression. If that's plotted for all times it makes sense physically.

Of course as you say there are other ways to look at this, one way being your second expression

V =  0                t<0                   (b)
V = V0*exp(-t/tau)    t>=0

where V(0) = V0 is an initial condition. If that's plotted for all times there is a big unphysical jump in voltage at t=0, but if one knows what is going on Laplacewise, it makes sense mathematically.

Contrast that with the situation of an inductor connected to a resistor in parallel with a switch, with the switch initially closed. For t<0 , a constant current I0 flows. When the switch is opened at t=0, the current is forced through the resistor. For a time constant L/R = tau, the voltage across the inductor is

V = 0                      t<0             (c)
V = I0*R*exp(-t/tau)       t>=0 

which has the same form as (b) but here the jump at t = 0 really is there physically.

Anyway I have to take back my previous comment objecting to f(t) = 0 for t<0 in general. I agree with you that one can always have that condition, although I also think it's misleading sometimes.

[Paul]

After reflection, I agree that whether or not the heaviside should be multiplied onto the solution of the non-zero initial value problem depends on what the underlying equations and the solution is intended to physically represent.

Answer 1

[Star Strider]

We do not know what you did to solve this. (Always post what you’ve already attempted.) However, if you got frustrated, that is easily understandable.

After a few frustrating minutes dealing with the Symbolic Math Toolbox Laplace transforms...

Try this:

syms s t T Y y(t) y0 Y(s) 
Dy = diff(y);
D2y = diff(Dy);
Eqn = D2y + 4*Dy + 3*y;
LapEqn = laplace(Eqn);
LapEqn = subs(LapEqn, {laplace(y(t), t, s), subs(diff(y(t), t), t, 0), y(0)},  {Y(s), 0, 0})
Ys(s) = simplify(LapEqn/Y)

producing finally (after all that):

Ys(s) =
s^2 + 4*s + 3

The time-domain solution (such as it is) is then:

yt = ilaplace(Ys,s,t)

yt =
3*dirac(t) + 4*dirac(1, t) + dirac(2, t)

If you are just beginning to use the MATLAB Symbolic Math Toolbox, I would not expect that you would be able to work with it in this context. This should quite definitely not be as difficult as it is.

[Gerard Nagle]

Hi Star Strider,

Great solution, i was spending a lot of time trying to work out the exactly the same question, and the MATLAB help and documented example is not really great. I do have a question on your answer.

In the question posed, where did the =1 on the RHS of the differential go in the answer?

What is the use of the y0 in the syms line.

I took your solution, applied it to the equation with inital conditions and and got out the answer.

syms s t T Y y(t) y0 Y(s) 
Dy = diff(y);
D2y = diff(Dy);
Eqn = D2y - Dy - 2*y == 5;
LapEqn = laplace(Eqn);
LapEqn = subs(LapEqn, {laplace(y(t), t, s), subs(diff(y(t), t), t, 0), y(0)},  {Y(s), 0, 0});
Ys(s) = simplify(LapEqn/Y);
pretty(Ys(s))

answer is

           2
s Y(s) (- s  + s + 2) == -5 and Y(s) ~= 0

which is or rearranged to which is the other option for the soluton, as done by hand.

My question is, do we know why MATLAB gives the warning ? something to do with zero inital conditions?

When I then try the inverse Lapalce,

yt = ilaplace(Ys,s,t)

I get the errors

Error using symengine
First argument must be of type 'Type::Arithmetical'.

Error in transform (line 74)
F = mupadmex('symobj::vectorizeSpecfunc', f.s, x.s, w.s, trans_func, 'infinity');

Error in sym/ilaplace (line 31)
F = transform('ilaplace', 's', 't', 'x', L, varargin{:});

the actual answer by hand is

[Konard Adams]

%% Solving 2nd order ODE using laplace transform

clc, clear , close all

%% Consider the initial value problem : y" - y' - 2 y = 5 , y(0) = 0 , y'(0) = 0

syms s t Y

% Defining the right-hand side function and find its Laplace transform

f = 5;

F = laplace(f,t,s);

% Finding the Laplace transform of y'(t) : Y1 = s Y - y(0)

Y1 = s*Y - 0;

% Finding the Laplace transform of y''(t) : Y2 = s Y1 - y'(0)

Y2 = s*Y1 - 0;

% Setting the Laplace transform of the left hand side minus the right hand side to zero and solve for Y:

Sol = solve(Y2 -Y1 - 2*Y - F, Y);

Solexp1= expand(Sol);

disp('Expanded Laplace Solution = '), disp(Solexp1)

Solsim1= simplify(Sol);

disp('Simplified Laplace Solution = '), disp(Solsim1)

% Finding the inverse Laplace transform of the solution

solin = ilaplace(Sol,s,t);

Solexp2= expand(solin);

disp('Expanded Inverse Laplace Solution = '), disp(Solexp2)

Solsim2= simplify(solin);

disp('Simplified Inverse Laplace Solution = '), disp(Solsim2)

Answer 2

[Konard Adams]

@Gerard Nagle I tried this code if it's any good and you can make it any better plesae let me know.

%% Solving 2nd order ODE using laplace transform

clc, clear , close all

syms x y t s X F

eqn = str2sym('diff(x(t),t,t)-diff(x(t),t)-2*y(t)= 5'); % Converting from string to an equation

F = laplace(eqn,s); % solving using laplace transform

eqn2 = str2sym('laplace(x(t),t,s)');

F = subs(F,{eqn2},{X}); % substituting the initial values then solve Laplace

eqn3 = str2sym('x(0)'); % Converting from string to an equation (initial values)

eqn4 = str2sym('Dx(0)'); % Converting from string to an equation (initial values)

F = subs(F,{eqn3,eqn4},{0,0}); % initial values

solx = solve(F,X);

solx = ilaplace(solx);

solxep = expand(solx); %pretty(x);

disp('Expanded Solution = '),disp(solxep)

solxsim = simplify(solx); %pretty(x);

disp('Simplified Solution = '),disp(solxsim)

Answer 3

[Jokita harini]

y''-6y'+8y=0,y(0)=-1,y'(0)=2

Answer 4

[Maadh]

120/s^6-1/(s-4)