Slow inverse_laplace_transformation

83 views
Skip to first unread message

Sten Sogaard

unread,
Feb 7, 2015, 5:06:46 PM2/7/15
to sy...@googlegroups.com
I am trying to use inverse_laplace_transformation for some relative simple filter calculations.

I works fine with very simple 1st and 2nd order filters, but if it is more complex it takes forever.

Even this very simple polynomial takes a long time (30 minutes and still no solution).

from sympy import *
from sympy.abc import s,t

H = 1/(0.07071*s**2 + 258.0*s + 1100000.)
inverse_laplace_transform(1/s*H,s,t)


Any ideas?

Aaron Meurer

unread,
Feb 7, 2015, 11:18:00 PM2/7/15
to sy...@googlegroups.com
In SymPy 0.7.6 as well as the git master for me it immediately gives nan (which I guess is a bug).

Aaron Meurer

--
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/bb31e487-58e6-48cc-b682-3e62ae660224%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

Kalevi Suominen

unread,
Feb 8, 2015, 3:57:26 PM2/8/15
to sy...@googlegroups.com


On Sunday, February 8, 2015 at 12:06:46 AM UTC+2, Sten Sogaard wrote

H = 1/(0.07071*s**2 + 258.0*s + 1100000.)
inverse_laplace_transform(1/s*H,s,t)


Any ideas?
 
The current implementation leads to floating point numbers with very large positive and negative exponents. After the small numbers are rounded to zero, nan results.

However, the computation succeeds with suitably normalized coefficients:
```
In [38]: H1 = H.subs(s, 10000.*s)
In [40]: print(inverse_laplace_transform(1/s*H1, s, t))
-1.0*(-9.09090909090909e-7*exp(0.182435299109037*t) + 4.74279450417556e-7*sin(0.349688926583877*t) + 9.09090909090909e-7*cos(0.349688926583877*t))*exp(-0.182435299109037*t)*Heaviside(t)
```
(Of course, the result has to be normalized again.)

For a more satisfactory solution, cancellation of exponents in the intermediate results would be required.

Aaron Meurer

unread,
Mar 11, 2015, 12:59:33 AM3/11/15
to sy...@googlegroups.com
So perhaps you would have better results if you started with a Float
with higher precision, like Float("0.07071", 100).

Aaron Meurer
> --
> 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/83a765a9-c984-4fa2-bbc2-00010e02ed83%40googlegroups.com.

Kalevi Suominen

unread,
Mar 11, 2015, 3:38:22 AM3/11/15
to sy...@googlegroups.com


On Wednesday, March 11, 2015 at 6:59:33 AM UTC+2, Aaron Meurer wrote:
So perhaps you would have better results if you started with a Float
with higher precision, like Float("0.07071", 100).

The floats are of the order  10**(-19000)  in both the numerator and the denominator before rounding.
 

Aaron Meurer

unread,
Mar 16, 2015, 3:54:34 PM3/16/15
to sy...@googlegroups.com
On Wed, Mar 11, 2015 at 2:38 AM, Kalevi Suominen <jks...@gmail.com> wrote:
>
>
> On Wednesday, March 11, 2015 at 6:59:33 AM UTC+2, Aaron Meurer wrote:
>>
>> So perhaps you would have better results if you started with a Float
>> with higher precision, like Float("0.07071", 100).
>
>
> The floats are of the order 10**(-19000) in both the numerator and the
> denominator before rounding.

Oh, that's a bit high. Is it possible to do the stabilization automatically?

Aaron Meurer
> https://groups.google.com/d/msgid/sympy/4b49a4ae-ec3b-4ba5-acec-425e65865c3f%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages