Python and Laplace transforms of an IVP

[Staffan Lundberg]

I am working with a project to replace Matlab with Python, in a calculus course. Explicitly, following problem is solved in Matlab. It is an inhomogenous IVP with constant coefficients together with  relevant IC's. In Matlab Code: Dy=diff(y(t),t); D2y=diff(Dy,t); g=8*t*(heaviside(t)-heaviside(t-5))+40*heaviside(t-5); ode=D2y+4*y-g ==0; My goal is to Laplace transform the equation, obtain a solution in frequency domain, and finally transform it back to time domain, obtaining the final solution. This procedure is successfully done by Matlab.

I address this issue to the developer team. Will it, in future releases of SymPy, be possible to solve this problem in Python/SymPy?

[Oscar Benjamin]

I haven't used sympy's Laplace transform code but I can tell you that
dsolve is able to solve the differential equation:

In [14]: y = Function('y')

In [15]: t = Symbol('t')

In [16]: g = 8*t*(Heaviside(t) - Heaviside(t-5)) + 40*Heaviside(t-5)

In [17]: eq = Eq(y(t).diff(t, 2) + 4*y(t), g)

In [18]: str(eq)
Out[18]: 'Eq(4*y(t) + Derivative(y(t), (t, 2)), 8*t*(Heaviside(t) -
Heaviside(t - 5)) + 40*Heaviside(t - 5))'

In [19]: str(dsolve(eq))
Out[19]: 'Eq(y(t), 2*(t*Heaviside(t) + 5*Heaviside(t - 5))*cos(2*t)**2
+ (C1 - 10*cos(10)*Heaviside(t - 5) +
4*Integral(t*sin(2*t)*Heaviside(t - 5), t))*cos(2*t) +
2*(t*Heaviside(t) - t*Heaviside(t - 5) + 5*Heaviside(t -
5))*sin(2*t)**2 + (C2 - cos(2*t)*Heaviside(t - 5) - Heaviside(t) +
cos(10)*Heaviside(t - 5))*sin(2*t))'

There is an unevaluated integral in there and it's quite slow because
of computing the integrals involved so there can be improvements to
made to the integration routines here.

Not having used sympy's laplace transforms before I'm not sure how
they would be used for a case like this:

In [24]: s = Symbol('s')

In [25]: str(laplace_transform(eq.lhs - eq.rhs, t, s))
Out[25]: '(4*LaplaceTransform(y(t), t, s) +
LaplaceTransform(Derivative(y(t), (t, 2)), t, s) - 40*exp(-5*s)/s +
8*(5*s - exp(5*s) + 1)*exp(-5*s)/s**2, 0, True)'

I guess that `LaplaceTransform(Derivative(y(t), (t, 2)), t, s)` is an
unevaluated Laplace transform of the derivative. How exactly does
Matlab represent that?


Oscar

> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/10568891-8a24-43aa-8f93-e1a6b7de72e6n%40googlegroups.com.

[Staffan Lundberg]

Yes Oscar, I am fully aware of dsolve.

Here is the successful Matlab code (Note: IC y(0)=y’(0)=0)

syms s y(t) Y g(t)

Dy=diff(y(t),t);

D2y=diff(Dy,t);

g=8*t*(heaviside(t)-heaviside(t-5))+40*heaviside(t-5);

ode=D2y+4*y-g ==0;

ltode=laplace(ode);

ltode=subs(ltode,laplace(y),Y);

ltode=subs(ltode,y(0),0);

ltode=subs(ltode,subs(diff(y(t), t), t, 0),0);

Y=solve(ltode,Y);

y=ilaplace(Y)

Otput from Matlab on input ltode=laplace(ode);

tode =

 

(8*exp(-5*s))/s^2 - s*y(0) + s^2*laplace(y(t), t, s) - subs(diff(y(t), t), t, 0) - 8/s^2 + 4*laplace(y(t), t, s) == 0

ltode=subs(ltode,laplace(y),Y)  gives output

ltode =

 

4*Y - s*y(0) + (8*exp(-5*s))/s^2 - subs(diff(y(t), t), t, 0) + Y*s^2 - 8/s^2 == 0

Matlab output on  commands

ltode=subs(ltode,y(0),0)

ltode=subs(ltode,subs(diff(y(t), t), t, 0),0)

gives

ltode =

 

4*Y + (8*exp(-5*s))/s^2 - subs(diff(y(t), t), t, 0) + Y*s^2 - 8/s^2 == 0

 

ltode =

 

4*Y + (8*exp(-5*s))/s^2 + Y*s^2 - 8/s^2 == 0

Y=solve(ltode,Y)  produces output

Y =

 

-((8*exp(-5*s))/s^2 - 8/s^2)/(s^2 + 4)

and finally

y=ilaplace(Y)

gives

y =

 

2*t - sin(2*t) + 8*heaviside(t - 5)*(sin(2*t - 10)/8 - t/4 + 5/4)

 

So, chis procedure is not possible to process in current version of SymPy.

When will it be possible? Any estimation?

Cheers,

Staffan

To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxRJUFhiF3t9qtsMBTqhMBqg8JPAOaV2gyTjdHCYdSs0MA%40mail.gmail.com.

[Oscar Benjamin]

As I said I'm not familiar with Laplace transform capabilities of
sympy myself so I'm not sure if there already is a way to do that.

It looks from the matlab output you show compared to the sympy output
I showed like there is a missing piece in sympy which is that the
Laplace transform of a derivative should expand in terms of the
initial conditions. If that's the case then it should be pretty easy
to add that. Specifically I mean that
LaplaceTransform(Derivative(y(t), t, 2)) should expand to something
involving y(0) and Subs(Derivative(y(t), t), t, 0) as it seems to in
the Matlab output.

As for "when" that will happen the answer is just that it will happen
when someone adds the feature. I'd be happy to answer specific
questions to help someone implement this if anyone is interested but
it should be fixed by someone with more experience of Laplace
techniques than me so that they are clear about what they are aiming
for (or at least the fix should be reviewed by someone like that).

There are quite a few enthusiastic newer contributors on github right
now who are generally looking for things to work on. If you were to
open an issue that explains in (very) clear terms what should be done
and what should be made to work then I expect that some of them would
be interested to do this. Sympy is not generally short of enthusiastic
contributors: what we need more of is people with experience either
with sympy or in a given problem domain who can guide and review their
contributions.

Alternatively you could implement this yourself. I don't expect this
to be a difficult feature to add but I'm not sure if there would then
be other changes needed to make your full calculation work.

Oscar

[emanuel.c...@gmail.com]

with 1.7.1, I get simpler results (no unsolved integrals) :

>>> y=Function("y")
>>> t=Symbol("t")
>>> g=8*(Heaviside(t)-Heaviside(t-5))+40*Heaviside(t-5)
>>> eq=Eq(y(t).diff(t,2)+4*y(t),g)
>>> dsolve(eq,y(t))
Eq(y(t), C1*sin(2*t) + C2*cos(2*t) + 4*sin(t)**2*Heaviside(t) + 16*sin(t)*sin(t - 10)*Heaviside(t - 5) - 8*cos(10)*Heaviside(t - 5) + 8*Heaviside(t - 5))
>>> dsolve(eq,y(t),ics={y(0):0,y(t).diff(t).subs(t,0):0})
Eq(y(t), 4*sin(t)**2*Heaviside(t) + 16*sin(t)*sin(t - 10)*Heaviside(t - 5) - 8*cos(10)*Heaviside(t - 5) + 8*Heaviside(t - 5))

And the result is consonant with whatever Mathematica returns…

HTH,

​

[Oscar Benjamin]

On Fri, 15 Jan 2021 at 09:30, emanuel.c...@gmail.com
<emanuel.c...@gmail.com> wrote:
>
> with 1.7.1, I get simpler results (no unsolved integrals) :

Oh yes. Thanks for pointing that out. Actually I checked now and I
also get the same as you with both 1.7.1 and the current master
branch.

I'm not sure what I was using that gave the unevaluated integrals.
Perhaps I was testing on some development branch that wasn't fully
working...


Oscar

[Aaron Meurer]

There's a more general issue with all the integral transforms in SymPy
(Laplace, Fourier, etc.) that they don't do transforms based on basic
table lookups. Right now they only do transforms based on computing
the integral using integrate(), meaning if the integrator cannot
compute an integral, the transform fails. Here's a discussion of the
idea https://groups.google.com/u/1/g/sympy/c/WkgwYAFC_P4/m/5tA2ZsOCCAAJ.
As I noted in that discussion, it would be best if the improvements
could be done to integrate() itself rather than the transform code,
that way the transform will work even if someone types it as a direct
integral (and even if they don't realise their integral is a common
integral transform it will work). Probably the best place to start
with this is manualintegrate, although improving meijerg would also be
good (but much more complicated as that requires a deep understanding
of that algorithm).

A good place to start, as Oscar mentioned, would be to open an issue
for this specific problem (Laplace transform of a derivative). As far
as I can tell one doesn't already exist.

If we can make Laplace transforms more robust in SymPy we should be
able to add an ODE solver that uses them, so that what you are trying
to do would be equivalent to dsolve(hint='laplace').

Aaron Meurer

> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxTBS6OoK9FtWySKuF88ZLrb-Xm6f8egJ%2B0U3-QGxCK_2w%40mail.gmail.com.