# Laplace

Notebook to perform Laplace transformations

Author: Lucas Schneider

---

## Initialization


```python
import sympy as sym
import numpy as np

from sympy.integrals import laplace_transform
from sympy.integrals import inverse_laplace_transform

from IPython.display import display
from IPython.display import Math
from sympy.interactive import printing

import scipy.signal as sig

# sym.init_printing()
```


```python
t_var = sym.symbols('t', real=True)
s_var = sym.symbols('s')
d_domain = 'RR'
d_digits = 5
```


```python
def round_expr(expr, num_digits = d_digits):
    return expr.xreplace({n : round(n, num_digits) for n in expr.atoms(sym.Number)})
```

---

## Direct Transform


```python
# Other symbols declaration
phi = sym.symbols('φ', real=True)
```


```python
# Input
expression = sym.cos(t_var)

expression *= sym.Heaviside(t_var)
print(str(expression).replace('**','^'))
expression
```


```python
U = laplace_transform(expression, t_var, s_var)
U[0]
```

---

## Inverse Transform

### Functions definition


```python
def poly_from_list(num_list, den_list, domain = d_domain, print_latex=True, print_string=False):
    '''
    Transforms the list of coefficients into sympy polynomial objects

    Parameters
    ----------
    num_list : list of floats
    den_list : list of floats
    domain : string
    print_latex : bool
    print_string : bool

    Returns
    -------
    num_pol : sympy.Poly
    den_pol : sympy.Poly
    
    '''

    num_pol = sym.Poly(num_list, s_var, domain=domain)
    den_pol = sym.Poly(den_list, s_var, domain=domain)

    F = num_pol / den_pol
    F_str = str(F).replace('**','^').replace('I', 'i').replace('exp', 'e^')
    
    if print_latex:
        display(Math(f'F(s) = {printing.default_latex(F)}'))

    if print_string:
        print(f'String for Wolfram: {F_str}\n')

    return num_pol, den_pol
```


```python
def residue(num_pol, den_pol, print_latex=True, print_string=False):
    '''
    Calulates the partial fraction and residues for the given rational polynomial

    Parameters
    ----------
    num_pol : sympy.Poly
    den_pol : sympy.Poly
    print_latex : bool
    print_string : bool

    Returns
    -------
    residues: list of floats
        Resiudes of the fraction part
    poles: list of floats
        Poles of the fraction part
    multi: list of floats
        Poles's multiplicity
    complete : list of floats
        Coefficients of complete part
    '''

    r, p, k = sig.residue(N_list, D_list)

    m=[]
    last_pole = None
    last_s = 1
    for pole in p:
        if pole == last_pole:
            m.append(last_multi + 1)
            last_multi += 1
        else:
            m.append(1)
            last_multi = 1
            last_pole = pole

    # TODO: Print residues table (with Pandas?)
    
    return r, p, m, k

```


```python
def partial_fractions(residues, poles, multi, complete, domain=d_domain):
    '''
    Display the partial fraction representation from the residues

    Parameters
    ----------
    residues: list of floats
    poles: list of floats
    multi: list of floats
    complete : list of floats
    domain : string
    '''

    terms = []
    terms.append(sym.Poly(complete, s_var, domain=domain).as_expr())
    for i in range(len(residues)):
        residue = residues[i]
        pole = poles[i]
        mult = multi[i]

        term = sym.Mul(residue, 1/(s_var - pole)**mult)
        terms.append(term)

    display(Math(f'F(t) = {printing.default_latex(round_expr(sum(terms)))}'))

```


```python
def conjugate_seen(pole, multi, poles_seen):
    '''
    Checks if an specific pole or its conjugate has already been computed before, with the same multiplicity

    Parameters
    ----------
    pole : copmlex
    multi : int
    poles_seen : list of tuples, with pole and its multiplicity

    Returns
    -------
    True or False
    '''
    for pole_seen, multi_seen in poles_seen:
        if (pole == pole_seen.conjugate()) and (multi == multi_seen):
            return True

    return False


def ILP_from_residues(residues, poles, multi, complete, complex_simplify=True, print_latex=True, print_string=False):
    '''
    Perform the Inverse Laplace Transform for the given residues

    Parameters
    ----------
    residues: list of floats
    poles: list of floats
    multi: list of floats
    complete : list of floats

    complex_simplify : bool
    print_latex : bool
    print_string : bool

    Returns
    -------
    f : sympy.Expr
    '''

    f_terms_complete = [complete[i] * sym.diff(sym.DiracDelta(t_var), t_var, len(complete) -1 -i) for i in range(len(complete))]

    f_terms_fraction = []

    complex_poles_pairs = []
    for i in range(len(residues)):
        residue = complex(residues[i])
        pole = complex(poles[i])
        mult = multi[i]

        if (abs(pole.imag) > 0.) and complex_simplify:
            term = 0
            if not conjugate_seen(pole, mult, complex_poles_pairs):
                if pole.imag < 0.:
                    Ak = residue.conjugate()
                    pk = pole.conjugate()
                else:
                    Ak = residue
                    pk = pole

                term = sym.Mul(t_var ** (mult-1) / sym.factorial(mult-1), 2 * sym.sqrt(Ak.real ** 2 + Ak.imag ** 2), sym.exp(pk.real * t_var), sym.cos(sym.Add(pk.imag*t_var, sym.arg(Ak)), evaluate=False), evaluate=False)

                complex_poles_pairs.append((pole, mult))
            else:
                continue
        else:
            term = sym.Mul(residues[i] / sym.factorial(mult - 1), t_var ** (mult - 1), sym.exp(poles[i] * t_var), evaluate=False)
 
        f_terms_fraction.append(term)

    f = sum(f_terms_complete) + sym.Mul(sum(f_terms_fraction), sym.Heaviside(t_var), evaluate=False)
    f = round_expr(f)
    f_str = str(f).replace('**','^').replace('I', 'i').replace('exp', 'e^')

    if print_latex:
        display(Math(f'f(t) = {printing.default_latex(f)}'))

    if print_string:
        print(f'String for Wolfram: {f_str}\n')

    return f

```


```python
def evaluate_f(f, t, print_latex = True):
    f_t = f.evalf(subs={t_var: t})

    if print_latex:
        display(Math(f'f({t}) = {f_t}'))
        
    return f_t

```

### Rational functions

Define the polynomials


```python
N_list = [4,0]
D_list = [1,2,16,32]

polys = poly_from_list(N_list, D_list, print_string=True)
```


$\displaystyle F(s) = \frac{4.0 s}{1.0 s^{3} + 2.0 s^{2} + 16.0 s + 32.0}$


    String for Wolfram: 4.0*s/(1.0*s^3 + 2.0*s^2 + 16.0*s + 32.0)
    



```python
res = residue(N_list, D_list)
partial_fractions(*res)
f_res = ILP_from_residues(*res, complex_simplify=True, print_string=False)
```


$\displaystyle F(t) = \frac{0.2 + 0.4 i}{s + 4.0 i} + \frac{0.2 - 0.4 i}{s - 4.0 i} - \frac{0.4}{s + 2.0}$



$\displaystyle f(t) = \left(0.89443 \cos{\left(4.0 t - 1.10715 \right)} - 0.4 e^{- 2.0 t}\right) \theta\left(t\right)$


### Any function


```python
# Use the polynomials from above

polys = poly_from_list(N_list, D_list, domain='QQ', print_latex=False)
F_func = polys[0] / polys[1]

# Define function with sympy expression

# F_func = (sym.cos(sym.pi/4) * s_var -2*sym.sin(sym.pi/4))/(s_var**2 + 4)

display(Math(f'F(t) = {printing.default_latex(F_func)}'))
```


$\displaystyle F(t) = \frac{s + 3}{s^{3} + 3 s^{2} + 2 s}$



```python
f_func = inverse_laplace_transform(F_func, s_var, t_var).simplify()
display(Math(f'f(t) = {printing.default_latex(f_func)}'))

# result = evaluate_f(f_res, 0.5)
```


$\displaystyle f(t) = \frac{\left(3 e^{2 t} - 4 e^{t} + 1\right) e^{- 2 t} \theta\left(t\right)}{2}$



```python

```

33622bbefe3803a04a12beb48e3d17ac2b62488b
ipynb
Jupyter Notebook
laplace.ipynb
lucastrschneider/laplace_sym
8f019fa1f05de5da6e0dae9a672496b227d30efc
MIT
2021-02-14T14:32:10.000Z
2021-02-14T14:32:10.000Z
laplace.ipynb
lucastrschneider/laplace_sym
8f019fa1f05de5da6e0dae9a672496b227d30efc
MIT
laplace.ipynb
lucastrschneider/laplace_sym
8f019fa1f05de5da6e0dae9a672496b227d30efc
MIT
2020-09-10T03:15:27.000Z
2020-09-10T03:15:27.000Z
Qwen/Qwen-72B
1. YES
2. YES