From f3b9477ffc22c1f72abfe31cf0e7daf6fef8c792 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Wed, 25 Jun 2025 22:16:54 -0400 Subject: [PATCH 01/19] Created basic solvers for 1st order separable eqs and linear systems --- src/diffeqs.jl | 145 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 145 insertions(+) create mode 100644 src/diffeqs.jl diff --git a/src/diffeqs.jl b/src/diffeqs.jl new file mode 100644 index 000000000..e9b962300 --- /dev/null +++ b/src/diffeqs.jl @@ -0,0 +1,145 @@ +using Groebner, Nemo + +""" +Solve first order separable ODE + +(mostly me getting used to Symbolics, not super useful in practice) + +For example, dx/dt + p(t)x ~ 0 +""" +function firstorder_separable_ode_solve(ex, x, t) + x, t = Symbolics.value(x), Symbolics.value(t) + p = Symbolics.coeff(ex.lhs, x) # function of t + P = sympy_integrate(p, t) + @variables C + return simplify(C * exp(-P)) +end + +""" +Returns evolution matrix e^(tD) +""" +evo_mat(D::Matrix{<:Number}, t::Num) = diagm(exp.(t .* diag(D))) + +""" +Solve linear continuous dynamical system of differential equations of the form Ax = x' with initial condition x0 + +# Arguments +- `A`: matrix of coefficients +- `x0`: intial conditions vector +- `t`: independent variable + +# Returns +- vector of symbolic solutions +""" +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) * 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 = Symbolics.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{Num} = 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 = 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{Num} = 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 + +# tests +@variables x, t + +Dt = Differential(t) +ex = Dt(x) + 2 * t * x ~ 0 +println(firstorder_separable_ode_solve(ex, x, t)) +println() + +A = [1 0; 0 -1] +x0 = [1, -1] +println(solve_uncoupled_system(A, x0, t)) +println() + +# commented out below because currently can't handle complex eigenvalues +# A = [-1 -2; 2 -1] +# x0 = [1, -1] +# println(solve_linear_system(A, x0, t)) + +# println() +# println(symbolic_eigen([-3 4; -2 3])) # should be [1, -1] and [1 2; 1 1] (or equivalent) +# println() +# println(symbolic_eigen([4 -3; 8 -6])) # should be [-2, 0] and [1 2; 3 4] (or equivalent) +# println() +# println(symbolic_eigen([1 -1 0; 1 2 1; -2 1 -1])) # should be [-1, 1, 2] and [-1 -1 -1; -2 0 1; 7 1 1] (or equivalent) +# println() + +println(solve_linear_system([-3 4; -2 3], [7, 2], t)) +println() +println(solve_linear_system([4 -3; 8 -6], [7, 2], t)) +println() +# println(inv(diagonalize([1 -1 0; 1 2 1; -2 1 -1])[1])) --- shows that inv isn't maintaining symbolics +println(solve_linear_system([1 -1 0; 1 2 1; -2 1 -1], [7, 2, 3], t)) +println() From e55dc3481d8f103fb9c8c2610ed0a725c24a215a Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Thu, 26 Jun 2025 03:26:52 -0400 Subject: [PATCH 02/19] added rationalization to inverted matrix in solve_linear_system --- src/diffeqs.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/diffeqs.jl b/src/diffeqs.jl index e9b962300..c0da19a05 100644 --- a/src/diffeqs.jl +++ b/src/diffeqs.jl @@ -49,7 +49,7 @@ function solve_linear_system(A::Matrix{<:Number}, x0::Vector{<:Number}, t::Num) S, D = diagonalize(A) - return simplify.(S * evo_mat(D, t) * S^-1 * x0) + return simplify.(S * evo_mat(D, t) * rationalize.(S^-1) * x0) end """ @@ -91,7 +91,7 @@ function symbolic_eigen(A::Matrix{<:Number}) values = symbolic_solve(p, λ) # solve polynomial # then, find eigenvectors - S::Matrix{Num} = Matrix(I, size(A, 1), 0) # matrix storing vertical 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 @@ -103,7 +103,7 @@ function symbolic_eigen(A::Matrix{<:Number}) if sol[1] isa Dict sol = [sol[1][var] for var in v[2:end]] end - vec::Vector{Num} = prepend!(sol, [1]) # add back the 1 (representing v_1) from substitution + vec::Vector{Number} = prepend!(sol, [1]) # add back the 1 (representing v_1) from substitution S = [S vec] # add vec to matrix end From 129803b9ac4d498b03e525b7106798bd3885db6d Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Fri, 27 Jun 2025 10:15:36 -0400 Subject: [PATCH 03/19] integrated into Symbolics, restructured files --- src/Symbolics.jl | 5 +++ src/diffeqs/diffeqs.jl | 31 +++++++++++++++ src/{diffeqs.jl => diffeqs/systems.jl} | 55 ++------------------------ 3 files changed, 39 insertions(+), 52 deletions(-) create mode 100644 src/diffeqs/diffeqs.jl rename src/{diffeqs.jl => diffeqs/systems.jl} (64%) diff --git a/src/Symbolics.jl b/src/Symbolics.jl index f192eb77c..9012b7316 100644 --- a/src/Symbolics.jl +++ b/src/Symbolics.jl @@ -218,6 +218,11 @@ include("solver/main.jl") include("solver/special_cases.jl") export symbolic_solve +# Diff Eq Solver +include("diffeqs/diffeqs.jl") +include("diffeqs/systems.jl") +export firstorder_separable_ode_solve, solve_linear_system + # Sympy Functions """ diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl new file mode 100644 index 000000000..70b6efbff --- /dev/null +++ b/src/diffeqs/diffeqs.jl @@ -0,0 +1,31 @@ +using Symbolics +import Symbolics: value, coeff, sympy_integrate + +struct LinearODE + # dⁿx/dtⁿ + pₙ(t)(dⁿ⁻¹x/dtⁿ⁻¹) + ... + p₂(t)(dx/dt) + p₁(t)x = q(t) + + x::SymbolicUtils.Symbolic # dependent variable + t::SymbolicUtils.Symbolic # independent variable + p::AbstractArray # coefficient functions of t ordered in increasing order (p₁, p₂, ...) + q # right hand side function of t, without any x + + LinearODE(x::Num, t::Num, p, q) = new(value(x), value(t), p, q) +end + + +is_homogeneous(eq::LinearODE) = isempty(Symbolics.get_variables(eq.q)) + +""" +Solve first order separable ODE + +(mostly me getting used to Symbolics, not super useful in practice) + +For example, dx/dt + p(t)x ~ 0 +""" +function firstorder_separable_ode_solve(ex, x, t) + x, t = Symbolics.value(x), Symbolics.value(t) + p = Symbolics.coeff(ex.lhs, x) # function of t + P = Symbolics.sympy_integrate(p, t) + @variables C + return simplify(C * exp(-P)) +end \ No newline at end of file diff --git a/src/diffeqs.jl b/src/diffeqs/systems.jl similarity index 64% rename from src/diffeqs.jl rename to src/diffeqs/systems.jl index c0da19a05..521cfdf06 100644 --- a/src/diffeqs.jl +++ b/src/diffeqs/systems.jl @@ -1,19 +1,4 @@ -using Groebner, Nemo - -""" -Solve first order separable ODE - -(mostly me getting used to Symbolics, not super useful in practice) - -For example, dx/dt + p(t)x ~ 0 -""" -function firstorder_separable_ode_solve(ex, x, t) - x, t = Symbolics.value(x), Symbolics.value(t) - p = Symbolics.coeff(ex.lhs, x) # function of t - P = sympy_integrate(p, t) - @variables C - return simplify(C * exp(-P)) -end +# import Symbolics: symbolic_solve """ Returns evolution matrix e^(tD) @@ -84,7 +69,7 @@ Replacement for `LinearAlgebra.eigen` function that uses symbolic functions to a """ function symbolic_eigen(A::Matrix{<:Number}) @variables λ # eigenvalue - v = Symbolics.variables(:v, 1:size(A, 1)) # vector of subscripted variables to represent eigenvector + 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 @@ -108,38 +93,4 @@ function symbolic_eigen(A::Matrix{<:Number}) end return Eigen(values, S) -end - -# tests -@variables x, t - -Dt = Differential(t) -ex = Dt(x) + 2 * t * x ~ 0 -println(firstorder_separable_ode_solve(ex, x, t)) -println() - -A = [1 0; 0 -1] -x0 = [1, -1] -println(solve_uncoupled_system(A, x0, t)) -println() - -# commented out below because currently can't handle complex eigenvalues -# A = [-1 -2; 2 -1] -# x0 = [1, -1] -# println(solve_linear_system(A, x0, t)) - -# println() -# println(symbolic_eigen([-3 4; -2 3])) # should be [1, -1] and [1 2; 1 1] (or equivalent) -# println() -# println(symbolic_eigen([4 -3; 8 -6])) # should be [-2, 0] and [1 2; 3 4] (or equivalent) -# println() -# println(symbolic_eigen([1 -1 0; 1 2 1; -2 1 -1])) # should be [-1, 1, 2] and [-1 -1 -1; -2 0 1; 7 1 1] (or equivalent) -# println() - -println(solve_linear_system([-3 4; -2 3], [7, 2], t)) -println() -println(solve_linear_system([4 -3; 8 -6], [7, 2], t)) -println() -# println(inv(diagonalize([1 -1 0; 1 2 1; -2 1 -1])[1])) --- shows that inv isn't maintaining symbolics -println(solve_linear_system([1 -1 0; 1 2 1; -2 1 -1], [7, 2, 3], t)) -println() +end \ No newline at end of file From a5b7c51f9ca2d6b0f354ecced891147b27683a5f Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Fri, 27 Jun 2025 10:20:06 -0400 Subject: [PATCH 04/19] added unit tests --- test/diffeqs.jl | 26 ++++++++++++++++++++++++++ test/runtests.jl | 1 + 2 files changed, 27 insertions(+) create mode 100644 test/diffeqs.jl diff --git a/test/diffeqs.jl b/test/diffeqs.jl new file mode 100644 index 000000000..619f883b6 --- /dev/null +++ b/test/diffeqs.jl @@ -0,0 +1,26 @@ +using Symbolics +using Symbolics: firstorder_separable_ode_solve, solve_linear_system, LinearODE, is_homogeneous +import Groebner, Nemo, SymPy +using Test + +@variables x, y, t, C +Dt = Differential(t) +@test_broken isequal(firstorder_separable_ode_solve(Dt(x) - x ~ 0, x, t), C*exp(t)) +@test isequal(firstorder_separable_ode_solve(Dt(x) + 2*t*x ~ 0, x, t), C*exp(-(t^2))) +@test isequal(firstorder_separable_ode_solve(Dt(x) + (t^2-3)*x ~ 0, x, t), C*exp(3t - (1//3)*t^3)) + +# Systems +@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))]) +@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_broken isequal(solve_linear_system([4 -3; 8 -6], [7, 2], t), [18 - 11exp(-2t), 24 - 22exp(-2t)]) +@test_broken isequal(solve_linear_system([-1 -2; 2 -1], [1, -1], t), [exp(-t)*(cos(2t) + sin(2t)), exp(-t)*(sin(2t) - cos(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]) + +# LinearODEs +@test is_homogeneous(LinearODE(x, t, [1, 1], 0)) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 2df5bd8ff..909739d8c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -63,6 +63,7 @@ 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 end end From d6daa7b6e5ab9b965adec20631cfad90e0434052 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Fri, 27 Jun 2025 11:59:43 -0400 Subject: [PATCH 05/19] added display functionality for LinearODE --- src/diffeqs/diffeqs.jl | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 70b6efbff..732b70a32 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -8,10 +8,16 @@ struct LinearODE t::SymbolicUtils.Symbolic # independent variable p::AbstractArray # coefficient functions of t ordered in increasing order (p₁, p₂, ...) q # right hand side function of t, without any x + Dt::Differential + order::Int LinearODE(x::Num, t::Num, p, q) = new(value(x), value(t), p, q) end +get_expression(eq::LinearODE) = (eq.Dt^eq.order)(eq.x) + sum([(eq.p[n])*(eq.Dt^n)(eq.x) for n = 1:length(eq.p)]) ~ eq.q + +Base.print(io::IO, eq::LinearODE) = print(io, "(D$(eq.t)^$(eq.order))$(eq.x) + " * join(["($(eq.p[length(eq.p)-n]))(D$(eq.t)^$(length(eq.p)-n))$(eq.x)" for n = 0:(eq.order-2)], " + ") * " ~ $(eq.q)") +Base.show(io::IO, eq::LinearODE) = print(io, eq) is_homogeneous(eq::LinearODE) = isempty(Symbolics.get_variables(eq.q)) From 871fff4ad9d2010e80ce40174af8fd345405c7fd Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Fri, 27 Jun 2025 12:00:27 -0400 Subject: [PATCH 06/19] implemented solver for homogeneous solutions (of degree 5 or less with no repeated roots) --- src/diffeqs/diffeqs.jl | 25 ++++++++++++++++++++++++- test/diffeqs.jl | 17 +++++++++++++++-- 2 files changed, 39 insertions(+), 3 deletions(-) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 732b70a32..1139f9d00 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -11,7 +11,8 @@ struct LinearODE Dt::Differential order::Int - LinearODE(x::Num, t::Num, p, q) = new(value(x), value(t), p, q) + LinearODE(x::Num, t::Num, p, q) = new(value(x), value(t), p, q, Differential(t), length(p)+1) + LinearODE(x, t, p, q) = new(x, t, p, q, Differential(t), length(p)+1) end get_expression(eq::LinearODE) = (eq.Dt^eq.order)(eq.x) + sum([(eq.p[n])*(eq.Dt^n)(eq.x) for n = 1:length(eq.p)]) ~ eq.q @@ -20,6 +21,28 @@ Base.print(io::IO, eq::LinearODE) = print(io, "(D$(eq.t)^$(eq.order))$(eq.x) + " Base.show(io::IO, eq::LinearODE) = print(io, eq) is_homogeneous(eq::LinearODE) = isempty(Symbolics.get_variables(eq.q)) +has_const_coeffs(eq::LinearODE) = all(isempty.(Symbolics.get_variables.(eq.p))) + +to_homogeneous(eq::LinearODE) = LinearODE(eq.x, eq.t, eq.p, 0) + +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 + +function homogeneous_solve(eq::LinearODE) + @variables r + p = characteristic_polynomial(eq, r) + roots = symbolic_solve(p, r, dropmultiplicity=false) + @variables C[1:degree(p, r)] + return sum(Symbolics.scalarize(C) .* exp.(roots*eq.t)) +end """ Solve first order separable ODE diff --git a/test/diffeqs.jl b/test/diffeqs.jl index 619f883b6..c280ddbcb 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -1,5 +1,5 @@ using Symbolics -using Symbolics: firstorder_separable_ode_solve, solve_linear_system, LinearODE, is_homogeneous +using Symbolics: firstorder_separable_ode_solve, solve_linear_system, LinearODE, is_homogeneous, has_const_coeffs, to_homogeneous, homogeneous_solve import Groebner, Nemo, SymPy using Test @@ -13,14 +13,27 @@ Dt = Differential(t) @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))]) + @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_broken isequal(solve_linear_system([4 -3; 8 -6], [7, 2], t), [18 - 11exp(-2t), 24 - 22exp(-2t)]) + @test_broken isequal(solve_linear_system([-1 -2; 2 -1], [1, -1], t), [exp(-t)*(cos(2t) + sin(2t)), exp(-t)*(sin(2t) - cos(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]) # LinearODEs -@test is_homogeneous(LinearODE(x, t, [1, 1], 0)) \ No newline at end of file +@test is_homogeneous(LinearODE(x, t, [1, 1], 0)) +@test !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 is_homogeneous(to_homogeneous(LinearODE(x, t, [t, 1], t^2))) +@variables C[1:3] +@test isequal(homogeneous_solve(LinearODE(x, t, [-1], 0)), C[1]*exp(t)) +@test isequal(homogeneous_solve(LinearODE(x, t, [-4, 3], 0)), C[1]*exp(-4t) + C[2]*exp(t)) \ No newline at end of file From b7a6c2e32a51abaf22f464300c25ff5ae3d50a14 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Fri, 27 Jun 2025 16:48:15 -0400 Subject: [PATCH 07/19] added solving via integrating factor (1st order ODEs) --- src/Symbolics.jl | 2 +- src/diffeqs/diffeqs.jl | 57 +++++++++++++++++++++++++++++------------- test/diffeqs.jl | 25 +++++++++++------- 3 files changed, 56 insertions(+), 28 deletions(-) diff --git a/src/Symbolics.jl b/src/Symbolics.jl index 9012b7316..d06ca5759 100644 --- a/src/Symbolics.jl +++ b/src/Symbolics.jl @@ -221,7 +221,7 @@ export symbolic_solve # Diff Eq Solver include("diffeqs/diffeqs.jl") include("diffeqs/systems.jl") -export firstorder_separable_ode_solve, solve_linear_system +export symbolic_solve_ode, solve_linear_system # Sympy Functions diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 1139f9d00..a14156fe3 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -10,14 +10,15 @@ struct LinearODE q # right hand side function of t, without any x Dt::Differential order::Int + C::Vector{Num} # constants - LinearODE(x::Num, t::Num, p, q) = new(value(x), value(t), p, q, Differential(t), length(p)+1) - LinearODE(x, t, p, q) = new(x, t, p, q, Differential(t), length(p)+1) + LinearODE(x::Num, t::Num, p, q) = new(value(x), value(t), p, q, Differential(t), length(p), variables(:C, 1:length(p))) + LinearODE(x, t, p, q) = new(x, t, p, q, Differential(t), length(p), variables(:C, 1:length(p))) end -get_expression(eq::LinearODE) = (eq.Dt^eq.order)(eq.x) + sum([(eq.p[n])*(eq.Dt^n)(eq.x) for n = 1:length(eq.p)]) ~ eq.q +get_expression(eq::LinearODE) = (eq.Dt^eq.order)(eq.x) + sum([(eq.p[n])*(eq.Dt^(n-1))(eq.x) for n = 1:length(eq.p)]) ~ eq.q -Base.print(io::IO, eq::LinearODE) = print(io, "(D$(eq.t)^$(eq.order))$(eq.x) + " * join(["($(eq.p[length(eq.p)-n]))(D$(eq.t)^$(length(eq.p)-n))$(eq.x)" for n = 0:(eq.order-2)], " + ") * " ~ $(eq.q)") +Base.print(io::IO, eq::LinearODE) = print(io, "(D$(eq.t)^$(eq.order))$(eq.x) + " * join(["($(eq.p[length(eq.p)-n]))(D$(eq.t)^$(length(eq.p)-n-1))$(eq.x)" for n = 0:(eq.order-1)], " + ") * " ~ $(eq.q)") Base.show(io::IO, eq::LinearODE) = print(io, eq) is_homogeneous(eq::LinearODE) = isempty(Symbolics.get_variables(eq.q)) @@ -36,25 +37,45 @@ function characteristic_polynomial(eq::LinearODE, r) return poly end -function homogeneous_solve(eq::LinearODE) +""" +Symbolically solve a linear ODE + +Cases handled: +- ☑ first order +- ☑ homogeneous with constant coefficients +- ▢ nonhomogeneous with constant coefficients +- ▢ particular solutions (variation of parameters? undetermined coefficients?) +- ▢ [Differential transform method](https://www.researchgate.net/publication/267767445_A_New_Algorithm_for_Solving_Linear_Ordinary_Differential_Equations) +""" +function symbolic_solve_ode(eq::LinearODE) + if eq.order == 1 + return integrating_factor_solve(eq) + end + + if is_homogeneous(eq) + if has_const_coeffs(eq) + return const_coeff_solve + end + end +end + +function const_coeff_solve(eq::LinearODE) @variables r p = characteristic_polynomial(eq, r) roots = symbolic_solve(p, r, dropmultiplicity=false) - @variables C[1:degree(p, r)] - return sum(Symbolics.scalarize(C) .* exp.(roots*eq.t)) + return sum(eq.C .* exp.(roots*eq.t)) end """ -Solve first order separable ODE - -(mostly me getting used to Symbolics, not super useful in practice) - -For example, dx/dt + p(t)x ~ 0 +Solve almost any first order ODE using an integrating factor """ -function firstorder_separable_ode_solve(ex, x, t) - x, t = Symbolics.value(x), Symbolics.value(t) - p = Symbolics.coeff(ex.lhs, x) # function of t - P = Symbolics.sympy_integrate(p, t) - @variables C - return simplify(C * exp(-P)) +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 + return Symbolics.sympy_simplify((1/v) * (sympy_integrate(eq.q*v, eq.t) + eq.C[1])) end \ No newline at end of file diff --git a/test/diffeqs.jl b/test/diffeqs.jl index c280ddbcb..ff066f36f 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -1,13 +1,9 @@ using Symbolics -using Symbolics: firstorder_separable_ode_solve, solve_linear_system, LinearODE, is_homogeneous, has_const_coeffs, to_homogeneous, homogeneous_solve +using Symbolics: solve_linear_system, LinearODE, is_homogeneous, has_const_coeffs, to_homogeneous, symbolic_solve_ode import Groebner, Nemo, SymPy using Test -@variables x, y, t, C -Dt = Differential(t) -@test_broken isequal(firstorder_separable_ode_solve(Dt(x) - x ~ 0, x, t), C*exp(t)) -@test isequal(firstorder_separable_ode_solve(Dt(x) + 2*t*x ~ 0, x, t), C*exp(-(t^2))) -@test isequal(firstorder_separable_ode_solve(Dt(x) + (t^2-3)*x ~ 0, x, t), C*exp(3t - (1//3)*t^3)) +@variables x, y, t # Systems @test isapprox(solve_linear_system([1 0; 0 -1], [1, -1], t), [exp(t), -exp(-t)]) @@ -34,6 +30,17 @@ Dt = Differential(t) @test !has_const_coeffs(LinearODE(x, t, [t^2, 1], 0)) @test is_homogeneous(to_homogeneous(LinearODE(x, t, [t, 1], t^2))) -@variables C[1:3] -@test isequal(homogeneous_solve(LinearODE(x, t, [-1], 0)), C[1]*exp(t)) -@test isequal(homogeneous_solve(LinearODE(x, t, [-4, 3], 0)), C[1]*exp(-4t) + C[2]*exp(t)) \ No newline at end of file + +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 +@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.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))*((@syms erf(z))[1])(t)/2 + C[1]*exp(t^2))) +@test_broken isequal(symbolic_solve_ode(LinearODE(x, t, [1], 2sin(t))), C[1]*exp(-t) + sin(t) - cos(t)) \ No newline at end of file From ee8bc8e6764626cb5ee2230b81346766f43f3085 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Fri, 27 Jun 2025 20:31:52 -0400 Subject: [PATCH 08/19] added handling of repeated characteristic roots --- src/diffeqs/diffeqs.jl | 17 ++++++++++++++--- test/diffeqs.jl | 10 ++++++++-- 2 files changed, 22 insertions(+), 5 deletions(-) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index a14156fe3..8863278e5 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -54,7 +54,7 @@ function symbolic_solve_ode(eq::LinearODE) if is_homogeneous(eq) if has_const_coeffs(eq) - return const_coeff_solve + return const_coeff_solve(eq) end end end @@ -63,7 +63,18 @@ function const_coeff_solve(eq::LinearODE) @variables r p = characteristic_polynomial(eq, r) roots = symbolic_solve(p, r, dropmultiplicity=false) - return sum(eq.C .* exp.(roots*eq.t)) + + # Handle repeated roots + solutions = exp.(roots*eq.t) + for i in eachindex(solutions)[1:end-1] + j = i+1 + while j <= length(solutions) && isequal(solutions[i], solutions[j]) + solutions[j] *= eq.t^(j-i) # multiply by t for each repetition + j+=1 + end + end + + return sum(eq.C .* solutions) end """ @@ -77,5 +88,5 @@ function integrating_factor_solve(eq::LinearODE) else v = exp(sympy_integrate(p, eq.t)) end - return Symbolics.sympy_simplify((1/v) * (sympy_integrate(eq.q*v, eq.t) + eq.C[1])) + return Symbolics.sympy_simplify((1/v) * ((isequal(eq.q, 0) ? 0 : sympy_integrate(eq.q*v, eq.t)) + eq.C[1])) end \ No newline at end of file diff --git a/test/diffeqs.jl b/test/diffeqs.jl index ff066f36f..154378df8 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -42,5 +42,11 @@ C = Symbolics.variables(:C, 1:5) @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.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))*((@syms erf(z))[1])(t)/2 + C[1]*exp(t^2))) -@test_broken isequal(symbolic_solve_ode(LinearODE(x, t, [1], 2sin(t))), C[1]*exp(-t) + sin(t) - cos(t)) \ No newline at end of file +@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))) +@test_broken 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(LinearODE(x, t, [1, 2], 0)), 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)) + From 2376fa943964e12017a52bd44f79331d7f399d3e Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Sat, 28 Jun 2025 16:29:45 -0400 Subject: [PATCH 09/19] implemented exponential response formula / resonant response formula --- src/diffeqs/diffeqs.jl | 80 +++++++++++++++++++++++++++++++++++++++++- test/diffeqs.jl | 1 + 2 files changed, 80 insertions(+), 1 deletion(-) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 8863278e5..7123d1e36 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -43,7 +43,8 @@ Symbolically solve a linear ODE Cases handled: - ☑ first order - ☑ homogeneous with constant coefficients -- ▢ nonhomogeneous with constant coefficients +- ◩ nonhomogeneous with constant coefficients + - ☑ ERF + RRF - ▢ particular solutions (variation of parameters? undetermined coefficients?) - ▢ [Differential transform method](https://www.researchgate.net/publication/267767445_A_New_Algorithm_for_Solving_Linear_Ordinary_Differential_Equations) """ @@ -57,6 +58,13 @@ function symbolic_solve_ode(eq::LinearODE) return const_coeff_solve(eq) end end + + if has_const_coeffs(eq) + rrf = resonant_response_formula(eq) + if rrf !== nothing + return const_coeff_solve(to_homogeneous(eq)) + rrf + end + end end function const_coeff_solve(eq::LinearODE) @@ -89,4 +97,74 @@ function integrating_factor_solve(eq::LinearODE) v = exp(sympy_integrate(p, eq.t)) end return Symbolics.sympy_simplify((1/v) * ((isequal(eq.q, 0) ? 0 : sympy_integrate(eq.q*v, eq.t)) + eq.C[1])) +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 = @rule cos(~b*t) / sin(~b*t) => ~b + if length(facs) == 1 && !isequal(imag(facs[1]), 0) && get_b(real(facs[1])/imag(facs[1])) !== nothing + 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 + +""" +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 s + p = characteristic_polynomial(eq, s) + Ds = Differential(s) + while isequal(substitute(expand_derivatives((Ds^k)(p)), Dict(s => r)), 0) + k += 1 + end + + return (eq.q*eq.t^k) / (substitute(expand_derivatives((Ds^k)(p)), Dict(s => r))) end \ No newline at end of file diff --git a/test/diffeqs.jl b/test/diffeqs.jl index 154378df8..16bf24e0d 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -50,3 +50,4 @@ C = Symbolics.variables(:C, 1:5) @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)) +@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)) \ No newline at end of file From b4b0f1e908cca532b3afd596e9a73155eae68471 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Mon, 30 Jun 2025 00:20:34 -0400 Subject: [PATCH 10/19] added solving particular solutions when q(t) includes sin or cos and an exponential --- src/diffeqs/diffeqs.jl | 137 +++++++++++++++++++++++++++++++---------- test/diffeqs.jl | 5 +- 2 files changed, 110 insertions(+), 32 deletions(-) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 7123d1e36..6b3a22e80 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -4,21 +4,35 @@ import Symbolics: value, coeff, sympy_integrate struct LinearODE # dⁿx/dtⁿ + pₙ(t)(dⁿ⁻¹x/dtⁿ⁻¹) + ... + p₂(t)(dx/dt) + p₁(t)x = q(t) - x::SymbolicUtils.Symbolic # dependent variable - t::SymbolicUtils.Symbolic # independent variable + x::Num # dependent variable + t::Num # independent variable p::AbstractArray # coefficient functions of t ordered in increasing order (p₁, p₂, ...) - q # right hand side function of t, without any x + q::Any # right hand side function of t, without any x Dt::Differential order::Int C::Vector{Num} # constants - LinearODE(x::Num, t::Num, p, q) = new(value(x), value(t), p, q, Differential(t), length(p), variables(:C, 1:length(p))) - LinearODE(x, t, p, q) = new(x, t, p, q, Differential(t), length(p), variables(:C, 1:length(p))) + function LinearODE(x::Num, t::Num, p, q) + new(value(x), value(t), p, q, Differential(t), + length(p), variables(:C, 1:length(p))) + end + function LinearODE(x, t, p, q) + new(x, t, p, q, Differential(t), length(p), variables(:C, 1:length(p))) + end end -get_expression(eq::LinearODE) = (eq.Dt^eq.order)(eq.x) + sum([(eq.p[n])*(eq.Dt^(n-1))(eq.x) for n = 1:length(eq.p)]) ~ eq.q +function get_expression(eq::LinearODE) + (eq.Dt^eq.order)(eq.x) + sum([(eq.p[n]) * (eq.Dt^(n - 1))(eq.x) for n in 1:length(eq.p)]) ~ eq.q +end -Base.print(io::IO, eq::LinearODE) = print(io, "(D$(eq.t)^$(eq.order))$(eq.x) + " * join(["($(eq.p[length(eq.p)-n]))(D$(eq.t)^$(length(eq.p)-n-1))$(eq.x)" for n = 0:(eq.order-1)], " + ") * " ~ $(eq.q)") +function Base.print(io::IO, eq::LinearODE) + print(io, + "(D$(eq.t)^$(eq.order))$(eq.x) + " * + join( + ["($(eq.p[length(eq.p)-n]))(D$(eq.t)^$(length(eq.p)-n-1))$(eq.x)" + for n in 0:(eq.order - 1)], + " + ") * " ~ $(eq.q)") +end Base.show(io::IO, eq::LinearODE) = print(io, eq) is_homogeneous(eq::LinearODE) = isempty(Symbolics.get_variables(eq.q)) @@ -31,7 +45,7 @@ function characteristic_polynomial(eq::LinearODE, r) @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) + poly += p[i] * r^(i - 1) end return poly @@ -43,10 +57,12 @@ Symbolically solve a linear ODE Cases handled: - ☑ first order - ☑ homogeneous with constant coefficients -- ◩ nonhomogeneous with constant coefficients +- ◩ particular solutions (variation of parameters? undetermined coefficients?) - ☑ ERF + RRF -- ▢ particular solutions (variation of parameters? undetermined coefficients?) + - ☑ complex ERF + RRF to handle sin/cos - ▢ [Differential transform method](https://www.researchgate.net/publication/267767445_A_New_Algorithm_for_Solving_Linear_Ordinary_Differential_Equations) +- ▢ Laplace Transform +- ▢ Expression parsing """ function symbolic_solve_ode(eq::LinearODE) if eq.order == 1 @@ -64,21 +80,25 @@ function symbolic_solve_ode(eq::LinearODE) if rrf !== nothing return const_coeff_solve(to_homogeneous(eq)) + rrf end + rrf_trig = exp_trig_particular_solution(eq) + if rrf_trig !== nothing + return const_coeff_solve(to_homogeneous(eq)) + rrf_trig + end end end function const_coeff_solve(eq::LinearODE) @variables r p = characteristic_polynomial(eq, r) - roots = symbolic_solve(p, r, dropmultiplicity=false) + roots = symbolic_solve(p, r, dropmultiplicity = false) # Handle repeated roots - solutions = exp.(roots*eq.t) - for i in eachindex(solutions)[1:end-1] - j = i+1 + solutions = exp.(roots * eq.t) + for i in eachindex(solutions)[1:(end - 1)] + j = i + 1 while j <= length(solutions) && isequal(solutions[i], solutions[j]) - solutions[j] *= eq.t^(j-i) # multiply by t for each repetition - j+=1 + solutions[j] *= eq.t^(j - i) # multiply by t for each repetition + j += 1 end end @@ -92,11 +112,12 @@ 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) + v = exp(p * eq.t) else v = exp(sympy_integrate(p, eq.t)) end - return Symbolics.sympy_simplify((1/v) * ((isequal(eq.q, 0) ? 0 : sympy_integrate(eq.q*v, eq.t)) + eq.C[1])) + return expand(Symbolics.sympy_simplify((1 / v) * ((isequal(eq.q, 0) ? 0 : + sympy_integrate(eq.q * v, eq.t)) + eq.C[1]))) end """ @@ -104,7 +125,7 @@ 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) @@ -113,19 +134,21 @@ function get_rrf_coeff(q, t) # 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 = @rule cos(~b*t) / sin(~b*t) => ~b - if length(facs) == 1 && !isequal(imag(facs[1]), 0) && get_b(real(facs[1])/imag(facs[1])) !== nothing - 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 + 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 + 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 @@ -133,13 +156,64 @@ function get_rrf_coeff(q, t) der = expand_derivatives(Differential(t)(not_a[1])) r = simplify(der / not_a[1]) - if !isempty(Symbolics.get_variables(r, t)) + if !isempty(Symbolics.get_variables(r, [t])) return nothing end return a, r 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 + +""" +For finding particular solution when q(t) = a*e^(rt)*cos(bt) (or sin(bt)) +""" +function exp_trig_particular_solution(eq::LinearODE) + facs = 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 + + combined_eq = LinearODE(eq.x, eq.t, eq.p, a * exp((r + b * im)eq.t)) + rrf = resonant_response_formula(combined_eq) + + return is_sin ? imag(rrf) : real(rrf) +end + """ Returns a particular solution to a constant coefficient ODE with q(t) = a*e^(rt) @@ -156,7 +230,7 @@ function resonant_response_formula(eq::LinearODE) 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 s @@ -166,5 +240,6 @@ function resonant_response_formula(eq::LinearODE) k += 1 end - return (eq.q*eq.t^k) / (substitute(expand_derivatives((Ds^k)(p)), Dict(s => r))) + return expand(simplify(a * exp(r * eq.t) * eq.t^k / + (substitute(expand_derivatives((Ds^k)(p)), Dict(s => r))))) end \ No newline at end of file diff --git a/test/diffeqs.jl b/test/diffeqs.jl index 16bf24e0d..fff732e8a 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -50,4 +50,7 @@ C = Symbolics.variables(:C, 1:5) @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)) -@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)) \ No newline at end of file +## 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)) From 9292ef4dff3df49bf8ccc9c080799b371e9206b9 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Mon, 30 Jun 2025 09:37:16 -0400 Subject: [PATCH 11/19] added handling of complex characteristic roots --- src/diffeqs/diffeqs.jl | 26 +++++++++++++++++++------- test/diffeqs.jl | 4 ++++ 2 files changed, 23 insertions(+), 7 deletions(-) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 6b3a22e80..7630560ff 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -92,17 +92,28 @@ function const_coeff_solve(eq::LinearODE) p = characteristic_polynomial(eq, r) roots = symbolic_solve(p, r, dropmultiplicity = false) - # Handle repeated roots + # Handle complex + repeated roots solutions = exp.(roots * eq.t) for i in eachindex(solutions)[1:(end - 1)] j = i + 1 - while j <= length(solutions) && isequal(solutions[i], solutions[j]) + + if imag(roots[i]) != 0 && 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^(j - i) # multiply by t for each repetition j += 1 end end - return sum(eq.C .* solutions) + solution = sum(eq.C .* solutions) + if solution isa Complex && isequal(imag(solution), 0) + solution = real(solution) + end + + return solution end """ @@ -117,7 +128,7 @@ function integrating_factor_solve(eq::LinearODE) v = exp(sympy_integrate(p, eq.t)) end return expand(Symbolics.sympy_simplify((1 / v) * ((isequal(eq.q, 0) ? 0 : - sympy_integrate(eq.q * v, eq.t)) + eq.C[1]))) + sympy_integrate(eq.q * v, eq.t)) + eq.C[1]))) end """ @@ -191,13 +202,14 @@ function exp_trig_particular_solution(eq::LinearODE) 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 + 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 && @@ -241,5 +253,5 @@ function resonant_response_formula(eq::LinearODE) end return expand(simplify(a * exp(r * eq.t) * eq.t^k / - (substitute(expand_derivatives((Ds^k)(p)), Dict(s => r))))) + (substitute(expand_derivatives((Ds^k)(p)), Dict(s => r))))) end \ No newline at end of file diff --git a/test/diffeqs.jl b/test/diffeqs.jl index fff732e8a..c481f9ad7 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -50,6 +50,10 @@ C = Symbolics.variables(:C, 1:5) @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 From 1d26ec5ef92ca7e92c1f684f655fee5727959bb1 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Mon, 30 Jun 2025 12:10:43 -0400 Subject: [PATCH 12/19] slight refactoring and improved documentation --- src/Symbolics.jl | 2 +- src/diffeqs/diffeqs.jl | 132 ++++++++++++++++++++++++++++++----------- test/diffeqs.jl | 2 +- 3 files changed, 98 insertions(+), 38 deletions(-) diff --git a/src/Symbolics.jl b/src/Symbolics.jl index d06ca5759..ac146a7c9 100644 --- a/src/Symbolics.jl +++ b/src/Symbolics.jl @@ -221,7 +221,7 @@ export symbolic_solve # Diff Eq Solver include("diffeqs/diffeqs.jl") include("diffeqs/systems.jl") -export symbolic_solve_ode, solve_linear_system +export LinearODE, symbolic_solve_ode, solve_linear_system # Sympy Functions diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 7630560ff..901bf94d0 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -1,45 +1,71 @@ using Symbolics import Symbolics: value, coeff, sympy_integrate +""" +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 - # dⁿx/dtⁿ + pₙ(t)(dⁿ⁻¹x/dtⁿ⁻¹) + ... + p₂(t)(dx/dt) + p₁(t)x = q(t) - - x::Num # dependent variable - t::Num # independent variable - p::AbstractArray # coefficient functions of t ordered in increasing order (p₁, p₂, ...) - q::Any # right hand side function of t, without any x - Dt::Differential - order::Int - C::Vector{Num} # constants - - function LinearODE(x::Num, t::Num, p, q) - new(value(x), value(t), p, q, Differential(t), - length(p), variables(:C, 1:length(p))) - end - function LinearODE(x, t, p, q) - new(x, t, p, q, Differential(t), length(p), variables(:C, 1:length(p))) - end + 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))) end +Dt(eq::LinearODE) = Differential(eq.t) +order(eq::LinearODE) = length(eq.p) + +"""Generates symbolic expression to represent `LinearODE`""" function get_expression(eq::LinearODE) - (eq.Dt^eq.order)(eq.x) + sum([(eq.p[n]) * (eq.Dt^(n - 1))(eq.x) for n in 1:length(eq.p)]) ~ eq.q + (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.print(io::IO, eq::LinearODE) - print(io, - "(D$(eq.t)^$(eq.order))$(eq.x) + " * - join( - ["($(eq.p[length(eq.p)-n]))(D$(eq.t)^$(length(eq.p)-n-1))$(eq.x)" - for n in 0:(eq.order - 1)], - " + ") * " ~ $(eq.q)") +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) +"""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" @@ -63,30 +89,64 @@ Cases handled: - ▢ [Differential transform method](https://www.researchgate.net/publication/267767445_A_New_Algorithm_for_Solving_Linear_Ordinary_Differential_Equations) - ▢ Laplace Transform - ▢ Expression parsing + +Uses methods: [`integrating_factor_solve`](@ref), [`find_homogeneous_solutions`](@ref), [`find_particular_solution`](@ref) + """ function symbolic_solve_ode(eq::LinearODE) - if eq.order == 1 + if order(eq) == 1 return integrating_factor_solve(eq) end + homogeneous_solutions = find_homogeneous_solutions(eq) + if is_homogeneous(eq) - if has_const_coeffs(eq) - return const_coeff_solve(eq) - end + return homogeneous_solutions + end + + return homogeneous_solutions + find_particular_solution(eq) +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 + return sum(find_particular_solution.(terms)) end if has_const_coeffs(eq) rrf = resonant_response_formula(eq) if rrf !== nothing - return const_coeff_solve(to_homogeneous(eq)) + rrf + return rrf end rrf_trig = exp_trig_particular_solution(eq) if rrf_trig !== nothing - return const_coeff_solve(to_homogeneous(eq)) + rrf_trig + return rrf_trig end 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 r p = characteristic_polynomial(eq, r) @@ -96,14 +156,14 @@ function const_coeff_solve(eq::LinearODE) solutions = exp.(roots * eq.t) for i in eachindex(solutions)[1:(end - 1)] j = i + 1 - + if imag(roots[i]) != 0 && 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)) + 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^(j - i) # multiply by t for each repetition + solutions[j] *= eq.t # multiply by t for each repetition j += 1 end end diff --git a/test/diffeqs.jl b/test/diffeqs.jl index c481f9ad7..5c8c48189 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -40,7 +40,7 @@ C = Symbolics.variables(:C, 1:5) ## first order @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.sympy_simplify(C[1]*exp((1//2)t^2 + t) - 1)) +@test isequal(symbolic_solve_ode(LinearODE(x, t, [-(1+t)], 1+t)), 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))) @test_broken isequal(symbolic_solve_ode(LinearODE(x, t, [1], 2sin(t))), C[1]*exp(-t) + sin(t) - cos(t)) From 38c59a3493ee63a3a3be2e0231acd7b0077b29a5 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Mon, 30 Jun 2025 20:28:58 -0400 Subject: [PATCH 13/19] added method of undetermined coefficients --- src/diffeqs/diffeqs.jl | 87 ++++++++++++++++++++++++++++++++++++++++-- test/diffeqs.jl | 6 ++- 2 files changed, 89 insertions(+), 4 deletions(-) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 901bf94d0..0fae30842 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -127,7 +127,7 @@ 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 - return sum(find_particular_solution.(terms)) + return sum(find_particular_solution(LinearODE(eq.x, eq.t, eq.p, term)) for term in terms) end if has_const_coeffs(eq) @@ -140,6 +140,11 @@ function find_particular_solution(eq::LinearODE) return rrf_trig end end + + undetermined_coeff = method_of_undetermined_coefficients(eq) + if undetermined_coeff !== nothing + return undetermined_coeff + end end """ @@ -187,8 +192,12 @@ function integrating_factor_solve(eq::LinearODE) else v = exp(sympy_integrate(p, eq.t)) end - return expand(Symbolics.sympy_simplify((1 / v) * ((isequal(eq.q, 0) ? 0 : - sympy_integrate(eq.q * v, eq.t)) + eq.C[1]))) + solution = (1 / v) * ((isequal(eq.q, 0) ? 0 : sympy_integrate(eq.q * v, eq.t)) + eq.C[1]) + @variables Integral + if !isempty(Symbolics.get_variables(solution, Integral)) + return nothing + end + return expand(Symbolics.sympy_simplify(solution)) end """ @@ -314,4 +323,76 @@ function resonant_response_formula(eq::LinearODE) return expand(simplify(a * exp(r * eq.t) * eq.t^k / (substitute(expand_derivatives((Ds^k)(p)), Dict(s => 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 = Symbolics.degree(eq.q, eq.t) # just a starting point + a = Symbolics.variables(:a, 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 = expand_derivatives(eq_subbed) + try + coeff_solution = symbolic_solve(eq_subbed, length(a) == 1 ? a[1] : a) + catch + coeff_solution = nothing + end + if degree > 0 && coeff_solution !== nothing && !isempty(coeff_solution) + return substitute(form, coeff_solution[1]) + end + + # exponential + @variables a + coeff = get_rrf_coeff(eq.q, eq.t) + if coeff !== nothing + r = coeff[2] + form = a*exp(r*eq.t) + eq_subbed = substitute(get_expression(eq), Dict(eq.x => form)) + eq_subbed = expand_derivatives(eq_subbed) + @show coeff_solution = symbolic_solve(eq_subbed, a) + + if coeff_solution !== nothing && !isempty(coeff_solution) + return substitute(form, coeff_solution[1]) + end + end + + # sin and cos + # this is a hacky way of doing things + @variables a, b + @variables cs, sn + parsed = _parse_trig(factors(eq.q)[end], eq.t) + if parsed !== nothing + ω = parsed[1] + form = a*cos(ω*eq.t) + b*sin(ω*eq.t) + eq_subbed = 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)=>cs, sin(ω*eq.t)=>sn))) + cos_eq = simplify(sum(terms_with(eq_subbed, cs))/cs) + sin_eq = simplify(sum(terms_with(eq_subbed, sn))/sn) + if !isempty(Symbolics.get_variables(cos_eq, [eq.t,sn,cs])) || !isempty(Symbolics.get_variables(sin_eq, [eq.t,sn,cs])) + coeff_solution = nothing + else + coeff_solution = symbolic_solve([cos_eq, sin_eq], [a,b]) + end + + if coeff_solution !== nothing && !isempty(coeff_solution) + return substitute(form, coeff_solution[1]) + end + end +end + +function is_solution(solution, eq) + if solution === nothing + return false + end + + expr = substitute(get_expression(eq), Dict(eq.x => solution)) + @show expr = expand(expand_derivatives(expr.lhs - expr.rhs)) + return isequal(expr, 0) end \ No newline at end of file diff --git a/test/diffeqs.jl b/test/diffeqs.jl index 5c8c48189..db52329f7 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -1,5 +1,5 @@ using Symbolics -using Symbolics: solve_linear_system, LinearODE, is_homogeneous, has_const_coeffs, to_homogeneous, symbolic_solve_ode +using Symbolics: solve_linear_system, LinearODE, is_homogeneous, has_const_coeffs, to_homogeneous, symbolic_solve_ode, find_particular_solution import Groebner, Nemo, SymPy using Test @@ -58,3 +58,7 @@ C = Symbolics.variables(:C, 1:5) @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) \ No newline at end of file From 77f3994bb572cbb040945432a3d6d3714256e16f Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Tue, 1 Jul 2025 02:19:06 -0400 Subject: [PATCH 14/19] added expression parsing --- src/diffeqs/diffeqs.jl | 54 ++++++++++++++++++++++++++++++++++++++++++ test/diffeqs.jl | 7 +++++- 2 files changed, 60 insertions(+), 1 deletion(-) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 0fae30842..49c0c0e06 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -33,6 +33,57 @@ struct LinearODE 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 + epxr = expr.lhs - expr.rhs + end + + expr = expand(simplify(expr)) + terms = Symbolics.terms(epxr) + p::Vector{Num} = [] + q = 0 + order = 0 + for term in terms + if isequal(Symbolics.get_variables(term, [x]), [x]) + facs = factors(term) + deriv = filter(fac -> Symbolics.hasderiv(Symbolics.value(fac)), facs) + p_n = prod(filter(fac -> !Symbolics.hasderiv(Symbolics.value(fac)), facs)) + if isempty(deriv) + p[1] = p_n/x + continue + end + + @assert length(deriv) == 1 "Expected linear term: $term" + n = _get_der_order(deriv[1], x, t) + if n+1 > length(p) + append!(p, zeros(Int, n-length(p) + 1)) + end + p[n + 1] = p_n + order = max(order, n) + + elseif isempty(Symbolics.get_variables(term, [x])) + q -= term + else + @error "Invalid term in LinearODE: $term" + end + end + + # normalize leading coefficient to 1 + leading_coeff = p[order + 1] + p = expand.(p .// leading_coeff) + q = expand(q // leading_coeff) + + new(x, t, p[1:order], q) + end +end + +function _get_der_order(expr, x, t) + if isequal(expr, x) + return 0 + end + + return _get_der_order(substitute(expr, Dict(Differential(t)(x) => x)), x, t) + 1 end Dt(eq::LinearODE) = Differential(eq.t) @@ -53,6 +104,9 @@ 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)) diff --git a/test/diffeqs.jl b/test/diffeqs.jl index db52329f7..0d521a70f 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -61,4 +61,9 @@ C = Symbolics.variables(:C, 1:5) ## 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) \ No newline at end of file +@test isequal(find_particular_solution(LinearODE(x, t, [1, 0], t^2)), t^2 - 2) + +# Parsing +Dt = Differential(t) +@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)) \ No newline at end of file From c1de7c18116299e494bf1338e6729f25059be2dd Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Tue, 1 Jul 2025 12:50:20 -0400 Subject: [PATCH 15/19] updated cases handled --- src/diffeqs/diffeqs.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 49c0c0e06..a2d991dcf 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -137,12 +137,14 @@ Symbolically solve a linear ODE Cases handled: - ☑ first order - ☑ homogeneous with constant coefficients -- ◩ particular solutions (variation of parameters? undetermined coefficients?) +- ◩ particular solutions - ☑ ERF + RRF - ☑ complex ERF + RRF to handle sin/cos + - ☑ method of undetermined coefficients + - ▢ variation of parameters - ▢ [Differential transform method](https://www.researchgate.net/publication/267767445_A_New_Algorithm_for_Solving_Linear_Ordinary_Differential_Equations) - ▢ Laplace Transform -- ▢ Expression parsing +- ☑ Expression parsing Uses methods: [`integrating_factor_solve`](@ref), [`find_homogeneous_solutions`](@ref), [`find_particular_solution`](@ref) From 11bd108e1b6af8bbf383cb1a10c4ef4a6a0a3a92 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Tue, 1 Jul 2025 16:34:58 -0400 Subject: [PATCH 16/19] added IVP --- src/Symbolics.jl | 2 +- src/diffeqs/diffeqs.jl | 38 ++++++++++++++++++++++++++++++++++++++ test/diffeqs.jl | 8 ++++++-- 3 files changed, 45 insertions(+), 3 deletions(-) diff --git a/src/Symbolics.jl b/src/Symbolics.jl index ac146a7c9..c049b395b 100644 --- a/src/Symbolics.jl +++ b/src/Symbolics.jl @@ -221,7 +221,7 @@ export symbolic_solve # Diff Eq Solver include("diffeqs/diffeqs.jl") include("diffeqs/systems.jl") -export LinearODE, symbolic_solve_ode, solve_linear_system +export LinearODE, IVP, symbolic_solve_ode, solve_linear_system, solve_IVP # Sympy Functions diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index a2d991dcf..dd4b26110 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -451,4 +451,42 @@ function is_solution(solution, eq) expr = substitute(get_expression(eq), Dict(eq.x => solution)) @show expr = expand(expand_derivatives(expr.lhs - expr.rhs)) return isequal(expr, 0) +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(substitute(general_solution, symbolic_solve(eqs, ivp.eq.C)[1]))) end \ No newline at end of file diff --git a/test/diffeqs.jl b/test/diffeqs.jl index 0d521a70f..8664f22d8 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -1,5 +1,5 @@ using Symbolics -using Symbolics: solve_linear_system, LinearODE, is_homogeneous, has_const_coeffs, to_homogeneous, symbolic_solve_ode, find_particular_solution +using Symbolics: solve_linear_system, LinearODE, is_homogeneous, has_const_coeffs, to_homogeneous, symbolic_solve_ode, find_particular_solution, IVP, solve_IVP import Groebner, Nemo, SymPy using Test @@ -66,4 +66,8 @@ C = Symbolics.variables(:C, 1:5) # Parsing Dt = Differential(t) @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)) \ No newline at end of file +@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)) \ No newline at end of file From 4c5400eefb6a98045083ac8b85539f25698f09bf Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Thu, 3 Jul 2025 13:15:52 -0400 Subject: [PATCH 17/19] added method to solve Clairaut's equation --- src/diffeqs/diffeqs.jl | 68 ++++++++++++++++++++++++++++++++++++++++-- test/diffeqs.jl | 7 ++++- 2 files changed, 71 insertions(+), 4 deletions(-) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index dd4b26110..f40463a70 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -65,7 +65,8 @@ struct LinearODE elseif isempty(Symbolics.get_variables(term, [x])) q -= term else - @error "Invalid term in LinearODE: $term" + # throw assertion error for invalid term so it can be easily caught + @assert false "Invalid term in LinearODE: $term" end end @@ -135,7 +136,9 @@ end Symbolically solve a linear ODE Cases handled: -- ☑ first order +- ☑ first order linear +- ☑ Clairaut's equation +- ◩ bernoulli equations - ☑ homogeneous with constant coefficients - ◩ particular solutions - ☑ ERF + RRF @@ -163,6 +166,23 @@ function symbolic_solve_ode(eq::LinearODE) return homogeneous_solutions + find_particular_solution(eq) end +function symbolic_solve_ode(expr::Equation, x, t) + if solve_clairaut(expr, x, t) !== nothing + return solve_clairaut(expr, x, t) + end + + try + eq = LinearODE(expr, x, t) + return symbolic_solve_ode(eq) + catch e + if e isa AssertionError + return nothing + else + throw(e) + end + end +end + """ Find homogeneous solutions of linear ODE `eq` with integration constants of `eq.C` @@ -489,4 +509,46 @@ function solve_IVP(ivp::IVP) end return expand(simplify(substitute(general_solution, symbolic_solve(eqs, ivp.eq.C)[1]))) -end \ No newline at end of file +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 + + @variables C + f = substitute(f, Dict(Dt(x) => C)) + if !isempty(Symbolics.get_variables(f, [x])) + return nothing + end + + return C*t + f +end diff --git a/test/diffeqs.jl b/test/diffeqs.jl index 8664f22d8..2eb4bdc76 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -70,4 +70,9 @@ Dt = Differential(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)) \ No newline at end of file +@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 +@variables C +@test isequal(symbolic_solve_ode(x ~ Dt(x)*t - ((Dt(x))^3), x, t), C*t - C^3) +@test isequal(symbolic_solve_ode(x ~ Dt(x)*t + (Dt(x))^2 - sin(Dt(x)) + 2, x, t), C*t + C^2 - sin(C) + 2) \ No newline at end of file From e1a8fb06962ed87c8d5460dc4a4c002053997378 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Thu, 3 Jul 2025 23:48:57 -0400 Subject: [PATCH 18/19] implemented solving bernoulli equations --- src/diffeqs/diffeqs.jl | 153 +++++++++++++++++++++++++++++++++++------ test/diffeqs.jl | 13 ++-- 2 files changed, 141 insertions(+), 25 deletions(-) diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index f40463a70..9dd203519 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -46,11 +46,16 @@ struct LinearODE order = 0 for term in terms if isequal(Symbolics.get_variables(term, [x]), [x]) - facs = factors(term) + facs = _true_factors(term) deriv = filter(fac -> Symbolics.hasderiv(Symbolics.value(fac)), facs) p_n = prod(filter(fac -> !Symbolics.hasderiv(Symbolics.value(fac)), facs)) if isempty(deriv) - p[1] = p_n/x + @assert Symbolics.degree(term, x) == 1 "Expected linear term: $term" + if isempty(p) + p = [p_n/x] + else + p[1] = p_n/x + end continue end @@ -153,22 +158,31 @@ Uses methods: [`integrating_factor_solve`](@ref), [`find_homogeneous_solutions`] """ function symbolic_solve_ode(eq::LinearODE) - if order(eq) == 1 - return integrating_factor_solve(eq) - end - homogeneous_solutions = find_homogeneous_solutions(eq) - - if is_homogeneous(eq) + + if is_homogeneous(eq) && homogeneous_solutions !== nothing return homogeneous_solutions end - - return homogeneous_solutions + find_particular_solution(eq) + + 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 function symbolic_solve_ode(expr::Equation, x, t) - if solve_clairaut(expr, x, t) !== nothing - return solve_clairaut(expr, 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 try @@ -176,6 +190,7 @@ function symbolic_solve_ode(expr::Equation, x, t) return symbolic_solve_ode(eq) catch e if e isa AssertionError + @warn e return nothing else throw(e) @@ -203,7 +218,11 @@ 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 - return sum(find_particular_solution(LinearODE(eq.x, eq.t, eq.p, term)) for term in terms) + 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) @@ -280,7 +299,7 @@ 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) + facs = _true_factors(q) # handle complex r # very convoluted, could probably be improved (possibly by making heavier use of @rule) @@ -338,7 +357,7 @@ end For finding particular solution when q(t) = a*e^(rt)*cos(bt) (or sin(bt)) """ function exp_trig_particular_solution(eq::LinearODE) - facs = factors(eq.q) + facs = _true_factors(eq.q) a = prod(filter(fac -> isempty(Symbolics.get_variables(fac, [eq.t])), facs)) @@ -431,7 +450,7 @@ function method_of_undetermined_coefficients(eq::LinearODE) form = a*exp(r*eq.t) eq_subbed = substitute(get_expression(eq), Dict(eq.x => form)) eq_subbed = expand_derivatives(eq_subbed) - @show coeff_solution = symbolic_solve(eq_subbed, a) + coeff_solution = symbolic_solve(eq_subbed, a) if coeff_solution !== nothing && !isempty(coeff_solution) return substitute(form, coeff_solution[1]) @@ -442,15 +461,15 @@ function method_of_undetermined_coefficients(eq::LinearODE) # this is a hacky way of doing things @variables a, b @variables cs, sn - parsed = _parse_trig(factors(eq.q)[end], eq.t) + parsed = _parse_trig(_true_factors(eq.q)[end], eq.t) if parsed !== nothing ω = parsed[1] form = a*cos(ω*eq.t) + b*sin(ω*eq.t) eq_subbed = 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)=>cs, sin(ω*eq.t)=>sn))) - cos_eq = simplify(sum(terms_with(eq_subbed, cs))/cs) - sin_eq = simplify(sum(terms_with(eq_subbed, sn))/sn) + cos_eq = simplify(sum(filter(term -> !isempty(Symbolics.get_variables(term, cs)), terms(eq_subbed)))/cs) + sin_eq = simplify(sum(filter(term -> !isempty(Symbolics.get_variables(term, sn)), terms(eq_subbed)))/sn) if !isempty(Symbolics.get_variables(cos_eq, [eq.t,sn,cs])) || !isempty(Symbolics.get_variables(sin_eq, [eq.t,sn,cs])) coeff_solution = nothing else @@ -469,7 +488,7 @@ function is_solution(solution, eq) end expr = substitute(get_expression(eq), Dict(eq.x => solution)) - @show expr = expand(expand_derivatives(expr.lhs - expr.rhs)) + expr = expand(expand_derivatives(expr.lhs - expr.rhs)) return isequal(expr, 0) end @@ -544,7 +563,7 @@ function solve_clairaut(expr, x, t) return nothing end - @variables C + C = Symbolics.variable(:C, 1) # constant of integration f = substitute(f, Dict(Dt(x) => C)) if !isempty(Symbolics.get_variables(f, [x])) return nothing @@ -552,3 +571,95 @@ function solve_clairaut(expr, x, t) 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)) + @assert _get_der_order(term//leading_coeff, x, t) == 1 "Expected linear term in $term" + elseif !isempty(Symbolics.get_variables(term, [x])) + facs = _true_factors(term) + x_fac = filter(fac -> !isempty(Symbolics.get_variables(fac, [x])), facs) + @assert length(x_fac) == 1 "Expected linear term in $term" + + 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 +""" +function solve_bernoulli(expr, x, t) + @variables v + eq, n = linearize_bernoulli(expr, x, t, v) + + solution = symbolic_solve_ode(eq) + if solution === nothing + return nothing + end + + return simplify(solution^(1//(1-n))) +end + +""" +Solve Bernoulli equations of the form dx/dt + p(t)x = q(t)x^n with initial condition x(0) = x0 +""" +function solve_bernoulli(expr, x, t, x0) + @variables v + eq, n = linearize_bernoulli(expr, x, t, v) + + v0 = x0^(1-n) # convert initial condition from x(0) to v(0) + + ivp = IVP(eq, [v0]) + solution = solve_IVP(ivp) + if solution === nothing + return nothing + end + + return symbolic_solve(solution ~ x^(1-n), x) +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 + + return true_facs +end \ No newline at end of file diff --git a/test/diffeqs.jl b/test/diffeqs.jl index 2eb4bdc76..cd4700282 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -43,7 +43,7 @@ C = Symbolics.variables(:C, 1:5) @test isequal(symbolic_solve_ode(LinearODE(x, t, [-(1+t)], 1+t)), 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))) -@test_broken isequal(symbolic_solve_ode(LinearODE(x, t, [1], 2sin(t))), C[1]*exp(-t) + sin(t) - cos(t)) +@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(LinearODE(x, t, [1, 2], 0)), C[1]*exp(-t) + C[2]*t*exp(-t)) @@ -73,6 +73,11 @@ Dt = Differential(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 -@variables C -@test isequal(symbolic_solve_ode(x ~ Dt(x)*t - ((Dt(x))^3), x, t), C*t - C^3) -@test isequal(symbolic_solve_ode(x ~ Dt(x)*t + (Dt(x))^2 - sin(Dt(x)) + 2, x, t), C*t + C^2 - sin(C) + 2) \ No newline at end of file + +## 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(x ~ Dt(x)*t + (Dt(x))^2 - sin(Dt(x)) + 2, x, t), C[1]*t + C[1]^2 - sin(C[1]) + 2) + +## Bernoulli equations +@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 isequal(symbolic_solve_ode(Dt(x) + (4//t)*x ~ t^3 * x^2, x, t), 1/(C[1]t^4 - t^4 * log(t))) \ No newline at end of file From 852ae5b6e48ffc983cba03021229926353382256 Mon Sep 17 00:00:00 2001 From: LordOfFrogs Date: Fri, 18 Jul 2025 16:15:04 -0400 Subject: [PATCH 19/19] updated documentation and fixed test imports --- docs/src/manual/ode.md | 15 +++++- src/diffeqs/diffeqs.jl | 111 +++++++++++++++++++++++++++++++++-------- src/diffeqs/systems.jl | 21 +++++++- test/diffeqs.jl | 21 ++++---- 4 files changed, 133 insertions(+), 35 deletions(-) diff --git a/docs/src/manual/ode.md b/docs/src/manual/ode.md index 33e1ba1a1..e120f404a 100644 --- a/docs/src/manual/ode.md +++ b/docs/src/manual/ode.md @@ -24,11 +24,22 @@ 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 +``` + ### SymPy ```@docs diff --git a/src/diffeqs/diffeqs.jl b/src/diffeqs/diffeqs.jl index 9dd203519..f6202b53f 100644 --- a/src/diffeqs/diffeqs.jl +++ b/src/diffeqs/diffeqs.jl @@ -1,6 +1,3 @@ -using Symbolics -import Symbolics: value, coeff, sympy_integrate - """ Represents a linear ordinary differential equation of the form: @@ -138,24 +135,47 @@ function characteristic_polynomial(eq::LinearODE, r) end """ -Symbolically solve a linear ODE + 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 -Cases handled: -- ☑ first order linear -- ☑ Clairaut's equation -- ◩ bernoulli equations -- ☑ homogeneous with constant coefficients -- ◩ particular solutions - - ☑ ERF + RRF - - ☑ complex ERF + RRF to handle sin/cos - - ☑ method of undetermined coefficients - - ▢ variation of parameters -- ▢ [Differential transform method](https://www.researchgate.net/publication/267767445_A_New_Algorithm_for_Solving_Linear_Ordinary_Differential_Equations) -- ▢ Laplace Transform -- ☑ Expression parsing +# Integrating Factor (note that SymPy is required for integration) +julia> symbolic_solve_ode(LinearODE(x, t, [5/t], 7t)) +(C₁ + t^7) / (t^5) -Uses methods: [`integrating_factor_solve`](@ref), [`find_homogeneous_solutions`](@ref), [`find_particular_solution`](@ref) +# 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) @@ -174,6 +194,46 @@ function symbolic_solve_ode(eq::LinearODE) 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 mentioend 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 @@ -592,11 +652,15 @@ function linearize_bernoulli(expr, x, t, v) if Symbolics.hasderiv(Symbolics.value(term)) facs = _true_factors(term) leading_coeff = prod(filter(fac -> !Symbolics.hasderiv(Symbolics.value(fac)), facs)) - @assert _get_der_order(term//leading_coeff, x, t) == 1 "Expected linear term in $term" + if _get_der_order(term//leading_coeff, x, t) != 1 + 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) - @assert length(x_fac) == 1 "Expected linear term in $term" + 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)) @@ -618,7 +682,12 @@ Solve Bernoulli equations of the form dx/dt + p(t)x = q(t)x^n """ function solve_bernoulli(expr, x, t) @variables v - eq, n = linearize_bernoulli(expr, x, t, v) + linearized = linearize_bernoulli(expr, x, t, v) + if linearized === nothing + return nothing + end + + eq, n = linearized solution = symbolic_solve_ode(eq) if solution === nothing diff --git a/src/diffeqs/systems.jl b/src/diffeqs/systems.jl index 521cfdf06..378b3c0e6 100644 --- a/src/diffeqs/systems.jl +++ b/src/diffeqs/systems.jl @@ -6,6 +6,7 @@ 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 @@ -14,7 +15,25 @@ Solve linear continuous dynamical system of differential equations of the form A - `t`: independent variable # Returns -- vector of symbolic solutions +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 diff --git a/test/diffeqs.jl b/test/diffeqs.jl index cd4700282..c37a06558 100644 --- a/test/diffeqs.jl +++ b/test/diffeqs.jl @@ -1,35 +1,34 @@ using Symbolics -using Symbolics: solve_linear_system, LinearODE, is_homogeneous, has_const_coeffs, to_homogeneous, symbolic_solve_ode, find_particular_solution, IVP, solve_IVP -import Groebner, Nemo, SymPy +using Symbolics: solve_linear_system, LinearODE, has_const_coeffs, to_homogeneous, symbolic_solve_ode, find_particular_solution, IVP, solve_IVP +using Groebner, Nemo, SymPy using Test @variables x, y, 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))]) +@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_broken isequal(solve_linear_system([4 -3; 8 -6], [7, 2], t), [18 - 11exp(-2t), 24 - 22exp(-2t)]) - -@test_broken isequal(solve_linear_system([-1 -2; 2 -1], [1, -1], t), [exp(-t)*(cos(2t) + sin(2t)), exp(-t)*(sin(2t) - cos(2t))]) +@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]) # LinearODEs -@test is_homogeneous(LinearODE(x, t, [1, 1], 0)) -@test !is_homogeneous(LinearODE(x, t, [t, 1], t^2)) +@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 is_homogeneous(to_homogeneous(LinearODE(x, t, [t, 1], t^2))) +@test Symbolics.is_homogeneous(to_homogeneous(LinearODE(x, t, [t, 1], t^2))) C = Symbolics.variables(:C, 1:5) @@ -40,7 +39,7 @@ C = Symbolics.variables(:C, 1:5) ## first order @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)), expand(Symbolics.sympy_simplify(C[1]*exp((1//2)t^2 + t) - 1))) +@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))) @test isequal(symbolic_solve_ode(LinearODE(x, t, [1], 2sin(t))), C[1]*exp(-t) + sin(t) - cos(t)) @@ -64,7 +63,7 @@ C = Symbolics.variables(:C, 1:5) @test isequal(find_particular_solution(LinearODE(x, t, [1, 0], t^2)), t^2 - 2) # Parsing -Dt = Differential(t) +Dt = Symbolics.Differential(t) @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))