Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "SymbolicNumericIntegration"
uuid = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
authors = ["Shahriar Iravanian <siravan@svtsim.com>"]
version = "1.1.0"
version = "1.2.0"

[deps]
DataDrivenDiffEq = "2445eb08-9709-466a-b3fc-47e12bd697a2"
Expand Down
17 changes: 10 additions & 7 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -21,16 +21,19 @@ the documentation, which contains the unreleased features.
## Example

```julia
using Symbolics
using SymbolicNumericIntegration
julia> using SymbolicNumericIntegration
julia> using Symbolics

@variables x
julia> @variables x a b

integrate(3x^3 + 2x - 5)
```
julia> integrate(3x^3 + 2x - 5)
(x^2 + (3//4)*(x^4) - (5//1)*x, 0, 0)

```
(x^2 + (3//4)*(x^4) - (5x), 0, 0)
julia> integrate(exp(a * x), x; symbolic = true)
(exp(a*x) / a, 0, 0)

julia> integrate(sin(a * x) * cos(b * x), x; symbolic = true, detailed = false)
(-a*cos(a*x)*cos(b*x) - b*sin(a*x)*sin(b*x)) / (a^2 - (b^2))
```

# Citation
Expand Down
25 changes: 9 additions & 16 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -6,21 +6,14 @@ using Documenter, SymbolicNumericIntegration
include("pages.jl")

makedocs(sitename = "SymbolicNumericIntegration.jl",
authors = "Shahriar Iravanian",
modules = [SymbolicNumericIntegration],
clean = true, doctest = false, linkcheck = true,
strict = [
:doctest,
:linkcheck,
:parse_error,
:example_block,
# Other available options are
# :autodocs_block, :cross_references, :docs_block, :eval_block, :example_block, :footnote, :meta_block, :missing_docs, :setup_block
],
format = Documenter.HTML(analytics = "UA-90474609-3",
assets = ["assets/favicon.ico"],
canonical = "https://docs.sciml.ai/SymbolicNumericIntegration/stable/"),
pages = pages)
authors = "Shahriar Iravanian",
modules = [SymbolicNumericIntegration],
clean = true, doctest = false, linkcheck = true,
warnonly = true,
format = Documenter.HTML(analytics = "UA-90474609-3",
assets = ["assets/favicon.ico"],
canonical = "https://docs.sciml.ai/SymbolicNumericIntegration/stable/"),
pages = pages)

deploydocs(repo = "github.com/SciML/SymbolicNumericIntegration.jl.git";
push_preview = true)
push_preview = true)
113 changes: 41 additions & 72 deletions docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -2,10 +2,13 @@

**SymbolicNumericIntegration.jl** is a hybrid symbolic/numerical integration package that works on the [Julia Symbolics](https://docs.sciml.ai/Symbolics/stable/) expressions.

**SymbolicNumericIntegration.jl** uses a randomized algorithm based on a hybrid of the *method of undetermined coefficients* and *sparse regression* and can solve a large subset of basic standard integrals (polynomials, exponential/logarithmic, trigonometric and hyperbolic, inverse trigonometric and hyperbolic, rational and square root).
The basis of how it works and the theory of integration using the Symbolic-Numeric methods refer to [Basis of Symbolic-Numeric Integration](https://github.com/SciML/SymbolicNumericIntegration.jl/blob/main/docs/theory.ipynb).
**SymbolicNumericIntegration.jl** uses a randomized algorithm based on a hybrid of the *method of undetermined coefficients* and *sparse regression* and can solve a large subset of basic standard integrals (polynomials, exponential/logarithmic, trigonometric and hyperbolic, inverse trigonometric and hyperbolic, rational and square root). The symbolic part of the algorithm is similar (but not identical) to Risch-Bronstein's poor man's integrator and generates a list of ansatzes (candidate terms). The numerical part uses sparse regression adopted from the Sparse identification of nonlinear dynamics (SINDy) algorithm to prune down the ansatzes and find the corresponding coefficients. The basis of how it works and the theory of integration using the Symbolic-Numeric methods refer to [Basis of Symbolic-Numeric Integration](https://github.com/SciML/SymbolicNumericIntegration.jl/blob/main/docs/theory.ipynb).

Function `integrate` returns the integral of a univariate expression with *constant* real or complex coefficients. `integrate` returns a tuple with three values. The first one is the solved integral, the second one is the sum of the unsolved terms, and the third value is the residual error. If `integrate` is successful, the unsolved portion is reported as 0.
[hyint](https://github.com/siravan/hyint) is the python counterpart of **SymbolicNumericIntegration.jl** that works with **sympy** expressions.

Originally, **SymbolicNumericIntegration.jl** expected only univariate expression with *constant* real or complex coefficients as input. As of version 1.2, `integrate` function accepts symbolic constants for a subset of solvable integrals.

`integrate` returns a tuple with three values. The first one is the solved integral, the second one is the sum of the unsolved terms, and the third value is the residual error. If `integrate` is successful, the unsolved portion is reported as 0. If we pass `detailed=false` to `integrate', the output is simplified to only the resulting integrals. In this case, `nothing` is returned if the integral cannot be solved. Note that the simplified output will become the default in a future version.

## Installation

Expand All @@ -19,57 +22,53 @@ Pkg.add("SymbolicNumericIntegration")
Examples:

```julia
using Symbolics
using SymbolicNumericIntegration
julia> using SymbolicNumericIntegration
julia> using Symbolics

@variables x
```
julia> @variables x a b

# if `detailed = true` (default), the output is a tuple of (solution, unsolved portion, err)

```julia
julia> integrate(3x^3 + 2x - 5)
(x^2 + (3//4)*(x^4) - (5x), 0, 0)

julia> integrate((5 + 2x)^-1)
((1//2)*log((5//2) + x), 0, 0.0)

julia> integrate(1 / (6 + x^2 - (5x)))
(log(x - 3) - log(x - 2), 0, 3.339372764128952e-16)
# `detailed = false` simplifies the output to just the resulting integral

julia> integrate(1 / (x^2 - 16))
((1//8)*log(x - 4) - ((1//8)*log(4 + x)), 0, 1.546926788028958e-16)
julia> integrate(x^2 / (16 + x^2); detailed = false)
x + 4atan((-1//4)*x)

julia> integrate(x^2 / (16 + x^2))
(x + 4atan((-1//4)*x), 0, 1.3318788420751984e-16)
julia> integrate(x^2 * log(x); detailed = false)
(1//3)*(x^3)*log(x) - (1//9)*(x^3)

julia> integrate(x^2 / sqrt(4 + x^2))
((1//2)*x*((4 + x^2)^0.5) - ((2//1)*log(x + sqrt(4 + x^2))), 0, 8.702422633074313e-17)
julia> integrate(sec(x) * tan(x); detailed = false)
sec(x)

julia> integrate(x^2 * log(x))
((1//3)*log(x)*(x^3) - ((1//9)*(x^3)), 0, 0)
# Symbolic integration. Here, a is a symbolic constant; therefore, we need
# to explicitly define the independent variable (say, x). Also, we set
# `symbolic = true` to force using the symbolic solver

julia> integrate(x^2 * exp(x))
(2exp(x) + exp(x)*(x^2) - (2x*exp(x)), 0, 0)
julia> integrate(sin(a * x), x; detailed = false, symbolic = true)
(-cos(a*x)) / a

julia> integrate(tan(2x))
((-1//2)*log(cos(2x)), 0, 0)
julia> integrate(x^2 * cos(a * x), x; detailed = false, symbolic = true)
((a^2)*(x^2)*sin(a*x) + 2.0a*x*cos(a*x) - 2.0sin(a*x)) / (a^3)

julia> integrate(sec(x) * tan(x))
(cos(x)^-1, 0, 0)
julia> integrate(log(log(a * x)) / x, x; detailed = false, symbolic = true)
log(a*x)*log(log(a*x)) - log(a*x)

julia> integrate(cosh(2x) * exp(x))
((2//3)*exp(x)*sinh(2x) - ((1//3)*exp(x)*cosh(2x)), 0, 7.073930088880992e-8)
# multiple symbolic constants

julia> integrate(cosh(x) * sin(x))
((1//2)*sin(x)*sinh(x) - ((1//2)*cos(x)*cosh(x)), 0, 4.8956233716268386e-17)
julia> integrate(cosh(a * x) * exp(b * x), x; detailed = false, symbolic = true)
(a*sinh(a*x)*exp(b*x) - b*cosh(a*x)*exp(b*x)) / (a^2 - (b^2))

julia> integrate(cosh(2x) * sin(3x))
(0.153845sinh(2x)*sin(3x) - (0.23077cosh(2x)*cos(3x)), 0, 4.9807620877373405e-6)
# definite integration, passing a tuple of (x, lower bound, higher bound) in the
# second argument

julia> integrate(log(log(x)) * (x^-1))
(log(x)*log(log(x)) - log(x), 0, 0)

julia> integrate(exp(x^2))
(0, exp(x^2), Inf) # as expected!
julia> integrate(x * sin(a * x), (x, 0, 1); symbolic = true, detailed = false)
(sin(a) - a*cos(a)) / (a^2)
```

SymbolicNumericIntegration.jl exports some special integral functions (defined over Complex numbers) and uses them in solving integrals:
Expand All @@ -81,15 +80,15 @@ SymbolicNumericIntegration.jl exports some special integral functions (defined o

For examples:

```
```julia
julia> integrate(exp(x + 1) / (x + 1))
(SymbolicNumericIntegration.Ei(1 + x), 0, 1.1796119636642288e-16)
(SymbolicNumericIntegration.Ei(1 + x), 0, 0.0)

julia> integrate(x * cos(x^2 - 1) / (x^2 - 1))
((1//2)*SymbolicNumericIntegration.Ci(x^2 - 1), 0, 2.7755575615628914e-17)
julia> integrate(x * cos(a*x^2 - 1) / (a*x^2 - 1), x; detailed=false, symbolic=true)
((1//2)*SymbolicNumericIntegration.Ci(a*(x^2) - 1)) / a

julia> integrate(1 / (x*log(log(x))))
(SymbolicNumericIntegration.Li(log(x)), 0, 1.1102230246251565e-16)
julia> integrate(1 / (x*log(log(x))), x; detailed=false, symbolic=true)
SymbolicNumericIntegration.Li(log(x))
```

```@docs
Expand All @@ -102,7 +101,7 @@ integrate(eq, x; kwargs...)

Additionally, 12 test suites from the *Rule-based Integrator* ([Rubi](https://rulebasedintegration.org/)) are included in the `/test` directory. For example, we can test the first one as below. *Axiom* refers to the format of the test files)

```julia
```
using SymbolicNumericIntegration
include("test/axiom.jl") # note, you may need to use the correct path

Expand Down Expand Up @@ -198,33 +197,3 @@ Pkg.status(; mode = PKGMODE_MANIFEST) # hide
```@raw html
</details>
```

```@raw html
You can also download the
<a href="
```

```@eval
using TOML
version = TOML.parse(read("../../Project.toml", String))["version"]
name = TOML.parse(read("../../Project.toml", String))["name"]
link = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
"/assets/Manifest.toml"
```

```@raw html
">manifest</a> file and the
<a href="
```

```@eval
using TOML
version = TOML.parse(read("../../Project.toml", String))["version"]
name = TOML.parse(read("../../Project.toml", String))["name"]
link = "https://github.com/SciML/" * name * ".jl/tree/gh-pages/v" * version *
"/assets/Project.toml"
```

```@raw html
">project</a> file.
```
2 changes: 1 addition & 1 deletion src/SymbolicNumericIntegration.jl
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@ using SymbolicUtils: exprtype, BasicSymbolic
using DataDrivenDiffEq, DataDrivenSparse

include("utils.jl")
include("tree.jl")
include("special.jl")

include("cache.jl")
Expand All @@ -29,6 +30,5 @@ include("integral.jl")
export integrate, generate_basis

include("symbolic.jl")
include("logger.jl")

end # module
66 changes: 0 additions & 66 deletions src/candidates.jl
Original file line number Diff line number Diff line change
Expand Up @@ -64,72 +64,6 @@ function closure(eq, x; max_terms = 50)
unique([[s for s in S]; [s * x for s in S]])
end

function candidate_pow_minus(p, k; abstol = 1e-8)
if isnan(poly_deg(p))
return [p^k, p^(k + 1), log(p)]
end

x = var(p)
d = poly_deg(p)

for j in 1:10 # will try 10 times to find the roots
r, s = find_roots(p, x)
if length(r) + length(s) >= d
break
end
end
s = s[1:2:end]
r = nice_parameter.(r)
s = nice_parameter.(s)

# ∫ 1 / ((x-z₁)(x-z₂)) dx = ... + c₁ * log(x-z₁) + c₂ * log(x-z₂)
q = Any[log(x - u) for u in r]
for i in eachindex(s)
β = s[i]
if abs(imag(β)) > abstol
push!(q, atan((x - real(β)) / imag(β)))
push!(q, (1 + x) * log(x^2 - 2 * real(β) * x + abs2(β)))
else
push!(q, log(x - real(β)))
end
end
q = unique(q)

# return [[p^k, p^(k+1)]; candidates(q₁, x)]
if k ≈ -1
return [[p^k]; q]
else
return [[p^k, p^(k + 1)]; q]
end
end

function candidate_sqrt(p, k)
x = var(p)

h = Any[p^k, p^(k + 1)]

if poly_deg(p) == 2
r, s = find_roots(p, x)
l = leading(p, x)
if length(r) == 2
if sum(r) ≈ 0
r₁ = abs(r[1])
if l > 0
push!(h, acosh(x / r₁))
else
push!(h, asin(x / r₁))
end
end
elseif real(s[1]) ≈ 0
push!(h, asinh(x / imag.(s[1])))
end
end

Δ = expand_derivatives(Differential(x)(p))
push!(h, log(0.5 * Δ + sqrt(p)))
return h
end

###############################################################################

enqueue_expr!(S, q, eq, x) = enqueue_expr!!(S, q, ops(eq)..., x)
Expand Down
Loading