diff --git a/Project.toml b/Project.toml index 3cdbdf8..1905d60 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,13 @@ name = "SymbolicNumericIntegration" uuid = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e" authors = ["Shahriar Iravanian "] -version = "1.9.0" +version = "1.9.1" [deps] DataDrivenDiffEq = "2445eb08-9709-466a-b3fc-47e12bd697a2" DataDrivenSparse = "5b588203-7d8b-4fab-a537-c31a7f73f46b" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" diff --git a/README.md b/README.md index a4b4deb..49bd353 100644 --- a/README.md +++ b/README.md @@ -38,10 +38,26 @@ integrate(sin(a * x) * cos(b * x), x; symbolic = true, detailed = false) # Citation -If you use SymbolicNumericIntegration.jl in your research, please cite [this paper](https://arxiv.org/abs/2201.12468): +If you use SymbolicNumericIntegration.jl in your research, please cite the [arxiv](https://arxiv.org/abs/2201.12468) and/or [ISSAC'24](https://dl.acm.org/doi/abs/10.1145/3666000.3669714) papers: ``` +@inproceedings{10.1145/3666000.3669714, +author = {Iravanian, Shahriar and Gowda, Shashi and Rackauckas, Christopher}, +title = {Hybrid Symbolic-Numeric and Numerically-Assisted Symbolic Integration}, +year = {2024}, +isbn = {9798400706967}, +publisher = {Association for Computing Machinery}, +address = {New York, NY, USA}, +url = {https://doi.org/10.1145/3666000.3669714}, +doi = {10.1145/3666000.3669714}, +pages = {410–418}, +numpages = {9}, +keywords = {sparse regression, symbolic integration, symbolic-numeric computation}, +location = {Raleigh, NC, USA}, +series = {ISSAC '24} +} + @article{Iravanian2022, author = {Shahriar Iravanian and Carl Julius Martensen and Alessandro Cheli and Shashi Gowda and Anand Jain and Yingbo Ma and Chris Rackauckas}, doi = {10.48550/arxiv.2201.12468}, diff --git a/src/SymbolicNumericIntegration.jl b/src/SymbolicNumericIntegration.jl index bb77edd..5db3426 100644 --- a/src/SymbolicNumericIntegration.jl +++ b/src/SymbolicNumericIntegration.jl @@ -1,5 +1,6 @@ module SymbolicNumericIntegration +using ForwardDiff using TermInterface: iscall using SymbolicUtils using SymbolicUtils: operation, arguments diff --git a/src/homotopy.jl b/src/homotopy.jl index 91ca12b..ad6882b 100644 --- a/src/homotopy.jl +++ b/src/homotopy.jl @@ -133,9 +133,9 @@ partial_int_rules = [ @rule 𝛷(~x, coth(~u)) => (log(sinh(~u)), ~u) # 1/trigonometric functions @rule 𝛷(~x, 1 / - sin(~u)) => (log(csc(~u) + cot(~u)) + log(sin(~u)), ~u) + sin(~u)) => (log(csc(~u) + cot(~u)) + log(sin(~u)), ~u) @rule 𝛷(~x, 1 / - cos(~u)) => (log(sec(~u) + tan(~u)) + log(cos(~u)), ~u) + cos(~u)) => (log(sec(~u) + tan(~u)) + log(cos(~u)), ~u) @rule 𝛷(~x, 1 / tan(~u)) => (log(sin(~u)) + log(tan(~u)), ~u) @rule 𝛷(~x, 1 / csc(~u)) => (cos(~u) + log(csc(~u)), ~u) @rule 𝛷(~x, 1 / sec(~u)) => (sin(~u) + log(sec(~u)), ~u) @@ -163,7 +163,7 @@ partial_int_rules = [ @rule 𝛷(~x, acoth(~u)) => (~u * acot(~u) + log(~u + 1), ~u) # logarithmic and exponential functions @rule 𝛷(~x, - log(~u)) => ( + log(~u)) => ( ~u + ~u * log(~u) + sum(pow_minus_rule(~u, ~x, -1); init = one(~u)), ~u) @@ -176,18 +176,18 @@ partial_int_rules = [ @rule 𝛷(~x, sqrt(~u)) => ( sum(sqrt_rule(~u, ~x, 0.5); init = one(~u)), ~u) @rule 𝛷(~x, 1 / - sqrt(~u)) => ( + sqrt(~u)) => ( sum(sqrt_rule(~u, ~x, -0.5); init = one(~u)), ~u) # rational functions @rule 𝛷(~x, - 1 / ^(~u::is_univar_poly, ~k::is_pos_int)) => ( + 1 / ^(~u::is_univar_poly, ~k::is_pos_int)) => ( sum(pow_minus_rule(~u, ~x, -~k); init = one(~u)), ~u) @rule 𝛷( - ~x, 1 / ~u::is_univar_poly) => ( + ~x, 1 / ~u::is_univar_poly) => ( sum(pow_minus_rule(~u, ~x, -1); init = one(~u)), ~u) @rule 𝛷(~x, ^(~u, -1)) => (log(~u) + ~u * log(~u), ~u) @@ -245,7 +245,7 @@ function pow_minus_rule(p, x, k; abstol = 1e-8) end function sqrt_rule(p, x, k) - h = Any[p ^ k, p ^ (k + 1)] + h = Any[p^k, p^(k + 1)] Δ = diff(p, x) push!(h, log(Δ / 2 + sqrt(p))) diff --git a/src/integral.jl b/src/integral.jl index 0c04c0e..04f9d04 100644 --- a/src/integral.jl +++ b/src/integral.jl @@ -49,7 +49,7 @@ julia> integrate(x * sin(a * x), (x, 0, 1); symbolic = true, detailed = false) Returns a tuple of (solved, unsolved, err) if `detailed == true`, where - solved: the solved integral + solved: the solved integral unsolved: the residual unsolved portion of the input err: the numerical error in reaching the solution @@ -122,7 +122,7 @@ function integrate(eq, xx::Tuple; kwargs...) sol = integrate(eq, x; kwargs...) if sol isa Tuple - if first(sol) != 0 && sol[2] == 0 + if !isequal(first(sol), 0) && sol[2] == 0 return substitute(first(sol), Dict(x => hi)) - substitute(first(sol), Dict(x => lo)) else @@ -200,7 +200,7 @@ function integrate_sum(eq, x; bypass = false, kwargs...) return expand(solved), unsolved, ε₀ end -# integrate_term is the central part of the univariate integration code that +# integrate_term is the central part of the univariate integration code that # tries different methods to integrate `eq`. # # note: this function will be replaced with solver(prob::Problem) in symbolic.jl @@ -295,8 +295,8 @@ function integrate_term(eq, x; kwargs...) end end -# try_integrate is the main dispatch point to call different sparse solvers. -# It tries to find a linear combination of the basis, whose derivative is +# try_integrate is the main dispatch point to call different sparse solvers. +# It tries to find a linear combination of the basis, whose derivative is # equal to eq function try_integrate(eq, x, basis; plan = default_plan()) if isempty(basis) diff --git a/src/rules.jl b/src/rules.jl index 0f8c60a..57b9f68 100644 --- a/src/rules.jl +++ b/src/rules.jl @@ -110,27 +110,27 @@ trigs_rules = [@rule tan(~x) => sin(~x) / cos(~x) @rule csc(~x) => one(~x) / sin(~x) @rule cot(~x) => cos(~x) / sin(~x) @rule sin(~n::is_int_gt_one * - ~x) => sin((~n - 1) * ~x) * cos(~x) + - cos((~n - 1) * ~x) * sin(~x) + ~x) => sin((~n - 1) * ~x) * cos(~x) + + cos((~n - 1) * ~x) * sin(~x) @rule cos(~n::is_int_gt_one * - ~x) => cos((~n - 1) * ~x) * cos(~x) - - sin((~n - 1) * ~x) * sin(~x) + ~x) => cos((~n - 1) * ~x) * cos(~x) - + sin((~n - 1) * ~x) * sin(~x) @rule tan(~n::is_int_gt_one * - ~x) => (tan((~n - 1) * ~x) + tan(~x)) / - (1 - tan((~n - 1) * ~x) * tan(~x)) + ~x) => (tan((~n - 1) * ~x) + tan(~x)) / + (1 - tan((~n - 1) * ~x) * tan(~x)) @rule csc(~n::is_int_gt_one * - ~x) => sec((~n - 1) * ~x) * sec(~x) * - csc((~n - 1) * ~x) * csc(~x) / - (sec((~n - 1) * ~x) * csc(~x) + - csc((~n - 1) * ~x) * sec(~x)) + ~x) => sec((~n - 1) * ~x) * sec(~x) * + csc((~n - 1) * ~x) * csc(~x) / + (sec((~n - 1) * ~x) * csc(~x) + + csc((~n - 1) * ~x) * sec(~x)) @rule sec(~n::is_int_gt_one * - ~x) => sec((~n - 1) * ~x) * sec(~x) * - csc((~n - 1) * ~x) * csc(~x) / - (csc((~n - 1) * ~x) * csc(~x) - - sec((~n - 1) * ~x) * sec(~x)) + ~x) => sec((~n - 1) * ~x) * sec(~x) * + csc((~n - 1) * ~x) * csc(~x) / + (csc((~n - 1) * ~x) * csc(~x) - + sec((~n - 1) * ~x) * sec(~x)) @rule cot(~n::is_int_gt_one * - ~x) => (cot((~n - 1) * ~x) * cot(~x) - 1) / - (cot((~n - 1) * ~x) + cot(~x)) + ~x) => (cot((~n - 1) * ~x) * cot(~x) - 1) / + (cot((~n - 1) * ~x) + cot(~x)) @rule 1 / sin(~x) => csc(~x) @rule 1 / cos(~x) => sec(~x) @rule 1 / tan(~x) => cot(~x) @@ -147,21 +147,21 @@ trigs_rules = [@rule tan(~x) => sin(~x) / cos(~x) @rule cos(~x + ~y) => cos(~x) * cos(~y) - sin(~x) * sin(~y) @rule tan(~x + ~y) => (tan(~x) + tan(~y)) / (1 - tan(~x) * tan(~y)) @rule csc(~x + - ~y) => sec(~x) * sec(~y) * csc(~x) * csc(~y) / - (sec(~x) * csc(~y) + csc(~x) * sec(~y)) + ~y) => sec(~x) * sec(~y) * csc(~x) * csc(~y) / + (sec(~x) * csc(~y) + csc(~x) * sec(~y)) @rule sec(~x + - ~y) => sec(~x) * sec(~y) * csc(~x) * csc(~y) / - (csc(~x) * csc(~y) - sec(~x) * sec(~y)) + ~y) => sec(~x) * sec(~y) * csc(~x) * csc(~y) / + (csc(~x) * csc(~y) - sec(~x) * sec(~y)) @rule cot(~x + ~y) => (cot(~x) * cot(~y) - 1) / (cot(~x) + cot(~y)) @rule sin(~x - ~y) => sin(~x) * cos(~y) - cos(~x) * sin(~y) @rule cos(~x - ~y) => cos(~x) * cos(~y) + sin(~x) * sin(~y) @rule tan(~x - ~y) => (tan(~x) - tan(~y)) / (1 + tan(~x) * tan(~y)) @rule csc(~x - - ~y) => sec(~x) * sec(~y) * csc(~x) * csc(~y) / - (sec(~x) * csc(~y) - csc(~x) * sec(~y)) + ~y) => sec(~x) * sec(~y) * csc(~x) * csc(~y) / + (sec(~x) * csc(~y) - csc(~x) * sec(~y)) @rule sec(~x - - ~y) => sec(~x) * sec(~y) * csc(~x) * csc(~y) / - (csc(~x) * csc(~y) + sec(~x) * sec(~y)) + ~y) => sec(~x) * sec(~y) * csc(~x) * csc(~y) / + (csc(~x) * csc(~y) + sec(~x) * sec(~y)) @rule cot(~x - ~y) => (cot(~x) * cot(~y) + 1) / (cot(~x) - cot(~y)) # @rule sin(2*~x) => 2*sin(~x)*cos(~x) @@ -250,8 +250,8 @@ rational_rules = [@rule Ω(+(~~xs)) => sum(map(Ω, ~~xs)) decompose_rational(~x), -~k)) @rule Ω(~x / - ^(~y::is_poly, ~k::is_pos_int)) => expand(~x * - ^( + ^(~y::is_poly, ~k::is_pos_int)) => expand(~x * + ^( decompose_rational(~y), ~k)) @rule Ω(~x / ~y::is_poly) => expand(~x * decompose_rational(~y))