diff --git a/docs/make.jl b/docs/make.jl index caf010d0b..58ea83467 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -53,6 +53,7 @@ makedocs( "Algebra" => [ "manual/solver.md", "manual/groebner.md", + "manual/partial_fractions.md", ], "Calculus" => [ diff --git a/docs/src/manual/ode.md b/docs/src/manual/ode.md index 33e1ba1a1..3a6f6b03b 100644 --- a/docs/src/manual/ode.md +++ b/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 diff --git a/docs/src/manual/partial_fractions.md b/docs/src/manual/partial_fractions.md new file mode 100644 index 000000000..a7a95581b --- /dev/null +++ b/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 +``` \ No newline at end of file diff --git a/src/Symbolics.jl b/src/Symbolics.jl index 077360e8d..cbd02b148 100644 --- a/src/Symbolics.jl +++ b/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 """ diff --git a/src/diffeqs/diffeq_helpers.jl b/src/diffeqs/diffeq_helpers.jl new file mode 100644 index 000000000..5d175609a --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl new file mode 100644 index 000000000..4c0674ad8 --- /dev/null +++ b/src/diffeqs/diffeqs.jl @@ -0,0 +1,667 @@ +""" +Represents a linear ordinary differential equation of the form: + +dⁿx/dtⁿ + pβ‚™(t)(dⁿ⁻¹x/dtⁿ⁻¹) + ... + pβ‚‚(t)(dx/dt) + p₁(t)x = q(t) + +# Fields +- `x`: dependent variable +- `t`: independent variable +- `p`: coefficient functions of `t` ordered in increasing order (p₁, pβ‚‚, ...) +- `q`: right hand side function of `t`, without any `x` + +# Examples +```jldoctest +julia> using Symbolics + +julia> @variables x, t +2-element Vector{Num}: + x + t + +julia> eq = LinearODE(x, t, [1, 2, 3], 3exp(4t)) +(Dt^3)x + (3)(Dt^2)x + (2)(Dt^1)x + (1)(Dt^0)x ~ 3exp(4t) +``` +""" +struct LinearODE + x::Num + t::Num + p::AbstractArray + q::Any + C::Vector{Num} + + LinearODE(x, t, p, q) = new(x, t, p, q, variables(:C, 1:length(p))) + + function LinearODE(expr, x, t) + if expr isa Equation + expr = expr.lhs - expr.rhs + end + + expr = expand(simplify(expr)) + + @assert is_linear_ode(expr, x, t) "Equation must be linear in $x and $t" + + n = _get_der_order(expr, x, t) + + ys = variables(:π“Ž, 1:n) + A, b, islinear = linear_expansion(reduce_order(expr, x, t, ys), ys) + + p = expand.(simplify.(-A[end, 1:end])) + q = expand(simplify(b[end])) + + new(x, t, p, q) + end +end + +function is_linear_ode(expr, x, t) + Dt = Differential(t) + ys = variables(:π“Ž, 1:_get_der_order(expr, x, t)) + n = _get_der_order(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 = fast_substitute(expr, y_sub) + + # isolate (Dt^n)(x) + f = symbolic_linear_solve(expr, variable(:𝒴), check=false) + + # couldn't isolate + if f === nothing + return false + end + + f = f[1] + system = [ys[2:n]; f] + + A, b, islinear = linear_expansion(system, ys) + return islinear && all(isempty.(get_variables.(A, x))) +end + +Dt(eq::LinearODE) = Differential(eq.t) +order(eq::LinearODE) = length(eq.p) + +"""Generates symbolic expression to represent `LinearODE`""" +function get_expression(eq::LinearODE) + (Dt(eq)^order(eq))(eq.x) + sum([(eq.p[n]) * (Dt(eq)^(n - 1))(eq.x) for n in 1:length(eq.p)]) ~ eq.q +end + +function Base.string(eq::LinearODE) + "(D$(eq.t)^$(order(eq)))$(eq.x) + " * + join( + ["($(eq.p[length(eq.p)-n]))(D$(eq.t)^$(length(eq.p)-n-1))$(eq.x)" + for n in 0:(order(eq) - 1)], + " + ") * " ~ $(eq.q)" +end + +Base.print(io::IO, eq::LinearODE) = print(io, string(eq)) +Base.show(io::IO, eq::LinearODE) = print(io, eq) +Base.isequal(eq1::LinearODE, eq2::LinearODE) = + isequal(eq1.x, eq2.x) && isequal(eq1.t, eq2.t) && + isequal(eq1.p, eq2.p) && isequal(eq1.q, eq2.q) + +"""Returns true if q(t) = 0 for linear ODE `eq`""" +is_homogeneous(eq::LinearODE) = isempty(Symbolics.get_variables(eq.q)) +"""Returns true if all coefficient functions p(t) of `eq` are constant""" +has_const_coeffs(eq::LinearODE) = all(isempty.(Symbolics.get_variables.(eq.p))) +"""Returns homgeneous version of `eq` where q(t) = 0""" +to_homogeneous(eq::LinearODE) = LinearODE(eq.x, eq.t, eq.p, 0) + +""" +Returns the characteristic polynomial p of `eq` (must have constant coefficients) in terms of variable `r` + +p(D) = Dⁿ + aₙ₋₁Dⁿ⁻¹ + ... + a₁D + aβ‚€I +""" +function characteristic_polynomial(eq::LinearODE, r) + poly = 0 + @assert has_const_coeffs(eq) "ODE must have constant coefficients to generate characteristic polynomial" + p = [eq.p; 1] # add implied coefficient of 1 to highest order + for i in eachindex(p) + poly += p[i] * r^(i - 1) + end + + return poly +end + +""" + symbolic_solve_ode(eq::LinearODE) +Symbolically solve a linear ordinary differential equation + +# Arguments +- eq: a `LinearODE` to solve + +# Returns +Symbolic solution to the ODE + +# Supported Methods +- first-order integrating factor +- constant coefficient homogeneous solutions (can handle repeated and complex characteristic roots) +- exponential and resonant response formula particular solutions (for any linear combination of `exp`, `sin`, `cos`, or `exp` times `sin` or `cos` (e.g. `e^2t * cos(-t) + e^-3t + sin(5t))`) +- method of undetermined coefficients particular solutions +- linear combinations of above particular solutions + +# Examples + +```jldoctest +julia> using Symbolics; import Nemo, SymPy + +julia> @variables x, t +2-element Vector{Num}: + x + t + +# Integrating Factor (note that SymPy is required for integration) +julia> symbolic_solve_ode(LinearODE(x, t, [5/t], 7t)) +(C₁ + t^7) / (t^5) + +# Constant Coefficients and RRF (note that Nemo is required to find characteristic roots) +julia> symbolic_solve_ode(LinearODE(x, t, [9, -6], 4exp(3t))) +C₁*exp(3t) + Cβ‚‚*t*exp(3t) + (2//1)*(t^2)*exp(3t) + +julia> symbolic_solve_ode(LinearODE(x, t, [6, 5], 2exp(-t)*cos(t))) +C₁*exp(-2t) + Cβ‚‚*exp(-3t) + (1//5)*cos(t)*exp(-t) + (3//5)*exp(-t)*sin(t) + +# Method of Undetermined Coefficients +julia> symbolic_solve_ode(LinearODE(x, t, [-3, 2], 2t - 5)) +(11//9) - (2//3)*t + C₁*exp(t) + Cβ‚‚*exp(-3t) +``` +""" +function symbolic_solve_ode(eq::LinearODE) + homogeneous_solutions = find_homogeneous_solutions(eq) + + if is_homogeneous(eq) && homogeneous_solutions !== nothing + return homogeneous_solutions + end + + particular_solution = find_particular_solution(eq) + if homogeneous_solutions !== nothing && particular_solution !== nothing + return homogeneous_solutions + particular_solution + end + + if order(eq) == 1 + return integrating_factor_solve(eq) + end +end + +""" + symbolic_solve_ode(expr::Equation, x, t) +Symbolically solve an ODE + +# Arguments +- expr: a symbolic ODE +- x: dependent variable +- t: independent variable + +# Supported Methods +- all methods of solving linear ODEs mentioned for `symbolic_solve_ode(eq::LinearODE)` +- Clairaut's equation +- Bernoulli equations + +# Examples + +```jldoctest +julia> using Symbolics; import Nemo + +julia> @variables x, t +2-element Vector{Num}: + x + t + +julia> Dt = Differential(t) +Differential(t) + +# LinearODE (via constant coefficients and RRF) +julia> symbolic_solve_ode(9t*x - 6*Dt(x) ~ 4exp(3t), x, t) +C₁*exp(3t) + Cβ‚‚*t*exp(3t) + (2//1)*(t^2)*exp(3t) + +# Clairaut's equation +julia> symbolic_solve_ode(x ~ Dt(x)*t - ((Dt(x))^3), x, t) +C₁*t - (C₁^3) + +# Bernoulli equations +julia> symbolic_solve_ode(Dt(x) + (4//t)*x ~ t^3 * x^2, x, t) +1 / (C₁*(t^4) - (t^4)*log(t)) +``` +""" +function symbolic_solve_ode(expr::Equation, x, t) + clairaut = solve_clairaut(expr, x, t) + if clairaut !== nothing + return clairaut + end + + bernoulli = solve_bernoulli(expr, x, t) + if bernoulli !== nothing + return bernoulli + end + + if is_linear_ode(expr, x, t) + eq = LinearODE(expr, x, t) + return symbolic_solve_ode(eq) + end +end + +""" +Find homogeneous solutions of linear ODE `eq` with integration constants of `eq.C` + +Currently only works for constant coefficient ODEs +""" +function find_homogeneous_solutions(eq::LinearODE) + if has_const_coeffs(eq) + return const_coeff_solve(to_homogeneous(eq)) + end +end + +""" +Find a particular solution to linear ODE `eq` + +Currently works for any linear combination of exponentials, sin, cos, or an exponential times sin or cos (e.g. e^2t * cos(-t) + e^-3t + sin(5t)) +""" +function find_particular_solution(eq::LinearODE) + # if q has multiple terms, find a particular solution for each and sum together + terms = Symbolics.terms(eq.q) + if length(terms) != 1 + solutions = find_particular_solution.(LinearODE.(Ref(eq.x), Ref(eq.t), Ref(eq.p), terms)) + if any(s -> s === nothing, solutions) + return nothing + end + return sum(solutions) + end + + if has_const_coeffs(eq) + rrf = resonant_response_formula(eq) + if rrf !== nothing + return rrf + end + rrf_trig = exp_trig_particular_solution(eq) + if rrf_trig !== nothing + return rrf_trig + end + end + + undetermined_coeff = method_of_undetermined_coefficients(eq) + if undetermined_coeff !== nothing + return undetermined_coeff + end +end + +""" +Returns homogeneous solutions to linear ODE `eq` with constant coefficients + +xβ‚•(t) = C₁e^(r₁t) + Cβ‚‚e^(rβ‚‚t) + ... + Cβ‚™e^(rβ‚™t) +""" +function const_coeff_solve(eq::LinearODE) + @variables 𝓇 + p = characteristic_polynomial(eq, 𝓇) + roots = symbolic_solve(p, 𝓇, dropmultiplicity = false) + + # Handle complex + repeated roots + solutions = exp.(roots * eq.t) + for i in eachindex(solutions)[1:(end - 1)] + j = i + 1 + + if !isequal(imag(roots[i]), 0) && isequal(roots[i], conj(roots[j])) + solutions[i] = exp(real(roots[i] * eq.t)) * cos(imag(roots[i] * eq.t)) + solutions[j] = exp(real(roots[i] * eq.t)) * sin(imag(roots[i] * eq.t)) + end + + while j <= length(solutions) && isequal(roots[i], roots[j]) + solutions[j] *= eq.t # multiply by t for each repetition + j += 1 + end + end + + solution = sum(eq.C .* solutions) + if solution isa Complex && isequal(imag(solution), 0) + solution = real(solution) + end + + return solution +end + +""" +Solve almost any first order ODE using an integrating factor. Requires SymPy! +""" +function integrating_factor_solve(eq::LinearODE) + p = eq.p[1] # only p + v = 0 # integrating factor + if isempty(Symbolics.get_variables(p)) + v = exp(p * eq.t) + else + v = exp(sympy_integrate(p, eq.t)) + end + solution = (1 / v) * ((isequal(eq.q, 0) ? 0 : sympy_integrate(eq.q * v, eq.t)) + eq.C[1]) + + if !isempty(Symbolics.get_variables(solution, variable(:Integral))) + return nothing + end + return expand(Symbolics.sympy_simplify(solution)) +end + +""" +Returns a, r from q(t)=a*e^(rt) if it is of that form. If not, returns `nothing` +""" +function get_rrf_coeff(q, t) + facs = factors(q) + + # handle complex r + # very convoluted, could probably be improved (possibly by making heavier use of @rule) + + # Description of process: + # only one factor of c*e^((a + bi)t) -> c*cos(bt)e^at + i*c*sin(bt)e^(at) + # real(factor) / imag(factor) = cos(bt)/sin(bt) - can extract imaginary part b from this + # then, divide real(factor) = c*cos(bt)e^at by cos(bt) to get c*e^at + # call self to get c and a, then add back in b + get_b = Symbolics.Chain([ + (@rule cos(t) / sin(t) => 1), (@rule cos(~b * t) / sin(~b * t) => ~b)]) + if length(facs) == 1 && !isequal(imag(facs[1]), 0) && + !isequal(get_b(real(facs[1]) / imag(facs[1])), real(facs[1]) / imag(facs[1])) + r_im = get_b(real(facs[1]) / imag(facs[1])) + real_q = real(facs[1]) / cos(r_im * t) + if isempty(Symbolics.get_variables(real_q, [t])) + return real_q, r_im * im + end + a, r_re = get_rrf_coeff(real(facs[1]) / cos(r_im * t), t) + return a, r_re + r_im * im + end + + a = prod(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 + return nothing + end + + der = expand_derivatives(Differential(t)(not_a[1])) + r = simplify(der / not_a[1]) + if !isempty(Symbolics.get_variables(r, [t])) + return nothing + end + + return a, r +end + +""" +For finding particular solution when q(t) = a*e^(rt)*cos(bt) (or sin(bt)) +""" +function exp_trig_particular_solution(eq::LinearODE) + facs = _true_factors(eq.q) + + a = prod(filter(fac -> isempty(Symbolics.get_variables(fac, [eq.t])), facs)) + + not_a = filter(fac -> !isempty(Symbolics.get_variables(fac, [eq.t])), facs) + + r = nothing + b = nothing + is_sin = false + + if length(not_a) == 1 && _parse_trig(not_a[1], eq.t) !== nothing + r = 0 + b, is_sin = _parse_trig(not_a[1], eq.t) + elseif length(not_a) != 2 + return nothing + elseif get_rrf_coeff(not_a[1], eq.t) !== nothing && + _parse_trig(not_a[2], eq.t) !== nothing + r = get_rrf_coeff(not_a[1], eq.t)[2] + b, is_sin = _parse_trig(not_a[2], eq.t) + elseif get_rrf_coeff(not_a[2], eq.t) !== nothing && + _parse_trig(not_a[1], eq.t) !== nothing + r = get_rrf_coeff(not_a[2], eq.t)[2] + b, is_sin = _parse_trig(not_a[1], eq.t) + else + return nothing + end + + # do complex rrf + # figure out how many times p needs to be differentiated before denominator isn't 0 + k = 0 + @variables π“ˆ + p = characteristic_polynomial(eq, π“ˆ) + Ds = Differential(π“ˆ) + 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 / + (fast_substitute(expand_derivatives((Ds^k)(p)), Dict(π“ˆ => r+b*im))))) + + return is_sin ? imag(rrf) : real(rrf) +end + +""" +Returns a particular solution to a constant coefficient ODE with q(t) = a*e^(rt) + +Exponential Response Formula: x_p(t) = a*e^(rt)/p(r) where p(r) is characteristic polynomial + +Resonant Response Formula: If r is a characteristic root, multiply by t and take the derivative of p (possibly multiple times) +""" +function resonant_response_formula(eq::LinearODE) + @assert has_const_coeffs(eq) + + # get a and r from q = a*e^(rt) + rrf_coeff = get_rrf_coeff(eq.q, eq.t) + if rrf_coeff === nothing + return nothing + end + a, r = rrf_coeff + + # figure out how many times p needs to be differentiated before denominator isn't 0 + k = 0 + @variables π“ˆ + p = characteristic_polynomial(eq, π“ˆ) + Ds = Differential(π“ˆ) + 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 / + (fast_substitute(expand_derivatives((Ds^k)(p)), Dict(π“ˆ => r))))) +end + +function method_of_undetermined_coefficients(eq::LinearODE) + # constant + p = eq.p[1] + if isempty(Symbolics.get_variables(p, eq.t)) && isempty(Symbolics.get_variables(eq.q, eq.t)) + return eq.q // p + end + + # polynomial + degree = max(Symbolics.degree(eq.q, eq.t), Symbolics.degree.(eq.p, eq.t)...) # just a starting point + a = Symbolics.variables(:𝒢, 1:degree+1) + form = sum(a[n]*eq.t^(n-1) for n = 1:degree+1) + eq_subbed = substitute(get_expression(eq), Dict(eq.x => form)) + eq_subbed = eq_subbed.lhs - eq_subbed.rhs + eq_subbed = expand_derivatives(eq_subbed) + + try + coeff_solution = solve_interms_ofvar(eq_subbed, eq.t) + catch + coeff_solution = nothing + end + + 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 + coeff = get_rrf_coeff(eq.q, eq.t) + if coeff !== nothing + a_form = form # use form from polynomial case + + r = coeff[2] + form = a_form*exp(r*eq.t) + 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 fast_substitute(form, coeff_solution[1]) + end + end + + # sin and cos + # this is a hacky way of doing things + @variables 𝒢, 𝒷 + @variables π’Έπ“ˆ, π“ˆπ“ƒ + parsed = _parse_trig(_true_factors(eq.q)[end], eq.t) + if parsed !== nothing + Ο‰ = parsed[1] + form = 𝒢*cos(Ο‰*eq.t) + 𝒷*sin(Ο‰*eq.t) + eq_subbed = fast_substitute(get_expression(eq), Dict(eq.x => form)) + eq_subbed = expand_derivatives(eq_subbed) + 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,π“ˆπ“ƒ,π’Έπ“ˆ])) + coeff_solution = nothing + else + coeff_solution = symbolic_solve([cos_eq, sin_eq], [𝒢,𝒷]) + end + + if coeff_solution !== nothing && !isempty(coeff_solution) + return fast_substitute(form, coeff_solution[1]) + end + end +end + +""" +Initial value problem (IVP) for a linear ODE +""" +struct IVP + eq::LinearODE + initial_conditions::Vector{Num} # values at t = 0 of nth derivative of x + + function IVP(eq::LinearODE, initial_conditions::Vector{<:Number}) + @assert length(initial_conditions) == order(eq) "# of Initial conditions must match order of ODE" + new(eq, initial_conditions) + end +end + +function solve_IVP(ivp::IVP) + general_solution = symbolic_solve_ode(ivp.eq) + if general_solution === nothing + return nothing + end + + eqs = [] + for i in eachindex(ivp.initial_conditions) + eq::Num = expand_derivatives((Dt(ivp.eq)^(i-1))(general_solution)) - ivp.initial_conditions[i] + + eq = substitute(eq, Dict(ivp.eq.t => 0), fold=false) + + # make sure exp, sin, and cos don't evaluate to floats + exp0 = substitute(exp(ivp.eq.t), Dict(ivp.eq.t => 0), fold=false) + sin0 = substitute(sin(ivp.eq.t), Dict(ivp.eq.t => 0), fold=false) + cos0 = substitute(cos(ivp.eq.t), Dict(ivp.eq.t => 0), fold=false) + + eq = expand(simplify(substitute(eq, Dict(exp0 => 1, sin0 => 0, cos0 => 1), fold=false))) + push!(eqs, eq) + end + + return expand(simplify(fast_substitute(general_solution, symbolic_solve(eqs, ivp.eq.C)[1]))) +end + +""" +Solve Clairaut's equation of the form x = x'*t + f(x'). + +Returns solution of the form x = C*t + f(C) where C is a constant. +""" +function solve_clairaut(expr, x, t) + Dt = Differential(t) + rhs = 0 + if isequal(expr.rhs, x) + rhs = expr.lhs + elseif isequal(expr.lhs, x) + rhs = expr.rhs + else + return nothing + end + + terms = Symbolics.terms(rhs) + matched = false # if expr contains term Dt(x)*t + f = 0 + for term in terms + if isequal(term, Dt(x)*t) + matched = true + elseif !isempty(Symbolics.get_variables(term, [t])) + return nothing + else + f += term + end + end + + if !matched + return nothing + end + + C = Symbolics.variable(:C, 1) # constant of integration + f = fast_substitute(f, Dict(Dt(x) => C)) + if !isempty(Symbolics.get_variables(f, [x])) + return nothing + end + + return C*t + f +end + +""" +Linearize a Bernoulli equation of the form dx/dt + p(t)x = q(t)x^n into a `LinearODE` of the form dv/dt + (1-n)p(t)v = (1-n)q(t) where v = x^(1-n) +""" +function linearize_bernoulli(expr, x, t, v) + Dt = Differential(t) + + if expr isa Equation + expr = expr.lhs - expr.rhs + end + + terms = Symbolics.terms(expr) + + p = 0 + q = 0 + n = 0 + leading_coeff = 1 + 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)) + if !isequal(term//leading_coeff, Dt(x)) + return nothing + end + elseif !isempty(Symbolics.get_variables(term, [x])) + facs = _true_factors(term) + x_fac = filter(fac -> !isempty(Symbolics.get_variables(fac, [x])), facs) + if length(x_fac) != 1 + return nothing + end + + if isequal(x_fac[1], x) + p = prod(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)) + end + end + end + + p //= leading_coeff + q //= leading_coeff + + return LinearODE(v, t, [p*(1-n)], q*(1-n)), n +end + +""" +Solve Bernoulli equations of the form dx/dt + p(t)x = q(t)x^n. May require SymPy to solve using integrating factor +""" +function solve_bernoulli(expr, x, t) + @variables 𝓋 + linearized = linearize_bernoulli(expr, x, t, 𝓋) + if linearized === nothing + return nothing + end + + eq, n = linearized + + solution = symbolic_solve_ode(eq) + if solution === nothing + return nothing + end + + return simplify(solution^(1//(1-n))) +end \ No newline at end of file diff --git a/src/diffeqs/laplace.jl b/src/diffeqs/laplace.jl new file mode 100644 index 000000000..a158ab920 --- /dev/null +++ b/src/diffeqs/laplace.jl @@ -0,0 +1,296 @@ +import DomainSets.ClosedInterval + +# from https://tutorial.math.lamar.edu/Classes/DE/Laplace_Table.aspx +transform_rules(f, t, F, s) = Symbolics.Chain([ + @rule 1 => 1/s + @rule exp(t) => 1/(s - 1) + @rule exp(~a * t) => 1/(-~a + s) + @rule t => 1/s^2 + @rule t^~n => factorial(~n)/s^(~n + 1) + @rule sqrt(t) => term(sqrt, pi)/(2 * s^(3/2)) + @rule sin(t) => 1/(1 + s^2) + @rule sin(~a * t) => ~a/((~a)^2 + s^2) + @rule cos(t) => s/(1 + s^2) + @rule cos(~a * t) => s/((~a)^2 + s^2) + @rule t*sin(t) => 1/(1 + s^2)^2 + @rule t*sin(~a * t) => 2*~a*s / ((~a)^2 + s^2)^2 + @rule t*cos(t) => (s^2 - 1) / (1 + s^2)^2 + @rule t*cos(~a * t) => (-(~a)^2 + s^2) / ((~a)^2 + s^2)^2 + @rule sin(t) - t*cos(t) => 2 / (1 + s^2)^2 + @rule sin(~a*t) - ~a*t*cos(~a*t) => 2*(~a)^3 / ((~a)^2 + s^2)^2 + @rule sin(t) + t*cos(t) => 2s^2 / (1 + s^2)^2 + @rule sin(~a*t) + ~a*t*cos(~a*t) => 2*~a*s^2 / ((~a)^2 + s^2)^2 + @rule cos(~a*t) - ~a*t*sin(~a*t) => s*((~a)^2 + s^2) / ((~a)^2 + s^2)^2 + @rule cos(~a*t) + ~a*t*sin(~a*t) => s*(s^2 + 3*(~a)^2) / ((~a)^2 + s^2)^2 + @rule sin(~b + ~a*t) => (s*sin(~b) + ~a*cos(~b)) / ((~a)^2 + s^2) + @rule cos(~b + ~a*t) => (s*cos(~b) - ~a*sin(~b)) / ((~a)^2 + s^2) + @rule sinh(~a * t) => ~a/(-(~a)^2 + s^2) + @rule cosh(~a * t) => s/(-(~a)^2 + s^2) + @rule exp(~a*t) * sin(~b * t) => ~b / ((~b)^2 + (-~a+s)^2) + @rule exp(~a*t) * cos(~b * t) => (-~a+s) / ((~b)^2 + (-~a+s)^2) + @rule exp(~a*t) * sinh(~b * t) => ~b / (-(~b)^2 + (-~a+s)^2) + @rule exp(~a*t) * cosh(~b * t) => (-~a+s) / (-(~b)^2 + (-~a+s)^2) + @rule t^~n * exp(~a * t) => factorial(~n) / (-~a + s)^(~n + 1) + @rule t*exp(~a * t) => 1 / (-~a + s)^(2) + @rule t^~n * exp(t) => factorial(~n) / (s)^(~n + 1) + @rule t*exp(t) => 1 / (s)^(2) + @rule exp(~c*t) * ~g => laplace(~g, f, t, F, s - ~c) # s-shift rule + @rule f(t)*t => -Differential(s)(F(s)) # s-derivative rule + @rule f(t)*t^(~n) => (-1)^(~n) * (Differential(s)^~n)(F(s)) # s-derivative rule + @rule f(~a + t) => exp(~a*s)*F(s) # t-shift rule + @rule f(t) => F(s) +]) + +""" + laplace(expr, f, t, F, s) + +Performs the Laplace transform of `expr` with respect to the variable `t`, where `f(t)` is a function in `expr` being transformed, and `F(s)` is the Laplace transform of `f(t)`. Returns the transformed expression in terms of `s`. + +Note that `f(t)` and `F(s)` should be defined using `@syms` + +Currently relies mostly on linearity and a rules table. When the rules table does not apply, it falls back to the integral definition of the Laplace transform. + +# Examples + +```jldoctest +julia> @variables t, s +2-element Vector{Num}: + t + s + +julia> @syms f(t)::Real F(s)::Real +(f, F) + +julia> laplace(exp(4t) + 5, f, t, F, s) +5 / s + 1 / (-4 + s) + +julia> laplace(10 + 4t - t^2, f, t, F, s) +10 / s + 4 / (s^2) + -2 / (s^3) + +julia> laplace(exp(-2t)*cos(3t) + 5exp(-2t)*sin(3t), f, t, F, s) +(2 + s) / (9 + (2 + s)^2) + 15 / (9 + (2 + s)^2) + +julia> laplace(t^2 * f(t), f, t, F, s) # s-derivative rule +Differential(s)(Differential(s)(F(s))) + +julia> laplace(5f(t-4), f, t, F, s) # t-shift rule +5F(s)*exp(-4s) + +julia> laplace(log(t), f, t, F, s) # fallback to definition +Integral(t, 0.0 .. Inf)(exp(-s*t)*log(t)) +``` +""" +function laplace(expr, f, t, F, s) + expr = expand(expr) + + if isequal(expr, 0) + return 0 + end + + Dt = Differential(t) + + transformed = transform_rules(f, t, F, s)(expr) + + # Check if transformation was successful + if !isequal(transformed, expr) + return transformed + end + + # t-derivative rule ((Dt^n)(f(t)) -> s^n*F(s) - s^(n-1)*f(0) - s^(n-2)*f'(0) - ... - f^(n-1)(0)) + n, expr = unwrap_der(expr, Dt) + if n != 0 && isequal(expr, f(t)) + f0 = Symbolics.variables(:f0, 0:(n-1)) + transformed = s^n*F(s) + for i = 1:n + transformed -= s^(n-i)*f0[i] + end + + return transformed + end + + terms = Symbolics.terms(expr) + result = 0 + + # unable to apply linearity, so return based on definition + if length(terms) == 1 && length(filter(x->isempty(Symbolics.get_variables(x)), _true_factors(terms[1]))) == 0 + return Integral(t in ClosedInterval(0, Inf))(expr*exp(-s*t)) + end + + # apply linearity by splitting into terms and factoring out constants + for term in terms + factors = _true_factors(wrap(term)) + constant = filter(x -> isempty(Symbolics.get_variables(x)), factors) + if !isempty(constant) + result += laplace(term / constant[1], f, t, F, s) * constant[1] + else + result += laplace(term, f, t, F, s) + end + end + + return result +end + +function laplace(expr::Equation, f, t, F, s) + return laplace(expr.lhs, f, t, F, s) ~ laplace(expr.rhs, f, t, F, s) +end + +# postprocess_root prevents automatic evaluation of sqrt to its floating point value +function processed_sqrt(x) + return postprocess_root(term(sqrt, x)) +end + +# F and f aren't used here, but are here for future-proofing +inverse_transform_rules(F, s, f, t) = Symbolics.Chain([ + @rule 1/s => 1 + @rule 1/(~a + s) => exp(-~a * t) + @rule 1/s^(~n) => t^(~n-1) / factorial(~n-1) + @rule 1/(2 * s^(3/2)) => sqrt(t)/term(term(sqrt, pi)) + @rule 1/(~a + s^2) => sin(processed_sqrt(~a) * t)/processed_sqrt(~a) + @rule s/(~a + s^2) => cos(processed_sqrt(~a) * t) + @rule s / (~a + s^2)^2 => t*sin(processed_sqrt(~a) * t)/(2*processed_sqrt(~a)) + @rule (-~a + s^2) / (~a + s^2)^2 => t*cos(processed_sqrt(~a) * t) + @rule 1 / (~a + s^2)^2 => (sin(processed_sqrt(~a)*t) - processed_sqrt(~a)*t*cos(processed_sqrt(~a)*t))/ (2*processed_sqrt(~a)^3) + @rule s^2 / (~a + s^2)^2 => (sin(processed_sqrt(~a)*t) + processed_sqrt(~a)*t*cos(processed_sqrt(~a)*t)) / (2*processed_sqrt(~a)) + @rule s*(~a + s^2) / (~a + s^2)^2 => cos(processed_sqrt(~a)*t) - processed_sqrt(~a)*t*sin(processed_sqrt(~a)*t) + @rule s*(3*~a + s^2) / (~a + s^2)^2 => cos(processed_sqrt(~a)*t) + processed_sqrt(~a)*t*sin(processed_sqrt(~a)*t) + @rule (s*sin(~b) + ~a*cos(~b)) / (~a + s^2) => sin(~b + processed_sqrt(~a)*t) + @rule (s*cos(~b) - ~a*sin(~b)) / ((~a)^2 + s^2) => cos(~b + ~a*t) + @rule 1/(s^2 - (~b)^2) => sinh(~b * t)/~b + @rule s/(s^2 - (~b)^2) => cosh(~b * t) + @rule 1 / ((~c+s)^2 + (~b)^2) => exp(-~c*t) * sin(~b * t) / ~b + @rule (~c+s) / ((~c+s)^2 + (~b)^2) => exp(-~c*t) * cos(~b * t) + @rule 1 / ((~c+s)^2 - (~b)^2) => exp(-~c*t) * sinh(~b * t) / ~b + @rule (~c+s) / ((~c+s)^2 - (~b)^2) => exp(-~c*t) * cosh(~b * t) + @rule 1 / (~a + s)^(~n) => t^(~n-1) * exp(-~a * t) / factorial(~n-1) +]) + +""" + inverse_laplace(expr, F, s, f, t) + +Performs the inverse Laplace transform of `expr` with respect to the variable `s`, where `F(s)` is the Laplace transform of `f(t)`. Returns the transformed expression in terms of `t`. + +Note that `f(t)` and `F(s)` should be defined using `@syms`. + +Will perform partial fraction decomposition and linearity before applying the inverse Laplace transform rules. When unable to find a result, returns `nothing`. + +# Examples + +```jldoctest +julia> @variables t, s +2-element Vector{Num}: + t + s + +julia> @syms f(t)::Real F(s)::Real +(f, F) + +julia> inverse_laplace(7/(s+3)^3, F, s, f, t) +(7//2)*(t^2)*exp(-3t) + +julia> inverse_laplace((s+2)/(s^2 - 3s - 4), F, s, f, t) # using partial fraction decomposition +-(1//5)*exp(-t) + (6//5)*exp(4t) + +julia> inverse_laplace(1/s^4, F, s, f, t) +(1//6)*(t^3) +``` +""" +function inverse_laplace(expr, F, s, f, t) + if isequal(expr, 0) + return 0 + end + + # check for partial fractions + partial_fractions = partial_frac_decomposition(expr, s) + if partial_fractions !== nothing && !isequal(partial_fractions, expr) + return inverse_laplace(partial_fractions, F, s, f, t) + end + + transformed = inverse_transform_rules(F, s, f, t)(expr) + + # Check if transformation was successful + if !isequal(transformed, expr) + return transformed + end + + _terms = terms(numerator(expr)) ./ denominator(expr) + + result = 0 + if length(_terms) == 1 && length(filter(x -> isempty(get_variables(x)), _true_factors(_terms[1]))) == 0 + @warn "Inverse laplace failed: $expr" + return nothing # no result + end + + # apply linearity by splitting into terms and factoring out constants + for term in _terms + factors = _true_factors(term) + constant = filter(x -> isempty(Symbolics.get_variables(x)), factors) + if !isempty(constant) + result += inverse_laplace(term / constant[1], F, s, f, t) * constant[1] + else + result += inverse_laplace(term, F, s, f, t) + end + end + + return result +end + +function inverse_laplace(expr::Equation, F, s, f, t) + return inverse_laplace(expr.lhs, F, s, f, t) ~ inverse_laplace(expr.rhs, F, s, f, t) +end + +""" + laplace_solve_ode(eq, f, t, f0) + +Solves the ordinary differential equation `eq` for the function `f(t)` using the Laplace transform method. + +`f0` is a vector of initial conditions evaluated at `t=0` (`[f(0), f'(0), f''(0), ...]`, must be same length as order of `eq`). + +# Examples + +```jldoctest +@variables t, s +@syms f(t)::Real F(s)::Real + +Dt = Differential(t) + +julia> laplace_solve_ode(Dt(f(t)) + 3f(t) ~ t^2*exp(-3t) + t*exp(-2t) + t, f, t, [1]) +-(1//9) + (1//3)*t + (19//9)*exp(-3t) - exp(-2t) + t*exp(-2t) + (1//3)*(t^3)*exp(-3t) + +julia> laplace_solve_ode((Dt^2)(f(t)) + f(t) ~ 2 + 2cos(t), f, t, [0, 0]) +(2//1) - (2//1)*cos(t) + t*sin(t) + +julia> laplace_solve_ode((Dt^3)(f(t)) - Dt(f(t)) ~ 6 - 3t^2, f, t, [1, 1, 1]) +exp(t) + t^3 +``` +""" +function laplace_solve_ode(eq, f, t, f0) + s = variable(:π“ˆ) + @syms 𝓕(s) + + # transform equation + transformed_eq = laplace(eq, f, t, 𝓕, s) + # substitute in initial conditions + transformed_eq = fast_substitute(transformed_eq, Dict(𝓕(s) => variable(:𝓕), [variable(:f0, i-1) => f0[i] for i=1:length(f0)]...)) + transformed_eq = expand(transformed_eq.lhs - transformed_eq.rhs) + + # solve for/isolate F(s) + F_terms = 0 + other_terms = [] + for term in terms(transformed_eq) + if isempty(get_variables(term, [variable(:𝓕)])) + push!(other_terms, -1*term) + else + F_terms += term/variable(:𝓕) # assumes term is something times F + end + end + + if isempty(other_terms) + other_terms = 0 + end + + # (a + b + ...)*F(s) = (c + d + ...) -> F(s) = (c + d + ...) / (a + b + ...) + transformed_soln = simplify(sum(other_terms ./ F_terms)) + + # perform inverse laplace transform to get f(t) + return expand(inverse_laplace(transformed_soln, 𝓕, s, f, t)) +end \ No newline at end of file diff --git a/src/diffeqs/systems.jl b/src/diffeqs/systems.jl new file mode 100644 index 000000000..df0efe281 --- /dev/null +++ b/src/diffeqs/systems.jl @@ -0,0 +1,115 @@ +# import Symbolics: symbolic_solve + +""" +Returns evolution matrix e^(tD) +""" +evo_mat(D::Matrix{<:Number}, t::Num) = diagm(exp.(t .* diag(D))) + +""" + solve_linear_system(A::Matrix{<:Number}, x0::Vector{<:Number}, t::Num) +Solve linear continuous dynamical system of differential equations of the form Ax = x' with initial condition x0 + +# Arguments +- `A`: matrix of coefficients +- `x0`: initial conditions vector +- `t`: independent variable + +# Returns +vector of symbolic solutions + +# Examples +!!! note uses method `symbolic_solve`, so packages `Nemo` and `Groebner` are often required +```jldoctest +julia> @variables t +1-element Vector{Num}: + t + +julia> solve_linear_system([1 0; 0 -1], [1, -1], t) # requires Nemo +2-element Vector{Num}: + exp(t) + -exp(-t) + +julia> solve_linear_system([-3 4; -2 3], [7, 2], t) # requires Groebner +2-element Vector{Num}: + (10//1)*exp(-t) - (3//1)*exp(t) + (5//1)*exp(-t) - (3//1)*exp(t) +``` +""" +function solve_linear_system(A::Matrix{<:Number}, x0::Vector{<:Number}, t::Num) + # Check A is square + if size(A, 1) != size(A, 2) + throw(ArgumentError("Matrix A must be square.")) + end + + # Check x0 matches size of A + if size(A, 1) != length(x0) + throw(ArgumentError("Initial condition vector x0 must match the size of matrix A.")) + end + + if isdiag(A) + # If A is diagonal, use uncoupled system solver + return solve_uncoupled_system(A, x0, t) + end + + S, D = diagonalize(A) + + return simplify.(S * evo_mat(D, t) * rationalize.(S^-1) * x0) +end + +""" +Solve a system of uncoupled ODEs of the form: + + dx/dt = A*x + +where A is a diagonal matrix. +""" +function solve_uncoupled_system(A::Matrix{<:Number}, x0::Vector{<:Number}, t::Num) + # Check A is diagonal + if !isdiag(A) + throw(ArgumentError("Matrix A must be diagonal.")) + end + + return evo_mat(A, t) * x0 +end + +""" +Diagonalize matrix A, returning matrix S with eigenvectors as columns and diagonal matrix D with eigenvalues +""" +function diagonalize(A::Matrix{<:Number})::Tuple{Matrix{<:Number},Matrix{<:Number}} + eigs::Eigen = symbolic_eigen(A) + S = eigs.vectors + D = diagm(eigs.values) + return S, D +end + + +""" +Replacement for `LinearAlgebra.eigen` function that uses symbolic functions to avoid floating-point inaccuracies +""" +function symbolic_eigen(A::Matrix{<:Number}) + @variables Ξ» # eigenvalue + v = variables(:v, 1:size(A, 1)) # vector of subscripted variables to represent eigenvector + + # find eigenvalues first + p = det(Ξ»*I - A) ~ 0 # polynomial to solve + values = symbolic_solve(p, Ξ») # solve polynomial + + # then, find eigenvectors + S::Matrix{Number} = Matrix(I, size(A, 1), 0) # matrix storing vertical eigenvectors + + for value in values + eqs = (value*I - A) * v# .~ zeros(size(A, 1)) # equations to give eigenvectors + eqs = fast_substitute(eqs, Dict(v[1] => 1)) # set first element to 1 to constrain solution space + + sol = symbolic_solve(eqs[1:end-1], v[2:end]) # solve all but one equation (because of constraining solutions above) + + # parse multivar solutions into Vector (in order) + if sol[1] isa Dict + sol = [sol[1][var] for var in v[2:end]] + end + vec::Vector{Number} = prepend!(sol, [1]) # add back the 1 (representing v_1) from substitution + S = [S vec] # add vec to matrix + end + + return Eigen(values, S) +end \ No newline at end of file diff --git a/src/partialfractions.jl b/src/partialfractions.jl new file mode 100644 index 000000000..31ba00171 --- /dev/null +++ b/src/partialfractions.jl @@ -0,0 +1,175 @@ +# used to represent linear or irreducible quadratic factors +struct Factor + expr + root + multiplicity + x +end + +function Factor(expr, multiplicity, x) + fac_rule = @rule x + ~r => -~r + return Factor(expr, fac_rule(expr), multiplicity, x) +end + +function Base.isequal(a::Factor, b::Factor) + return isequal(a.expr, b.expr) && a.multiplicity == b.multiplicity +end + +# https://math.mit.edu/~hrm/18.031/pf-coverup.pdf +""" + partial_frac_decomposition(expr, x) + +Performs partial fraction decomposition for expressions with linear, repeated, or irreducible quadratic factors in the denominator. Can't currently handle irrational roots. + +When leading coefficient of the denominator is not 1, it will be factored out and then put back in at the end, often leading to non-integer coefficients in the result. Will return `nothing` if the expression is not a valid polynomial fraction, or if it has irrational roots. + +# Examples + +```jldoctest +julia> @variables x +1-element Vector{Num}: + x + +julia> partial_frac_decomposition((3x-1) / (x^2 + x - 6), x) +(1//1) / (-2 + x) + (2//1) / (3 + x) + +julia> partial_frac_decomposition((4x^3 + 16x + 7)/(x^2 + 4)^2, x) # repeated irreducible quadratic factor +(4x) / (4 + x^2) + 7 / ((4 + x^2)^2) + +julia> partial_frac_decomposition((4x^2 - 22x + 7)/((2x+3)*(x-2)^2), x) # non-one leading coefficient +(-3//1) / ((-2 + x)^2) + (2//1) / ((3//2) + x) +``` + +!!! note that irreducible quadratic and repeated linear factors require the `Groebner` package to solve a system of equations +""" +function partial_frac_decomposition(expr, x) + A, B = numerator(expr), denominator(expr) + + # check if both numerator and denominator are polynomials + if !isequal(polynomial_coeffs(A, [x])[2], 0) || !isequal(polynomial_coeffs(B, [x])[2], 0) + return nothing + end + + if degree(A) >= degree(B) + return nothing + end + + facs = factorize(B, x) + if facs === nothing + return nothing + end + + leading_coeff = coeff_vector(expand(B), x)[end] # of denominator + + # already in partial fraction form + if length(facs) == 1 && only(facs).multiplicity == 1 && degree(A) <= 1 + return expr + end + + result = [] + c_idx = 0 # index to keep track of which C subscript to use + if length(facs) == 1 + fac = only(facs) + + if fac.root === nothing # irreducible quadratic factor + for i = 1:fac.multiplicity + push!(result, (variable(:π’ž, c_idx+=1)*x + variable(:π’ž, c_idx+=1))/(fac.expr^i)) # (Ax + B)/(x-r)^i + end + else + append!(result, variables(:π’ž, (c_idx+1):(c_idx+=fac.multiplicity)) ./ fac.expr.^(1:fac.multiplicity)) # C1/(x-r) + C2/(x-2)^2 ... + end + else + for fac in facs + if fac.root === nothing # irreducible quadratic factor + for i = 1:fac.multiplicity + push!(result, (variable(:π’ž, c_idx+=1)*x + variable(:π’ž, c_idx+=1))/(fac.expr^i)) # (Ax + B)/(x-r)^i + end + continue + end + + # cover up method + other_facs = filter(f -> !isequal(f, fac), facs) + + numerator = rationalize(unwrap(fast_substitute(A / prod((f -> f.expr^f.multiplicity).(other_facs)), Dict(x => fac.root)))) # plug in root to expression without its factor in denominator + push!(result, numerator / fac.expr^fac.multiplicity) + + if fac.multiplicity > 1 + append!(result, variables(:π’ž, (c_idx+1):(c_idx+=fac.multiplicity-1)) ./ fac.expr.^(1:fac.multiplicity-1)) # C1/(x-r) + C2/(x-2)^2 ... + end + end + end + + # no unknowns, so just return + if isequal(get_variables(sum(result)), [x]) + return sum(result ./ leading_coeff) + end + + lhs = numerator(expr) + rhs = expand(sum(simplify.(numerator.(result) .* ((B/leading_coeff) ./ denominator.(result))))) # multiply each numerator by the common denominator/its denominator, and sum to get numerator of whole expression + + solution = solve_interms_ofvar(lhs - rhs, x)[1] # solve for unknowns (C's) by looking at coefficients of the polynomial + + # single unknown + if !(solution isa Dict) + solution = Dict(variable(:π’ž, 1) => solution) + end + + return sum(fast_substitute.(result, Ref(solution)) ./ leading_coeff) # fast_substitute solutions back in and sum +end + +# increasing from 0 to degree n. doesn't skip powers of x like polynomial_coeffs +function coeff_vector(poly, x, n) + coeff_dict = polynomial_coeffs(poly, [x])[1] + vec = [] + for i = 0:n + if x^i in keys(coeff_dict) + push!(vec, coeff_dict[x^i]) + else + push!(vec, 0) + end + end + + return vec +end + +# increasing from 0 to degree of poly +function coeff_vector(poly, x) + return coeff_vector(poly, x, degree(poly)) +end + +function count_multiplicities(facs) + counts = Dict() + for fac in facs + if haskey(counts, fac) + counts[fac] += 1 + else + counts[fac] = 1 + end + end + + return counts +end + +# for partial fractions, into linear and irreducible quadratic factors +function factorize(expr, x) + roots = symbolic_solve(expr, x, dropmultiplicity=false) + + counts = count_multiplicities(roots) + facs = Set() + + for root in keys(counts) + if !isequal(abs(imag(root)), 0) + fac_expr = expand((x - root)*(x - conj(root))) + if !isequal(imag(fac_expr), 0) + @warn "Encountered issue with complex irrational roots. Returning nothing." + return nothing + end + push!(facs, Factor(real(fac_expr), counts[root], x)) + continue + end + + push!(facs, Factor(x - root, root, counts[root], x)) + end + + return facs +end \ No newline at end of file diff --git a/test/diffeqs.jl b/test/diffeqs.jl new file mode 100644 index 000000000..fd8a7173f --- /dev/null +++ b/test/diffeqs.jl @@ -0,0 +1,93 @@ +using Symbolics +using Symbolics: solve_linear_system, LinearODE, has_const_coeffs, to_homogeneous, symbolic_solve_ode, find_particular_solution, IVP, solve_IVP +using Groebner, Nemo +using Test + +@variables x, y, t +Dt = Symbolics.Differential(t) + +# Systems +# ideally, `isapprox` would all be `isequal`, but there seem to be some floating point inaccuracies +@test isapprox(solve_linear_system([1 0; 0 -1], [1, -1], t), [exp(t), -exp(-t)]) +@test isapprox(solve_linear_system([-3 4; -2 3], [7, 2], t), [10exp(-t) - 3exp(t), 5exp(-t) - 3exp(t)]) +@test isapprox(solve_linear_system([4 -3; 8 -6], [7, 2], t), [18 - 11exp(-2t), 24 - 22exp(-2t)]) + +@test_broken isapprox(solve_linear_system([-1 -2; 2 -1], [1, -1], t), [exp(-t)*(cos(2t) + sin(2t)), exp(-t)*(sin(2t) - cos(2t))]) # can't handle complex eigenvalues (though it should be able to) + +@test isapprox(solve_linear_system([1 -1 0; 1 2 1; -2 1 -1], [7, 2, 3], t), (5//3)*exp(-t)*[-1, -2, 7] - 14exp(t)*[-1, 0, 1] + (16//3)*exp(2t)*[-1, 1, 1]) + +@test isequal(solve_linear_system([1 0; 0 -1], [1, -1], t), [exp(t), -exp(-t)]) +@test isequal(solve_linear_system([-3 4; -2 3], [7, 2], t), [10exp(-t) - 3exp(t), 5exp(-t) - 3exp(t)]) +@test isapprox(solve_linear_system([4 -3; 8 -6], [7, 2], t), [18 - 11exp(-2t), 24 - 22exp(-2t)]) + +@test isequal(solve_linear_system([1 -1 0; 1 2 1; -2 1 -1], [7, 2, 3], t), (5//3)*exp(-t)*[-1, -2, 7] - 14exp(t)*[-1, 0, 1] + (16//3)*exp(2t)*[-1, 1, 1]) + +@test_throws ArgumentError solve_linear_system([1 2; 3 4], [1, 2, 3], t) # mismatch between A and x0 +@test_throws ArgumentError solve_linear_system([1 2 3; 4 5 6], [1, 2], t) # A isn't square +@test_throws ArgumentError Symbolics.solve_uncoupled_system([1 2; 3 4], [1, 2], t) # A isn't diagonal + +# LinearODEs +@test Symbolics.is_homogeneous(LinearODE(x, t, [1, 1], 0)) +@test !Symbolics.is_homogeneous(LinearODE(x, t, [t, 1], t^2)) + +@test has_const_coeffs(LinearODE(x, t, [1, 1], 0)) +@test !has_const_coeffs(LinearODE(x, t, [t^2, 1], 0)) + +@test Symbolics.is_homogeneous(to_homogeneous(LinearODE(x, t, [t, 1], t^2))) +@test !Symbolics.is_linear_ode(((Dt^2)(x))^2 ~ x^3, x, t) + +C = Symbolics.variables(:C, 1:5) + +## constant coefficients, nth-order +@test isequal(symbolic_solve_ode(LinearODE(x, t, [-1], 0)), C[1]*exp(t)) +@test isequal(symbolic_solve_ode(LinearODE(x, t, [-4, 3], 0)), C[1]*exp(-4t) + C[2]*exp(t)) + +## first order (solving via integrating factor can be found in test/sympy.jl) +@test isequal(symbolic_solve_ode(LinearODE(x, t, [1], 2sin(t))), C[1]*exp(-t) + sin(t) - cos(t)) + +## repeated characteristic roots +@test isequal(symbolic_solve_ode((Dt^2)(x) + 2(Dt^1)(x) + x ~ 0, x, t), C[1]*exp(-t) + C[2]*t*exp(-t)) +@test isequal(symbolic_solve_ode(LinearODE(x, t, [0, 0, 0, 4, -4], 0)), C[1] + C[2]*t + C[3]*t^2 + C[4]*exp(2t) + C[5]*t*exp(2t)) +@test isequal(symbolic_solve_ode(LinearODE(x, t, [8, 12, 6], 0)), C[1]*exp(-2t) + C[2]*t*exp(-2t) + C[3]*t^2*exp(-2t)) + +## complex characteristic roots +@test isequal(symbolic_solve_ode(LinearODE(x, t, [1, 0], 0)), C[1]*cos(t) + C[2]*sin(t)) +@test isequal(symbolic_solve_ode(LinearODE(x, t, [0, 1, 0], 0)), C[1] + C[2]*cos(t) + C[3]*sin(t)) + +## resonant response formula +@test isequal(symbolic_solve_ode(LinearODE(x, t, [9, -6], 4exp(3t))), C[1]*exp(3t) + C[2]*t*exp(3t) + 2(t^2)*exp(3t)) +### trig functions +@test isequal(symbolic_solve_ode(LinearODE(x, t, [6, 5], 2exp(-t)*cos(t))), C[1]*exp(-2t) + C[2]*exp(-3t) + (1//5)*exp(-t)*cos(t)+(3//5)*exp(-t)*sin(t)) + +## undetermined coefficients +@test isequal(symbolic_solve_ode(LinearODE(x, t, [-3, 2], 2t - 5)), C[1]exp(t) + C[2]exp(-3t) - (2//3)t + 11//9) +@test isequal(find_particular_solution(LinearODE(x, t, [1, 0], t^2)), t^2 - 2) + +# Parsing +@test isequal(LinearODE(x, t, [1], 0), LinearODE(Dt(x) + x ~ 0, x, t)) +@test isequal(LinearODE(x, t, [sin(t), 0, 3t^2], exp(2t) + 2cos(t)), LinearODE(6t^2*(Dt^2)(x) + 2sin(t)*x - 2exp(2t) + 2(Dt^3)(x) ~ 4cos(t), x, t)) + +# IVP +@test isequal(solve_IVP(IVP(LinearODE(x, t, [-3, 2], 0), [1, -1])), (1//2)exp(-3t) + (1//2)exp(t)) +@test isequal(solve_IVP(IVP(LinearODE(x, t, [9, -6], 4exp(3t)), [5, 6])), 5exp(3t) - 9t*exp(3t) + 2(t^2)*exp(3t)) + +# Other methods + +## Clairaut's equation +@test isequal(symbolic_solve_ode(x ~ Dt(x)*t - ((Dt(x))^3), x, t), C[1]*t - C[1]^3) +@test isequal(symbolic_solve_ode(Dt(x)*t + (Dt(x))^2 - sin(Dt(x)) + 2 ~ x, x, t), C[1]*t + C[1]^2 - sin(C[1]) + 2) +@test isnothing(symbolic_solve_ode(Dt(x) + (Dt(x))^2 ~ x, x, t)) +@test isnothing(symbolic_solve_ode(Dt(x)*t + 2t*(Dt(x))^2 ~ x, x, t)) +@test isnothing(symbolic_solve_ode(Dt(x) + x*(Dt(x))^2 ~ x, x, t)) + +## Bernoulli equations (integrating factor solve in test/sympy.jl) +@test isequal(symbolic_solve_ode(Dt(x) - 5x ~ exp(-2t)*x^(-2), x, t), (C[1]exp(15t) - (3//17)exp(-2t))^(1//3)) +@test isnothing(symbolic_solve_ode(sqrt((Dt^4)(x)) ~ log(x)^t, x, t)) + +# Helper function tests +ys = Symbolics.variables(:y, 1:2) +@test isequal(Symbolics.reduce_order((Dt^2)(x) + 3Dt(x) + 2x ~ 0, x, t, ys), [ys[2], -2ys[1] - 3ys[2]]) +@test isequal(Symbolics.unreduce_order([ys[1], ys[2]], x, t, ys), [x, Dt(x)]) + +@test Symbolics.is_solution(C[1]*exp(3t) + C[2]*t*exp(3t) + 2(t^2)*exp(3t), LinearODE(x, t, [9, -6], 4exp(3t))) +@test Symbolics.is_solution(C[1]*exp(-t) + C[2]*t*exp(-t), (Dt^2)(x) + 2(Dt^1)(x) + x ~ 0, x, t) \ No newline at end of file diff --git a/test/laplace.jl b/test/laplace.jl new file mode 100644 index 000000000..798e45a07 --- /dev/null +++ b/test/laplace.jl @@ -0,0 +1,37 @@ +using Test +using Symbolics +using Symbolics: laplace, inverse_laplace +import Nemo, Groebner + +@variables t, s +@syms f(t)::Real F(s)::Real + +# https://sites.math.washington.edu/~aloveles/Math307Fall2019/m307LaplacePractice.pdf +@test isequal(laplace(exp(4t) + 5, f, t, F, s), 1/(s-4) + 5/s) +@test isequal(laplace(cos(2t) + 7sin(2t), f, t, F, s), s/(s^2 + 4) + 14/(s^2 + 4)) +@test isequal(laplace(exp(-2t)*cos(3t) + 5exp(-2t)*sin(3t), f, t, F, s), (s+2)/((s+2)^2 + 9) + 15/((s+2)^2 + 9)) +@test isequal(laplace(10 + 5t + t^2 - 4t^3, f, t, F, s), expand(10/s + 5/s^2 + 2/s^3 - 24/s^4)) +@test isequal(laplace(exp(3t)*(t^2 + 4t + 2), f, t, F, s), 2/(s-3)^3 + 4/(s-3)^2 + 2/(s-3)) +@test isequal(laplace(6exp(5t)*cos(2t) - exp(7t), f, t, F, s), 6(s-5)/((s-5)^2 + 4) + expand(-1/(s-7))) + +# https://www.math.lsu.edu/~adkins/m2065/2065s08review2a.pdf +@test isequal(inverse_laplace(7/(s+3)^3, F, s, f, t), (7//2)t^2 * exp(-3t)) +@test isequal(inverse_laplace((s-9)/(s^2 + 9), F, s, f, t), cos(3t) - 3sin(3t)) +# partial fraction decomposition +@test isequal(inverse_laplace((s+2)/(s^2 - 3s - 4), F, s, f, t), (6//5)*exp(4t) - (1//5)*exp(-t)) +@test isequal(inverse_laplace(1/(s^2 - 10s + 9), F, s, f, t), (1//8)*exp(9t) - (1//8)*exp(t)) + +Dt = Differential(t) +@test isequal(laplace_solve_ode(Dt(f(t)) + 3f(t) ~ t^2*exp(-3t) + t*exp(-2t) + t, f, t, [1]), (1//3)*t^3*exp(-3t) + t*exp(-2t) + (1//3)*t + (19//9)*exp(-3t) - exp(-2t) - 1//9) +@test isequal(laplace_solve_ode((Dt^2)(f(t)) - 3Dt(f(t)) + 2f(t) ~ 4, f, t, [2, 3]), 2 - 3exp(t) + 3exp(2t)) +@test isequal(laplace_solve_ode((Dt^3)(f(t)) - Dt(f(t)) ~ 2, f, t, [4,4,4]), 5exp(t) - exp(-t) - 2t) +@test isequal(laplace_solve_ode((Dt^3)(f(t)) - Dt(f(t)) ~ 6 - 3t^2, f, t, [1, 1, 1]), exp(t) + t^3) +@test isequal(laplace_solve_ode((Dt^2)(f(t)) - f(t) ~ 2sin(t), f, t, [0, 0]), (1//2)exp(t) - (1//2)exp(-t) - sin(t)) +@test isequal(laplace_solve_ode((Dt^2)(f(t)) + 2Dt(f(t)) ~ 5f(t), f, t, [0, 0]), 0) +@test isequal(laplace_solve_ode((Dt^2)(f(t)) + f(t) ~ sin(4t), f, t, [0, 0]), (4//15)sin(t) - (1//15)sin(4t)) +@test isequal(laplace_solve_ode((Dt^2)(f(t)) + Dt(f(t)) ~ 1 + 2t, f, t, [0, 0]), 1 - exp(-t) + t^2 - t) +@test isequal(laplace_solve_ode((Dt^2)(f(t)) + 4Dt(f(t)) + 3f(t) ~ 6, f, t, [0, 0]), exp(-3t) - 3exp(-t) + 2) +@test isequal(laplace_solve_ode((Dt^2)(f(t)) - 2Dt(f(t)) ~ 3*(t + exp(2t)), f, t, [0, 0]), (3//8) - (3//4)t - (3//4)t^2 - (3//8)exp(2t) + (3//2)t*exp(2t)) +@test_broken isequal(laplace_solve_ode((Dt^2)(f(t)) - 2Dt(f(t)) ~ 20*exp(-t)*cos(t), f, t, [0, 0]), 3exp(2t) - 5 + 2exp(-t)*cos(t) - 4exp(-t)*sin(t)) # irreducible quadratic in inverse laplace +@test isequal(laplace_solve_ode((Dt^2)(f(t)) + f(t) ~ 2 + 2cos(t), f, t, [0, 0]), 2 - 2cos(t) + t*sin(t)) +@test isequal(laplace_solve_ode((Dt^2)(f(t)) - Dt(f(t)) ~ 30cos(3t), f, t, [0, 0]), 3exp(t) - 3cos(3t) - sin(3t)) \ No newline at end of file diff --git a/test/partialfractions.jl b/test/partialfractions.jl new file mode 100644 index 000000000..d6c75e814 --- /dev/null +++ b/test/partialfractions.jl @@ -0,0 +1,26 @@ +using Test +using Symbolics +import Nemo, Groebner +import Symbolics: partial_frac_decomposition + +@variables x + +# https://en.neurochispas.com/algebra/4-types-of-partial-fractions-decomposition-with-examples/ +@test isequal(partial_frac_decomposition((3x-1) / (x^2 + x - 6), x), 2/(x+3) + 1/(x-2)) +@test isequal(partial_frac_decomposition((9x^2 + 34x + 14) / ((x+2)*(x^2 - x - 12)), x), expand(3/(x+2) + 7/(x-4) - 1/(x+3))) + +# https://tutorial.math.lamar.edu/Problems/Alg/PartialFractions.aspx +@test isequal(partial_frac_decomposition((17x-53)/(x^2 - 2x - 15), x), expand(4/(x-5) + 13/(x+3))) +@test isequal(partial_frac_decomposition((34-12x)/(3x^2 - 10x - 8), x), (-3)/(2//3 + x) + -1/(-4 + x)) +@test isequal(partial_frac_decomposition((125 + 4x - 9x^2)/((x-1)*(x+3)*(x+4)), x), expand(6/(x-1) - 8/(x+3) - 7/(x+4))) +@test isequal(partial_frac_decomposition((10x+35)/((x+4)^2), x), 10/(x+4) + -5/(x+4)^2) +@test isequal(partial_frac_decomposition((6x+5)/((2x-1)^2), x), (3//2)/(x-1//2) + 2/(x-1//2)^2) +@test isequal(partial_frac_decomposition((7x^2-17x+38)/((x+6)*(x-1)^2), x), 8/(x+6) + -1/(x-1) + 4/(x-1)^2) +@test isequal(partial_frac_decomposition((4x^2 - 22x + 7)/((2x+3)*(x-2)^2), x), 2/(x+3//2) + -3/(x-2)^2) +@test_broken isequal(partial_frac_decomposition((3x^2 + 7x + 28)/(x*(x^2 + x + 7)), x), expand(4/x + (3-x)/(x^2+x+7))) # irrational roots +@test isequal(partial_frac_decomposition((4x^3 + 16x + 7)/(x^2 + 4)^2, x), 4x/(x^2+4) + 7/(x^2+4)^2) + +# check valid expressions +@test partial_frac_decomposition(sin(x), x) === nothing +@test partial_frac_decomposition(x^2/(x-1), x) === nothing +@test partial_frac_decomposition(1/(x^2 + 2), x) === nothing # irrational roots, should eventually be fixed \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index cc8d213f4..9e7271716 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -70,6 +70,9 @@ if GROUP == "All" || GROUP == "Core" @safetestset "Function inverses test" begin include("inverse.jl") end @safetestset "Taylor Series Test" begin include("taylor.jl") end @safetestset "Discontinuity registration test" begin include("discontinuities.jl") end + @safetestset "ODE solver test" begin include("diffeqs.jl") end + @safetestset "Laplace transform test" begin include("laplace.jl") end + @safetestset "Partial Fraction Decomposition Test" begin include("partialfractions.jl") end end end diff --git a/test/sympy.jl b/test/sympy.jl index 486872f52..3617b8520 100644 --- a/test/sympy.jl +++ b/test/sympy.jl @@ -80,6 +80,21 @@ expected_sol = C1 * exp(2 * x) canonical_sol_ode = Symbolics.substitute(sol_ode, Dict(const_sym => C1)) @test isequal(canonical_sol_ode, expected_sol) +## Native ODE solver, but using sympy_integrate +@variables x, t +Dt = Symbolics.Differential(t) +C = Symbolics.variables(:C, 1:5) + +@test isequal(symbolic_solve_ode(LinearODE(x, t, [5/t], 7t)), Symbolics.sympy_simplify(C[1]*t^(-5) + t^2)) +@test isequal(symbolic_solve_ode(LinearODE(x, t, [cos(t)], cos(t))), 1 + C[1]*exp(-sin(t))) +@test isequal(symbolic_solve_ode(LinearODE(x, t, [-(1+t)], 1+t)), Symbolics.expand(Symbolics.sympy_simplify(C[1]*exp((1//2)t^2 + t) - 1))) +# SymPy is being weird and not simplifying correctly (and some symbols are wrong, like pi and erf being syms), but these otherwise work +@test_broken isequal(symbolic_solve_ode(LinearODE(x, t, [-2t], 1)), Symbolics.sympy_simplify(exp(t^2)*sqrt(Symbolics.variable(:pi))*erf(t)/2 + C[1]*exp(t^2))) + +## Bernoulli equations +@test isequal(symbolic_solve_ode(Dt(x) + (4//t)*x ~ t^3 * x^2, x, t), 1/(C[1]t^4 - t^4 * log(t))) +@test isnothing(symbolic_solve_ode(Dt(x) + sqrt(t)*x ~ sin(4t)*x^3, x, t)) + # Test issue #1605: Power with symbolic exponent @variables x y expr_power = x^y