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

Laplace transform involving time-shifts #22422

Closed
mforets mannequin opened this issue Feb 23, 2017 · 50 comments
Closed

Laplace transform involving time-shifts #22422

mforets mannequin opened this issue Feb 23, 2017 · 50 comments

Comments

@mforets
Copy link
Mannequin

mforets mannequin commented Feb 23, 2017

Sage allows to compute the inverse Laplace transform through Maxima's ilt function,

    sage: var('s t')
    sage: inverse_laplace(1/s, s, t)
    1

An unevaluated expression is returned when no explicit inverse Laplace transform is computed, as in

    sage: inverse_laplace(exp(-s)/s, s, t)
    ilt(e^(-s)/s, s, t)

The result in this case is h(t-1), where h is the Heaviside step function. In Sage it is available as heaviside.

The problem in this ticket is to extend the current behavior of inverse_laplace to provide explicit expressions for proper real-rational functions with any number of real exponentials linear in the transform variable s (time-shifts) in the numerator. For consistency, the direct Laplace transform with a heaviside should also work as well.

These are some approaches:

(1) Implement an in-house solution, possibly in the lines of this answer.

(2) Add an algorithm flag that allows to choose sympy (similar to integration).

(3) Interface with Giac/XCAS. With this package installed, it is possible to do:

sage: giac('invlaplace(exp(-s)/s, s, t)')
Heaviside(t-1)

IMHO, a combination of (2)-(3) is the more robust approach. A small set of experiments show that (3) is, at the time of writing, more convenient than inverse_laplace_transform of SymPy in terms of quality of solution and execution time. Unfortunately, the giac interface does not currently support automatic translation back to the symbolic ring, as it does with SymPy objects via SR(..).

Any recommendations?

See also:

CC: @kcrisman @paulmasson @frederichan-IMJPRG @rwst

Component: calculus

Keywords: laplace, transform, symbolics, giac, heaviside

Author: Marcelo Forets

Branch/Commit: f845269

Reviewer: Paul Masson, Ralf Stephan

Issue created by migration from https://trac.sagemath.org/ticket/22422

@mforets mforets mannequin added this to the sage-7.6 milestone Feb 23, 2017
@mforets mforets mannequin added c: calculus labels Feb 23, 2017
@mforets
Copy link
Mannequin Author

mforets mannequin commented Feb 27, 2017

@mforets
Copy link
Mannequin Author

mforets mannequin commented Feb 27, 2017

Commit: 2096388

@mforets
Copy link
Mannequin Author

mforets mannequin commented Feb 27, 2017

comment:3

Needs info (please see commit message) to transform heaviside via _giac_(). Example:

    sage: f = heaviside(x); f, type(f)
    (heaviside(x), <type 'sage.symbolic.expression.Expression'>)
    sage: fg = f._giac_(); fg,  type(fg)
    (heaviside(x), <class 'sage.interfaces.giac.GiacElement'>)

But heaviside(x) doesn't seem to be understood by giac.


New commits:

2096388Added algorithm='sympy' (needs work) and algorithm='giac' (needs info).

@mforets

This comment has been minimized.

@mforets
Copy link
Mannequin Author

mforets mannequin commented Feb 27, 2017

Author: Marcelo Forets

@mforets mforets mannequin added the s: needs info label Feb 27, 2017
@rwst
Copy link

rwst commented Feb 27, 2017

comment:4

Code looking good so far but needs doctests for all code branches. Minor: you seem to change spacing, there are empty hunks in the patch, check your editor configuration.

@frederichan-IMJPRG
Copy link

comment:5

Replying to @mforets:

Needs info (please see commit message) to transform heaviside via _giac_(). Example:

    sage: f = heaviside(x); f, type(f)
    (heaviside(x), <type 'sage.symbolic.expression.Expression'>)
    sage: fg = f._giac_(); fg,  type(fg)
    (heaviside(x), <class 'sage.interfaces.giac.GiacElement'>)

But heaviside(x) doesn't seem to be understood by giac.

in giac it looks to be Heaviside, so you need to either translate the string or define it in giac:

giac("'heaviside:=Heaviside'")
sage: giac(f).integrate(x,-1,2)
2
sage: f
heaviside(x)
sage: type(f)
<type 'sage.symbolic.expression.Expression'>

similar with giacpy_sage which have a raw conversion to sage:

sage: from giacpy_sage import *
// Giac share root-directory:/home/fred-dev/sage/develop/sage.develop/local/share/giac/
// Giac share root-directory:/home/fred-dev/sage/develop/sage.develop/local/share/giac/
Help file /home/fred-dev/sage/develop/sage.develop/local/share/giac/doc/fr/aide_cas not found
Added 0 synonyms
sage: libgiac('heaviside:=Heaviside')
'Heaviside'
sage: f=heaviside(x)
sage: fg=libgiac(f)
sage: fg.integrate(x,-1,2)
2
sage: type(fg.integrate(x,-1,2))
<type 'giacpy_sage.Pygen'>
sage: SR(fg.integrate(x,-1,2))
2
sage: type(SR(fg.integrate(x,-1,2)))
<type 'sage.symbolic.expression.Expression'>
sage: (fg.integrate(x,-1,2)).sage()
2
sage: type((fg.integrate(x,-1,2)).sage())
<type 'sage.rings.integer.Integer'>

but with Heaviside will be sent to SR as a raw string so you will need some translation back to heaviside also.

@mforets
Copy link
Mannequin Author

mforets mannequin commented Mar 8, 2017

comment:6

@rwst: ok, I am learning such things as doctesting in days84. for the spacing issue good that you pointed this out, I'll be more careful. BTW I use Atom editor, and thanks for the prompt feedback.

@frederichan-IMJPRG: thanks for the insight! I see that I could use the conversion at the script level. but then each new Sage function which for some reason needs it, would have to handle the conversion. does it make sense to define this at the level of the giac.py interface?

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 10, 2017

Changed commit from 2096388 to a3a43d2

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 10, 2017

Branch pushed to git repo; I updated commit sha1. New commits:

a3a43d2A solution for heaviside <-> Heaviside.

@mforets
Copy link
Mannequin Author

mforets mannequin commented Mar 10, 2017

comment:8

this is a patch that seems to work :)

  • borrowed the un_camel method from the mathematica interface, hence conversion back to lowercase heaviside is not done in place
  • extended the conversion dictionary at functions/generalized.py, it does the job when calling giac()
  • to-do: add doctests

@mforets mforets mannequin added s: needs work and removed s: needs info labels Mar 10, 2017
@paulmasson
Copy link
Mannequin

paulmasson mannequin commented Mar 13, 2017

comment:9

Replying to @mforets:

this is a patch that seems to work :)

  • borrowed the un_camel method from the mathematica interface, hence conversion back to lowercase heaviside is not done in place

This works for heaviside but wil create problems for Sage functions with capital letters like bessel_J.

@mforets
Copy link
Mannequin Author

mforets mannequin commented Mar 13, 2017

comment:10

Replying to @paulmasson:

Replying to @mforets:

this is a patch that seems to work :)

  • borrowed the un_camel method from the mathematica interface, hence conversion back to lowercase heaviside is not done in place

This works for heaviside but wil create problems for Sage functions with capital letters like bessel_J.

ok, good. so i suppose the correct way is to implement a basic parser, as suggested in giac.py interface (my comments in #):

    def _sage_(self):
        r"""
        Convert a giac expression back to a Sage expression.

        This currently does not implement a parser for the Giac output language,
        therefore only very simple expressions will convert successfully.
        Warning: List conversion is slow.
        ...

        """
        result = repr(self) # this is a string representation
        if str(self.type()) != 'DOM_LIST' :
            try:
                from sage.symbolic.all import SR
                result = giac2sage(result) # pattern matching e.g. Heaviside -> heaviside
                return SR(result)
            except Exception:
                raise NotImplementedError("Unable to parse Giac output: %s" % result)
        else:
            return [entry.sage() for entry in self]

a first search for this kind of conversion in some other module didn't give me anything useful yet.

@mforets
Copy link
Mannequin Author

mforets mannequin commented Mar 14, 2017

comment:11

for the new keyword argument's name, i would vote for engine instead of algorithm.
the former is used for example in base.py. the keyword algorithm is misleading: a given computational engine may implement different algorithms.

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 15, 2017

Branch pushed to git repo; I updated commit sha1. New commits:

3732a8cadd new doctests and a very basic parse function at giac.py

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 15, 2017

Changed commit from a3a43d2 to 3732a8c

@mforets mforets mannequin added s: needs review and removed s: needs work labels Mar 15, 2017
@paulmasson
Copy link
Mannequin

paulmasson mannequin commented Mar 19, 2017

comment:14

Replying to @mforets:

for the new keyword argument's name, i would vote for engine instead of algorithm.
the former is used for example in base.py. the keyword algorithm is misleading: a given computational engine may implement different algorithms.

While the term engine is more accurate, algorithm is already established for a variety of commands, such as integrate and limit. Isn't consistency in Sage more important? Unless you're proposing we change all such instances to engine.

@paulmasson
Copy link
Mannequin

paulmasson mannequin commented Mar 19, 2017

comment:15

In the function _giac2sage, rather than recreating the conversion dictionary why not use symbol_table from sage.libs.pynac.pynac? That dictionary already contains all the defined conversions.

I'm also getting the Unable to parse Giac output error often. Partly that's because more Giac conversions need to be defined. Will that happen on this ticket or a separate one?

If _un_camel is no longer used, please remove it.

@mforets mforets mannequin added s: needs work and removed s: needs review labels Mar 20, 2017
@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 21, 2017

Branch pushed to git repo; I updated commit sha1. New commits:

1454679delete unused un_camel function
dc35047use keyword algorithm in calculus.py

@paulmasson
Copy link
Mannequin

paulmasson mannequin commented Mar 29, 2017

comment:26

Replying to @frederichan-IMJPRG:

NB: is the local dictionnary usefull?
or is it better to keep conversions in one place. It seems one can always do this

sage.libs.pynac.pynac.register_symbol(heaviside,{'giac':'Heaviside'})

if some keyword is missing.

I agree with Frederic on this: if there is an existing method for extending the Pynac dictionary, then we should use that instead of introducing an additional dictionary. This method should probably be added to the documentation in calculus.py.

FYI your documentation doesn't build right now because of unbalanced double backticks. Also, single backticks are what we use for inline LaTeX rendering in ReST, not dollar signs.

@frederichan-IMJPRG
Copy link

comment:27

Replying to @paulmasson:

Replying to @mforets:

Replying to @mforets:

Replying to @paulmasson:

How does the new commit look? In particular, do you suggest to add more functionality relevant for this ticket? Can you supply an example with the Unable to parse Giac output message? Thanks.

Where I'm seeing this error often is when an expression has extra powers of the variable of integration, as for example

inverse_laplace(s^2*exp(-s)/(s-1),s,t,algorithm='giac')

which returns a derivative of the Dirac delta. At some point we'll want to be able to parse an indefinite number of such derivatives: do you want to do that on this ticket or another?

I'm also getting the Unable to parse Giac output error often. Partly that's because more Giac conversions need to be defined. Will that happen on this ticket or a separate one?

I suggest that we move on with this ticket for the basic functionality described above, and to create separate tickets to enhance the symbolics with the giac interface in separate ones (also add limit, integrate, solve). If you agree we can create a new subsection "Giac interface" under https://trac.sagemath.org/wiki/symbolics

Sounds fine. I'd be willing to include more Giac translations for special functions on a separate ticket, if you can point me to where the official list of Giac names is located.

In the file:
local/share/giac/doc/aide_cas
lines starting by #
are keywords followed by synonyms, next lines are for small doc and see also.
for full doc it is there:

http://www-fourier.ujf-grenoble.fr/~parisse/giac/doc/en/cascmd_en/cascmd_en.html

@mforets
Copy link
Mannequin Author

mforets mannequin commented Mar 29, 2017

comment:28

Replying to @frederichan-IMJPRG:

giacpy_sage is an optional spkg, but you should be able to install it.

All the files in interfaces are pexpect interfaces for external programs. Among them some (ex: pari, gap, singular) also have a cython interface that works with the c or C++ library directly. There was a wish to prefer the cython interface when avaible.

But as giacpy_sage is optional while giac went standard recently, and as calculus have old entries with the pexepect interface of giac, I think that it is not a problem to do work like this.
But indeed it is not natural to work on the pexpect interface too much.

the cython interface to giac is giacpy_sage and it superseeds the pexpect one, but your translation table could be imported by both. So when things will be decided I will update giacpy_sage also if it doesn't slow down conversions of huge expressions.

great, thanks for the detailed answer!

i wasn't aware about register_symbol, yes we should use it.

on the other hand, i tried to see how locals dictionary is used in other interfaces. i have the impression that it may be handy in some use cases, in the interactive mode. this example is from mathematica.py:

compare

sage: ex = giac('myFun(x)')
sage: ex._sage_({'myFun': sin})
sin(x)

to

sage: ex = giac('myFun(x)')
sage: sage.libs.pynac.pynac.register_symbol(sin, {'giac':'myFun'})
sage: ex._sage_()
sin(x)

the long import list is.. complicated.

@mforets
Copy link
Mannequin Author

mforets mannequin commented Mar 29, 2017

comment:29

Replying to @paulmasson:

Replying to @mforets:

Replying to @mforets:

Replying to @paulmasson:

How does the new commit look? In particular, do you suggest to add more functionality relevant for this ticket? Can you supply an example with the Unable to parse Giac output message? Thanks.

Where I'm seeing this error often is when an expression has extra powers of the variable of integration, as for example

inverse_laplace(s^2*exp(-s)/(s-1),s,t,algorithm='giac')

which returns a derivative of the Dirac delta. At some point we'll want to be able to parse an indefinite number of such derivatives: do you want to do that on this ticket or another?

I'm also getting the Unable to parse Giac output error often. Partly that's because more Giac conversions need to be defined. Will that happen on this ticket or a separate one?

I suggest that we move on with this ticket for the basic functionality described above, and to create separate tickets to enhance the symbolics with the giac interface in separate ones (also add limit, integrate, solve). If you agree we can create a new subsection "Giac interface" under https://trac.sagemath.org/wiki/symbolics

Sounds fine. I'd be willing to include more Giac translations for special functions on a separate ticket, if you can point me to where the official list of Giac names is located. I'd also like to do some code cleanup on nearby lines pertaining to LaTeX representations in Sage, so if that doesn't bother you we can kill two birds with one stone.

OK. Please see the new ticket #22706. I also created a new subsection at symbolics wiki to group these tickets.

Good finding about the missing support for derivatives of dirac delta! In my opinion it is ok to solve it in a separate ticket, since originally we begun to discuss about transforming proper functions (in the jargon of control theory, meaning that deg numerator <= deg denominator).

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 29, 2017

Branch pushed to git repo; I updated commit sha1. New commits:

c1f20cbdocstring improvements
4f347e2fix another typo in docstring

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 29, 2017

Changed commit from 4aae4f5 to 4f347e2

@mforets mforets mannequin added s: needs review and removed s: needs work labels Mar 30, 2017
@mforets
Copy link
Mannequin Author

mforets mannequin commented Mar 30, 2017

comment:32

Oh, and we discussed previously on the locals dictionary at giac.py interface, and the consensus was to remove it. Could you please confirm? Thanks

@frederichan-IMJPRG
Copy link

comment:33

IF you like this locals then it is OK for me, but a comment in giac.py code to remind and recommend how to use the table would be nice.

by the way why not just:
result = giac.laplace(ex, s, t)
(it looks that the conversion is also done)

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 30, 2017

Branch pushed to git repo; I updated commit sha1. New commits:

5936013simplify calls to giac
0577b5badd register_symbol explanation and example to giac.py

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 30, 2017

Changed commit from 4f347e2 to 0577b5b

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 30, 2017

Branch pushed to git repo; I updated commit sha1. New commits:

f597902minor docstring modif

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Mar 30, 2017

Changed commit from 0577b5b to f597902

@paulmasson
Copy link
Mannequin

paulmasson mannequin commented Apr 1, 2017

comment:36

Doctests all pass and documentation builds.

Single backticks typeset math so they render regular words in italic. All of the code terms in giac.py should be inside double backticks. Similarly, in calculus.py the word "cond" should be inside double backticks and one instance of "F" needs single backticks.

The documentation added to giac.py doesn't appear in the built documents because _sage_ is a private method. Is this a deliberate choice? How does one access the documentation for this method?

@paulmasson paulmasson mannequin modified the milestones: sage-7.6, sage-8.0 Apr 1, 2017
@rwst
Copy link

rwst commented Apr 5, 2017

comment:37

Replying to @paulmasson:

Doctests all pass and documentation builds.

Single backticks typeset math so they render regular words in italic. All of the code terms in giac.py should be inside double backticks. Similarly, in calculus.py the word "cond" should be inside double backticks and one instance of "F" needs single backticks.

Apart from these cosmetics I think this ticket is fine.

The documentation added to giac.py doesn't appear in the built documents because _sage_ is a private method. Is this a deliberate choice? How does one access the documentation for this method?

As far as I know there is no way to display such docstrings (I have struggled myself with this in the past). If something important is there it should be moved to the global documentation at the top of the file.

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Apr 6, 2017

Changed commit from f597902 to f845269

@sagetrac-git
Copy link
Mannequin

sagetrac-git mannequin commented Apr 6, 2017

Branch pushed to git repo; I updated commit sha1. New commits:

e4fc9e0fix some double backticks
8267dfdadded conversions and cleanup to giac.py doc
caae8dfadd 2 seealso blocks
f845269fix authors string

@mforets
Copy link
Mannequin Author

mforets mannequin commented Apr 6, 2017

comment:39

Replying to @paulmasson:

Doctests all pass and documentation builds.

Single backticks typeset math so they render regular words in italic. All of the code terms in giac.py should be inside double backticks. Similarly, in calculus.py the word "cond" should be inside double backticks and one instance of "F" needs single backticks.

The documentation added to giac.py doesn't appear in the built documents because _sage_ is a private method. Is this a deliberate choice? How does one access the documentation for this method?

thanks for the feedback. for the proper string/latex formatting that's good to fix. the problem i have with ReST is that it is absurdly picky when whitespace messes up the proper format; i think it's too much for 2017 but that's way the way it is i guess..

anyway, in the new commits, both remarks have been addressed, taking into account the suggestion by rws, and improve a bit the tutorial.

@mforets
Copy link
Mannequin Author

mforets mannequin commented Apr 6, 2017

comment:40

Replying to @rwst:

As far as I know there is no way to display such docstrings (I have struggled myself with this in the past). If something important is there it should be moved to the global documentation at the top of the file.

sorry for the noise but is there a ticket for this already? just discovered that the algorithm description of solve_linear_de is great, but it is hidden (!)

@rwst
Copy link

rwst commented Apr 8, 2017

Reviewer: Paul Masson, Ralf Stephan

@rwst
Copy link

rwst commented Apr 8, 2017

comment:41

Patchbot has one error but I cannot confirm it, so we're good. Paul, I dared to add your name to reviewers---just change if you don't want that.

@vbraun
Copy link
Member

vbraun commented Apr 10, 2017

Changed branch from u/mforets/laplace_transform_involving_time_shifts to f845269

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants