Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

slow inverse_laplace_transform for easy poles in LH plane #13606

Open
cbm755 opened this issue Nov 14, 2017 · 11 comments

Comments

@cbm755
Copy link
Contributor

cbm755 commented Nov 14, 2017

Consider:

>>> f = 1/((p+1)*(p+8))
         1       
  ───────────────
  (p + 1)⋅(p + 8)
>>> inverse_laplace_transform(f, p, t)   % takes maybe 2 seconds
⎛ 7⋅t    ⎞  -8⋅t             
⎝ℯ    - 1⎠⋅ℯ    ⋅Heaviside(t)
─────────────────────────────
              7              

But change it to f = 1/((p+1)*(p+10)) and it takes maybe 10 seconds.

Then change it to have ...(p+15) and it takes at least 5 minutes (killed it then).


These functions have simple poles in the left have plane. They should be easy to invert.

@cbm755

This comment has been minimized.

Copy link
Contributor Author

cbm755 commented Nov 14, 2017

I labeled this "easy to fix" but... Well, its probably easy to get involved and potentially a nice project for someone who has learned their laplace transforms, but other than that I have no idea..., it might be a deep dive into the transform code.

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 14, 2017

This probably isn't easy to fix. The integral transforms are currently computed with the meijerg algorithm. If you don't know what's going on there, it could take some time to learn it.

What we should do is add some table lookups a la manualintegrate for the common integral transforms so that cases like this can be computed quickly.

@jksuom

This comment has been minimized.

Copy link
Member

jksuom commented Nov 15, 2017

It seems that most of the time is spent in hyperexpand. The inverse transform of 1/((p+1)*(p+10)) is quickly found to involve meijerg(((-9, 0), ()), ((), (-10, -1)), z). (Note the poles in bq and counterterms in ap.) Then the expression is simplified using hyperexpand:

In [15]: f = meijerg(((-9, 0), ()), ((), (-10, -1)), z)

In [16]: %time hyperexpand(f)
CPU times: user 28.46 s, sys: 0.08 s, total: 28.54 s
Wall time: 28.54 s

It looks like the expansion process is computing with large (thousands of ops) polynomial expressions. It should be investigated if those could be replaced by Polys.

@cbm755 cbm755 removed the Easy to Fix label Nov 27, 2017
@vishalg2235

This comment has been minimized.

Copy link
Contributor

vishalg2235 commented Nov 30, 2017

hey! I would like to work on this issue. @cbm755 can I work on this issue ?

@cbm755

This comment has been minimized.

Copy link
Contributor Author

cbm755 commented Nov 30, 2017

Anyone can work on any issue, no permission required :-)

@vishalg2235

This comment has been minimized.

Copy link
Contributor

vishalg2235 commented Nov 30, 2017

Thank you!
can you help me how to fix it ? which part of file meijerint.py I should look for ?

@cbm755

This comment has been minimized.

Copy link
Contributor Author

cbm755 commented Nov 30, 2017

Sorry, I cannot. Note also @asmeurer's note earlier comment #13606 (comment)

@vishalg2235

This comment has been minimized.

Copy link
Contributor

vishalg2235 commented Nov 30, 2017

ok @cbm755 .
@asmeurer how to add a lookup table ?

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 30, 2017

I would look at adding transform integrals to manualintegrate. However, I'm not sure if manualintegrate supports integration limits yet (maybe it does, I'm not sure). If it doesn't, that will need to be added first.

@vishalg2235

This comment has been minimized.

Copy link
Contributor

vishalg2235 commented Dec 1, 2017

@asmeurer is there any work here in which I could help ?

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Dec 1, 2017

I'm not aware of existing work but you might search the open issues and pull requests to see what is there.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
4 participants
You can’t perform that action at this time.