Equivalent to Sage's assume: inverse Laplace of second order transfer function

[Carl Sandrock]

I am trying to build a workbook to illustrate the effect of various parameters of second order transfer functions. The full workbook is on GitHub, but here is a minimal example of the problem:

import sympy

tau, zeta, t, w, K = sympy.symbols('tau, zeta, t, w, K', real=True, positive=True)
s = sympy.Symbol('s')

G = K/(tau**2*s**2 + 2*tau*zeta*s + 1)

The impulse response of the second order system is simply the inverse Laplace of G. However, the nature of the inverse depends on the parameter ζ. TAttempting directly to calculate the inverse results in 

TypeError: cannot determine truth value of
-zeta/tau - sqrt(zeta + 1)*cos(atan2(0, zeta - 1)/2)*sqrt(Abs(zeta - 1))/tau < oo

There are three cases of interest, ζ>1, ζ=1 and 0<ζ<1
In Sage, I would be able to use assume(zeta > 1) before calculating the inverse to obtain the correct version of the inverse, but I have not found a way to impose such constraints in SymPy. So, first question is whether I can find nice solutions for these cases to the inverse.

Failing that, I want at least to be able to calculate the inverse with known values of all the parameters so that I can animate the response using IPython notebook widgets. Here I have also been out of luck, as some cases result in special values which are not cleanly evaluated to 

knownbadvalues = [{K: 5.05, tau: 5., zeta: 1.},
                  {K: 5.05, tau: 5.05, zeta: 1.05}
                  ]
for values in knownbadvalues:
    print sympy.inverse_laplace_transform(G.subs(values), s, t, noconds=True)
(0.202*t - 0.404*EulerGamma - 0.404*polygamma(0, 1.0))*exp(-0.2*t)
0.198019801980198*meijerg(((0.728681938243239, 0.855476477598345), ()), ((), (-0.271318061756761, -0.144523522401655)), exp(t))
ts = numpy.linspace(0, tmax, 100)
sympy.lambdify(t, invL(G.subs(values)), ['numpy', 'sympy'])(ts).n()


fails with "ValueError: sequence too large; must be smaller than 32"
Any advice on getting either getting the closed forms or just finding a version which can be evaluated cleanly?

[Francesco Bonazzi]

You may assume variables to be positive/negative, for now.

Sage uses a lot of backends. In case of assume(zeta > 1), I believe it relies on Maxima, which can handle that kind of assumptions.

SymPy's aim, on the other hand, is to implement all the algorithms on its own.

Assumptions with inequalities look like they are not yet supported, but I suppose they will be supported by a syntax similar to this one, as soon as the algorithm handling them will be finished:

from sympy.assumptions.assume import global_assumptions
global_assumptions.add(Q.positive(zeta - 1))

(Note: this code doesn't currently work as expected)

[Aaron Meurer]

This error usually indicates a bug in library code. Can you open an
issue with a full reproducing example and a full traceback?

Aaron Meurer

>>
>> -zeta/tau - sqrt(zeta + 1)*cos(atan2(0, zeta - 1)/2)*sqrt(Abs(zeta -
>> 1))/tau < oo
>>
>>
>> There are three cases of interest, ζ>1, ζ=1 and 0<ζ<1
>>
>> In Sage, I would be able to use assume(zeta > 1) before calculating the
>> inverse to obtain the correct version of the inverse, but I have not found a
>> way to impose such constraints in SymPy. So, first question is whether I can
>> find nice solutions for these cases to the inverse.
>>
>>
>> Failing that, I want at least to be able to calculate the inverse with
>> known values of all the parameters so that I can animate the response using
>> IPython notebook widgets. Here I have also been out of luck, as some cases
>> result in special values which are not cleanly evaluated to
>>
>>
>> knownbadvalues = [{K: 5.05, tau: 5., zeta: 1.},
>> {K: 5.05, tau: 5.05, zeta: 1.05}
>> ]
>> for values in knownbadvalues:
>> print sympy.inverse_laplace_transform(G.subs(values), s, t,
>> noconds=True)
>>
>> (0.202*t - 0.404*EulerGamma - 0.404*polygamma(0, 1.0))*exp(-0.2*t)
>> 0.198019801980198*meijerg(((0.728681938243239, 0.855476477598345), ()),
>> ((), (-0.271318061756761, -0.144523522401655)), exp(t))
>>
>> ts = numpy.linspace(0, tmax, 100)
>>
>> sympy.lambdify(t, invL(G.subs(values)), ['numpy', 'sympy'])(ts).n()
>>
>>
>> fails with "ValueError: sequence too large; must be smaller than 32"
>> Any advice on getting either getting the closed forms or just finding a
>> version which can be evaluated cleanly?
>

> --
> 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 post to this group, send email to sy...@googlegroups.com.
> Visit this group at http://groups.google.com/group/sympy.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/sympy/c731194a-56ad-4385-ab6a-b0e5514d74e0%40googlegroups.com.
>
> For more options, visit https://groups.google.com/d/optout.

[Kalevi Suominen]

On Thursday, September 3, 2015 at 7:19:38 PM UTC+3, Carl Sandrock wrote:

I am trying to build a workbook to illustrate the effect of various parameters of second order transfer functions. The full workbook is on GitHub, but here is a minimal example of the problem:

import sympy

tau, zeta, t, w, K = sympy.symbols('tau, zeta, t, w, K', real=True, positive=True)
s = sympy.Symbol('s')

G = K/(tau**2*s**2 + 2*tau*zeta*s + 1)

The impulse response of the second order system is simply the inverse Laplace of G. However, the nature of the inverse depends on the parameter ζ. TAttempting directly to calculate the inverse results in 

TypeError: cannot determine truth value of
-zeta/tau - sqrt(zeta + 1)*cos(atan2(0, zeta - 1)/2)*sqrt(Abs(zeta - 1))/tau < oo

This error message does not appear in the current master. However, the result is very complicated and not necessarily correct for all values of the parameters.
 

There are three cases of interest, ζ>1, ζ=1 and 0<ζ<1
In Sage, I would be able to use assume(zeta > 1) before calculating the inverse to obtain the correct version of the inverse, but I have not found a way to impose such constraints in SymPy. So, first question is whether I can find nice solutions for these cases to the inverse.

Failing that, I want at least to be able to calculate the inverse with known values of all the parameters so that I can animate the response using IPython notebook widgets. Here I have also been out of luck, as some cases result in special values which are not cleanly evaluated to 

knownbadvalues = [{K: 5.05, tau: 5., zeta: 1.},
                  {K: 5.05, tau: 5.05, zeta: 1.05}
                  ]

The results will be better if exact (rational) values are substituted. The integrator can not deal with floats in many cases.
 

[Kalevi Suominen]

It appears possible to compute the transform for ``zeta > 1`` by replacing ``zeta`` with ``1 + zeta``.

In [22]: G = K/(tau**2*s**2 + 2*tau*(1 + zeta)*s + 1)

In [24]: print(inverse_laplace_transform(G, s, t))
K*(exp(t*sqrt(zeta)*sqrt(zeta + 2)/tau) - 1)*(exp(t*sqrt(zeta)*sqrt(zeta + 2)/tau) + 1)*exp(-t*(sqrt(zeta)*sqrt(zeta + 2) + zeta + 1)/tau)*Heaviside(t)/(2*tau*sqrt(zeta)*sqrt(zeta + 2))


The transform of the original is obtained by replacing ``zeta`` with ``zeta - 1``.
       

[Carl Sandrock]

I have created an issue: https://github.com/sympy/sympy/issues/9900

 

It appears that this is a regression - the inverse is calculated appropriately in 0.7.5

[Carl Sandrock]

On Friday, September 4, 2015 at 6:28:08 PM UTC+2, Kalevi Suominen wrote:

On Thursday, September 3, 2015 at 7:19:38 PM UTC+3, Carl Sandrock wrote:

I am trying to build a workbook to illustrate the effect of various parameters of second order transfer functions. The full workbook is on GitHub, but here is a minimal example of the problem:

import sympy

tau, zeta, t, w, K = sympy.symbols('tau, zeta, t, w, K', real=True, positive=True)
s = sympy.Symbol('s')

G = K/(tau**2*s**2 + 2*tau*zeta*s + 1)

The impulse response of the second order system is simply the inverse Laplace of G. However, the nature of the inverse depends on the parameter ζ. TAttempting directly to calculate the inverse results in 

TypeError: cannot determine truth value of
-zeta/tau - sqrt(zeta + 1)*cos(atan2(0, zeta - 1)/2)*sqrt(Abs(zeta - 1))/tau < oo

This error message does not appear in the current master. However, the result is very complicated and not necessarily correct for all values of the parameters.
 

There are three cases of interest, ζ>1, ζ=1 and 0<ζ<1
In Sage, I would be able to use assume(zeta > 1) before calculating the inverse to obtain the correct version of the inverse, but I have not found a way to impose such constraints in SymPy. So, first question is whether I can find nice solutions for these cases to the inverse.

Failing that, I want at least to be able to calculate the inverse with known values of all the parameters so that I can animate the response using IPython notebook widgets. Here I have also been out of luck, as some cases result in special values which are not cleanly evaluated to 

knownbadvalues = [{K: 5.05, tau: 5., zeta: 1.},
                  {K: 5.05, tau: 5.05, zeta: 1.05}
                  ]

The results will be better if exact (rational) values are substituted. The integrator can not deal with floats in many cases.

The problem seems to go away when using Sympy.Rational(value*1000, 1000). This is a decent workaround for now, thank you.