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

@@ -24,11 +24,32 @@ The analytical solution can be investigated symbolically using `observed(sys)`.
 
 ## Symbolically Solving ODEs
 
-Currently there is no native symbolic ODE solver. Though there are bindings to SymPy
-
 !!! note
     This area is currently under heavy development. More solvers will be available in the near future.
 
+```@docs
+Symbolics.LinearODE
+```
+
+```@docs
+Symbolics.symbolic_solve_ode
+```
+
+### Continuous Dynamical Systems
+```@docs
+Symbolics.solve_linear_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."))
@@ -220,6 +223,13 @@ include("solver/main.jl")
 include("solver/special_cases.jl")
 export symbolic_solve
 
+# Diff Eq Solver
+include("diffeqs/diffeqs.jl")
+include("diffeqs/systems.jl")
+include("diffeqs/laplace.jl")
+include("diffeqs/diffeq_helpers.jl")
+export LinearODE, IVP, symbolic_solve_ode, solve_linear_system, solve_IVP, laplace, inverse_laplace, laplace_solve_ode
+
 # Sympy Functions
 
 """

src/diffeqs/diffeq_helpers.jl

@@ -0,0 +1,113 @@
+# recursively find highest derivative order in `expr`
+function _get_der_order(expr, x, t)
+    if !hasderiv(unwrap(expr))
+        return 0
+    end
+
+    if length(terms(expr)) > 1
+        return maximum(_get_der_order.(terms(expr), Ref(x), Ref(t)))
+    end
+
+    if length(factors(expr)) > 1
+        return maximum(_get_der_order.(factors(expr), Ref(x), Ref(t)))
+    end
+
+    return _get_der_order(fast_substitute(expr, Dict(Differential(t)(x) => x)), x, t) + 1
+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)
+    reduce_rule = @rule Dt(~x) => ~x
+
+    if reduce_rule(expr) === nothing
+        return 0, expr
+    end
+
+    order, expr = unwrap_der(reduce_rule(expr), Dt)
+    return order + 1, expr
+end
+
+# takes into account fractions
+function _true_factors(expr)
+    facs = factors(expr)
+    true_facs::Vector{Real} = []
+    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
+
+    return convert(Vector{Num}, true_facs)
+end
+
+"""
+    reduce_order(eq, x, t, ys)
+
+Reduce order of an ODE by substituting variables for derivatives to form a system of first order ODEs
+"""
+function reduce_order(eq, x, t, ys)
+    Dt = Differential(t)
+    n = _get_der_order(eq, x, t)
+    @assert n >= 1 "ODE must have at least one derivative"
+
+    # reduction of order
+    y_sub = Dict([[(Dt^i)(x) => ys[i+1] for i=0:n-1]; (Dt^n)(x) => variable(:𝒴)])
+    eq = fast_substitute(eq, y_sub)
+
+    # isolate (Dt^n)(x)
+    f = symbolic_linear_solve(eq, variable(:𝒴), check=false)
+    @assert f !== nothing "Failed to isolate highest order derivative term"
+    f = f[1]
+    system = [ys[2:n]; f]
+
+    return system
+end
+
+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 fast_substitute(expr, rev_y_sub)
+end
+
+function is_solution(solution, eq::Equation, x, t)
+    is_solution(solution, eq.lhs - eq.rhs, x, t)
+end
+
+function is_solution(solution, eq::LinearODE)
+    is_solution(solution, get_expression(eq), eq.x, eq.t)
+end
+
+function is_solution(solution, eq, x, t)
+    if solution === nothing
+        return false
+    end
+
+    expr = substitute(eq, Dict(x => solution))
+    expr = expand(expand_derivatives(expr))
+    return isequal(expr, 0)
+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
+    end
+
+    if !isequal(parse_cos(expr), expr)
+        return parse_cos(expr), false
+    end
+
+    return nothing
+end