Laplace Transform

docs/make.jl

@@ -53,6 +53,7 @@ makedocs(
             "Algebra" => [
                 "manual/solver.md",
                 "manual/groebner.md",
+                "manual/partial_fractions.md",
             ],
 
             "Calculus" => [

docs/src/manual/ode.md

@@ -43,6 +43,16 @@ Symbolics.solve_symbolic_IVP
 Symbolics.solve_linear_ode_system
 ```
 
+### Laplace Transform
+
+The Laplace transform can be used to solve ODEs by transforming the whole equation, solving algebraically, then applying the inverse transform. The Laplace transform and inverse transform functionality is currently based on a rule table and applying linearity, so this method is limited in what expressions are able to be transformed and inverse transformed.
+
+```@docs
+Symbolics.laplace
+Symbolics.inverse_laplace
+Symbolics.laplace_solve_ode
+```
+
 ### SymPy 
 
 ```@docs

docs/src/manual/partial_fractions.md

@@ -0,0 +1,9 @@
+# Partial Fraction Decomposition
+
+Partial fraction decomposition is performed using the cover-up method. This involves "covering up" a factor in the denominator and substituting the root into the remaining expression. When the denominator can be completely factored into non-repeated linear factors, this produces the desired result. When there are repeated or irreducible quadratic factors, it produces terms with unknown coefficients in the numerator that is solved as a system of equations.
+
+It is often used when solving integrals or performing an inverse Laplace transform (see [`inverse_laplace`](@ref)).
+
+```docs
+Symbolics.partial_frac_decomposition
+```

src/Symbolics.jl

@@ -175,6 +175,9 @@ include("operators.jl")
 include("limits.jl")
 export limit
 
+include("partialfractions.jl")
+export partial_frac_decomposition
+
 # Hacks to make wrappers "nicer"
 const NumberTypes = Union{AbstractFloat,Integer,Complex{<:AbstractFloat},Complex{<:Integer}}
 (::Type{T})(x::SymbolicUtils.Symbolic) where {T<:NumberTypes} = throw(ArgumentError("Cannot convert Sym to $T since Sym is symbolic and $T is concrete. Use `substitute` to replace the symbolic unwraps."))
@@ -223,8 +226,9 @@ export symbolic_solve
 # Diff Eq Solver
 include("diffeqs/diffeqs.jl")
 include("diffeqs/systems.jl")
+include("diffeqs/laplace.jl")
 include("diffeqs/diffeq_helpers.jl")
-export SymbolicLinearODE, symbolic_solve_ode, solve_linear_ode_system, solve_symbolic_IVP
+export SymbolicLinearODE, symbolic_solve_ode, solve_linear_ode_system, solve_symbolic_IVP, laplace, inverse_laplace, laplace_solve_ode
 
 # Sympy Functions
 

src/diffeqs/diffeq_helpers.jl

@@ -12,25 +12,37 @@ function _get_der_order(expr, x, t)
         return maximum(_get_der_order.(factors(expr), Ref(x), Ref(t)))
     end
 
-    return _get_der_order(substitute(expr, Dict(Differential(t)(x) => x)), x, t) + 1
+    return _get_der_order(fast_substitute(expr, Dict(Differential(t)(x) => x)), x, t) + 1
+end
+
+function reduce_rule(expr, Dt)
+    iscall(expr) && isequal(operation(expr), Dt) ? wrap(arguments(expr)[1]) : nothing
+end
+
+"""
+    unwrap_der(expr, Dt)
+
+Helper function to unwrap derivatives of `f(t)` in `expr` with respect to the differential operator `Dt = Differential(t)`. Returns a tuple `(n, base_expr)`, where `n` is the order of the derivative and `base_expr` is the expression with the derivatives removed. If `expr` does not contain `f(t)` or its derivatives, returns `(0, expr)`.
+"""
+function unwrap_der(expr, Dt)
+
+    if reduce_rule(unwrap(expr), Dt) === nothing
+        return 0, expr
+    end
+
+    order, expr = unwrap_der(reduce_rule(unwrap(expr), Dt), Dt)
+    return order + 1, expr
 end
 
 # takes into account fractions
 function _true_factors(expr)
-    facs = factors(expr)
-    true_facs::Vector{Number} = []
-    frac_rule = @rule (~x)/(~y) => [~x, 1/~y]
-    for fac in facs
-        frac = frac_rule(fac)
-        if frac !== nothing && !isequal(frac[1], 1)
-            append!(true_facs, _true_factors(frac[1]))
-            append!(true_facs, _true_factors(frac[2]))
-        else
-            push!(true_facs, fac)
-        end
-    end
+    expr = flatten_fractions(unwrap(expr)) # flatten nested fractions
+
+    numerator_factors = SymbolicUtils.numerators(unwrap(expr))
+    denominator_factors = SymbolicUtils.denominators(unwrap(expr))
 
-    return convert(Vector{Num}, true_facs)
+    facs = filter(fac -> !isequal(fac, 1), [numerator_factors; 1 ./ denominator_factors])
+    return isempty(facs) ? [1] : facs
 end
 
 """
@@ -45,7 +57,7 @@ function reduce_order(eq, x, t, ys)
 
     # reduction of order
     y_sub = Dict([[(Dt^i)(x) => ys[i+1] for i=0:n-1]; (Dt^n)(x) => variable(:𝒴)])
-    eq = substitute(eq, y_sub)
+    eq = fast_substitute(eq, y_sub)
 
     # isolate (Dt^n)(x)
     f = symbolic_linear_solve(eq, variable(:𝒴), check=false)
@@ -60,7 +72,7 @@ function unreduce_order(expr, x, t, ys)
     Dt = Differential(t)
     rev_y_sub = Dict(ys[i] => (Dt^(i-1))(x) for i in eachindex(ys))
 
-    return substitute(expr, rev_y_sub)
+    return fast_substitute(expr, rev_y_sub)
 end
 
 function is_solution(solution, eq::Equation, x, t)
@@ -82,15 +94,12 @@ function is_solution(solution, eq, x, t)
 end
 
 function _parse_trig(expr, t)
-    parse_sin = Symbolics.Chain([(@rule sin(t) => 1), (@rule sin(~x * t) => ~x)])
-    parse_cos = Symbolics.Chain([(@rule cos(t) => 1), (@rule cos(~x * t) => ~x)])
-
-    if !isequal(parse_sin(expr), expr)
-        return parse_sin(expr), true
+    if iscall(expr) && isequal(operation(expr), sin) && any(isequal.(t, factors(arguments(expr)[1])))
+        return arguments(expr)[1]/t, true
     end
 
-    if !isequal(parse_cos(expr), expr)
-        return parse_cos(expr), false
+    if iscall(expr) && isequal(operation(expr), cos) && any(isequal.(t, factors(arguments(expr)[1])))
+        return arguments(expr)[1]/t, false
     end
 
     return nothing

src/diffeqs/diffeqs.jl

@@ -59,7 +59,7 @@ function is_linear_ode(expr, x, t)
     @assert n >= 1 "ODE must have at least one derivative"
 
     y_sub = Dict([[(Dt^i)(x) => ys[i+1] for i=0:n-1]; (Dt^n)(x) => variable(:𝒴)])
-    expr = substitute(expr, y_sub)
+    expr = fast_substitute(expr, y_sub)
 
     # isolate (Dt^n)(x)
     f = symbolic_linear_solve(expr, variable(:𝒴), check=false)
@@ -362,7 +362,7 @@ function get_rrf_coeff(q, t)
         return a, r_re + r_im * im
     end
 
-    a = prod(filter(fac -> isempty(Symbolics.get_variables(fac, [t])), facs))
+    a = prod([1; filter(fac -> isempty(Symbolics.get_variables(fac, [t])), facs)])
 
     not_a = filter(fac -> !isempty(Symbolics.get_variables(fac, [t])), facs) # should just be e^(rt)
     if length(not_a) != 1
@@ -384,7 +384,7 @@ For finding particular solution when q(t) = a*e^(rt)*cos(bt) (or sin(bt))
 function exp_trig_particular_solution(eq::SymbolicLinearODE)
     facs = _true_factors(eq.q)
 
-    a = prod(filter(fac -> isempty(Symbolics.get_variables(fac, [eq.t])), facs))
+    a = prod([1; filter(fac -> isempty(Symbolics.get_variables(fac, [eq.t])), facs)])
 
     not_a = filter(fac -> !isempty(Symbolics.get_variables(fac, [eq.t])), facs)
 
@@ -415,12 +415,12 @@ function exp_trig_particular_solution(eq::SymbolicLinearODE)
     @variables 𝓈
     p = characteristic_polynomial(eq, 𝓈)
     Ds = Differential(𝓈)
-    while isequal(substitute(expand_derivatives((Ds^k)(p)), Dict(𝓈 => r+b*im)), 0)
+    while isequal(fast_substitute(expand_derivatives((Ds^k)(p)), Dict(𝓈 => r+b*im)), 0)
         k += 1
     end
 
     rrf = expand(simplify(a * exp((r + b * im) * eq.t) * eq.t^k /
-                           (substitute(expand_derivatives((Ds^k)(p)), Dict(𝓈 => r+b*im)))))
+                           (fast_substitute(expand_derivatives((Ds^k)(p)), Dict(𝓈 => r+b*im)))))
 
     return is_sin ? imag(rrf) : real(rrf)
 end
@@ -447,12 +447,12 @@ function resonant_response_formula(eq::SymbolicLinearODE)
     @variables 𝓈
     p = characteristic_polynomial(eq, 𝓈)
     Ds = Differential(𝓈)
-    while isequal(substitute(expand_derivatives((Ds^k)(p)), Dict(𝓈 => r)), 0)
+    while isequal(fast_substitute(expand_derivatives((Ds^k)(p)), Dict(𝓈 => r)), 0)
         k += 1
     end
 
     return expand(simplify(a * exp(r * eq.t) * eq.t^k /
-                           (substitute(expand_derivatives((Ds^k)(p)), Dict(𝓈 => r)))))
+                           (fast_substitute(expand_derivatives((Ds^k)(p)), Dict(𝓈 => r)))))
 end
 
 function method_of_undetermined_coefficients(eq::SymbolicLinearODE)
@@ -476,8 +476,8 @@ function method_of_undetermined_coefficients(eq::SymbolicLinearODE)
         coeff_solution = nothing
     end
 
-    if degree > 0 && coeff_solution !== nothing && !isempty(coeff_solution) && isequal(expand(substitute(eq_subbed, coeff_solution[1])), 0)
-        return substitute(form, coeff_solution[1])
+    if degree > 0 && coeff_solution !== nothing && !isempty(coeff_solution) && isequal(expand(fast_substitute(eq_subbed, coeff_solution[1])), 0)
+        return fast_substitute(form, coeff_solution[1])
     end
 
     # exponential
@@ -487,13 +487,13 @@ function method_of_undetermined_coefficients(eq::SymbolicLinearODE)
 
         r = coeff[2]
         form = a_form*exp(r*eq.t)
-        eq_subbed = substitute(get_expression(eq), Dict(eq.x => form))
+        eq_subbed = fast_substitute(get_expression(eq), Dict(eq.x => form))
         eq_subbed = expand_derivatives(eq_subbed)
         eq_subbed = simplify(expand((eq_subbed.lhs - eq_subbed.rhs) / exp(r*eq.t)))
         coeff_solution = solve_interms_ofvar(eq_subbed, eq.t)
 
         if coeff_solution !== nothing && !isempty(coeff_solution)
-            return substitute(form, coeff_solution[1])
+            return fast_substitute(form, coeff_solution[1])
         end
     end
 
@@ -505,9 +505,9 @@ function method_of_undetermined_coefficients(eq::SymbolicLinearODE)
     if parsed !== nothing
         ω = parsed[1]
         form = 𝒶*cos(ω*eq.t) + 𝒷*sin(ω*eq.t)
-        eq_subbed = substitute(get_expression(eq), Dict(eq.x => form))
+        eq_subbed = fast_substitute(get_expression(eq), Dict(eq.x => form))
         eq_subbed = expand_derivatives(eq_subbed)
-        eq_subbed = expand(substitute(eq_subbed.lhs - eq_subbed.rhs, Dict(cos(ω*eq.t)=>𝒸𝓈, sin(ω*eq.t)=>𝓈𝓃)))
+        eq_subbed = expand(fast_substitute(eq_subbed.lhs - eq_subbed.rhs, Dict(cos(ω*eq.t)=>𝒸𝓈, sin(ω*eq.t)=>𝓈𝓃)))
         cos_eq = simplify(sum(filter(term -> !isempty(Symbolics.get_variables(term, 𝒸𝓈)), terms(eq_subbed)))/𝒸𝓈)
         sin_eq = simplify(sum(filter(term -> !isempty(Symbolics.get_variables(term, 𝓈𝓃)), terms(eq_subbed)))/𝓈𝓃)
         if !isempty(Symbolics.get_variables(cos_eq, [eq.t,𝓈𝓃,𝒸𝓈])) || !isempty(Symbolics.get_variables(sin_eq, [eq.t,𝓈𝓃,𝒸𝓈]))
@@ -517,7 +517,7 @@ function method_of_undetermined_coefficients(eq::SymbolicLinearODE)
         end
 
         if coeff_solution !== nothing && !isempty(coeff_solution)
-            return substitute(form, coeff_solution[1])
+            return fast_substitute(form, coeff_solution[1])
         end
     end
 end
@@ -556,7 +556,7 @@ function solve_symbolic_IVP(ivp::IVP)
         push!(eqs, eq)
     end
 
-    return expand(simplify(substitute(general_solution, symbolic_solve(eqs, ivp.eq.C)[1])))
+    return expand(simplify(fast_substitute(general_solution, symbolic_solve(eqs, ivp.eq.C)[1])))
 end
 
 """
@@ -622,7 +622,7 @@ function solve_clairaut(expr, x, t)
     end
 
     C = Symbolics.variable(:C, 1) # constant of integration
-    f = substitute(f, Dict(Dt(x) => C))
+    f = fast_substitute(f, Dict(Dt(x) => C))
     if !isempty(Symbolics.get_variables(f, [x]))
         return nothing
     end
@@ -649,7 +649,7 @@ function linearize_bernoulli(expr, x, t, v)
     for term in terms
         if Symbolics.hasderiv(Symbolics.value(term))
             facs = _true_factors(term)
-            leading_coeff = prod(filter(fac -> !Symbolics.hasderiv(Symbolics.value(fac)), facs))
+            leading_coeff = prod([1; filter(fac -> !Symbolics.hasderiv(Symbolics.value(fac)), facs)])
             if !isequal(term//leading_coeff, Dt(x))
                 return nothing
             end
@@ -661,10 +661,10 @@ function linearize_bernoulli(expr, x, t, v)
             end
 
             if isequal(x_fac[1], x)
-                p = prod(filter(fac -> isempty(Symbolics.get_variables(fac, [x])), facs))
+                p = prod([1; filter(fac -> isempty(Symbolics.get_variables(fac, [x])), facs)])
             else
                 n = degree(x_fac[1])
-                q = -prod(filter(fac -> isempty(Symbolics.get_variables(fac, [x])), facs))
+                q = -prod([1; filter(fac -> isempty(Symbolics.get_variables(fac, [x])), facs)])
             end
         end
     end