merging

Author: shahriariravanian

Project.toml

@@ -1,30 +1,28 @@
 name = "SymbolicNumericIntegration"
 uuid = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
 authors = ["Shahriar Iravanian <siravan@svtsim.com>"]
-version = "1.0.2"
+version = "1.1.0"
 
 [deps]
 DataDrivenDiffEq = "2445eb08-9709-466a-b3fc-47e12bd697a2"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 DataDrivenDiffEq = "0.8"
 DataStructures = "0.18"
-PyCall = "1"
-SymPy = "1"
+Optim = "1"
 SymbolicUtils = "0.19"
 Symbolics = "4"
 julia = "1.6"
 
 [extras]
-PyCall = "438e738f-606a-5dbb-bf0a-cddfbfd45ab0"
-SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["PyCall", "SymPy", "Test"]
+test = ["Test"]

docs/src/index.md

@@ -70,6 +70,26 @@ julia> integrate(exp(x^2))
 (0, exp(x^2), Inf)    # as expected!
 ```
 
+SymbolicNumericIntegration.jl exports some special integral functions (defined over Complex numbers) and uses them in solving integrals:
+
+  - `Ei`: exponential integral (define as ∫ exp(x) / x dx)
+  - `Si`: sine integral (define as ∫ sin(x) / x dx)
+  - `Ci`: cosine integral (define as ∫ cos(x) / x dx)
+  - `Li`: logarithmic integral (define as ∫ 1 / log(x) dx)
+
+For examples:
+
+```
+julia> integrate(exp(x + 1) / (x + 1))
+(SymbolicNumericIntegration.Ei(1 + x), 0, 1.1796119636642288e-16)
+
+julia> integrate(x * cos(x^2 - 1) / (x^2 - 1))
+((1//2)*SymbolicNumericIntegration.Ci(x^2 - 1), 0, 2.7755575615628914e-17)
+
+julia> integrate(1 / (x*log(log(x))))
+(SymbolicNumericIntegration.Li(log(x)), 0, 1.1102230246251565e-16)
+```
+
 `integrate` has the form `integrate(y; kw...)` or `integrate(y, x; kw...)`, where `y` is the integrand and the optional `x` is the variable of integration. The keyword parameters are:
 
   - `abstol` (default `1e-6`): the error tolerance to accept a solution.

src/SymbolicNumericIntegration.jl

@@ -9,13 +9,17 @@ using SymbolicUtils.Rewriters
 using DataDrivenDiffEq
 
 include("utils.jl")
+include("special.jl")
+
 include("cache.jl")
 include("roots.jl")
 
 include("rules.jl")
 include("candidates.jl")
 include("homotopy.jl")
 
+export Ei, Si, Ci, Li
+
 include("numeric_utils.jl")
 include("sparse.jl")
 include("optim.jl")

src/cache.jl

@@ -41,7 +41,7 @@ function deriv_fun!(c::ExprCache, xs...)
 end
 
 function deriv_fun!(c, xs...)
-    return build_function(deriv!(c, xs), xs...; expression = false)
+    return build_function(deriv!(c, xs...), xs...; expression = false)
 end
 
 Base.show(io::IO, c::ExprCache) = show(expr(c))

src/candidates.jl

@@ -1,10 +1,15 @@
 using DataStructures
 
 # this is the main heurisctic used to find the test fragments
-function generate_basis(eq, x; homotopy = true)
-    # if homotopy return generate_homotopy(eq, x) end
+function generate_basis(eq, x, try_kernel = true; homotopy = true)
+    if homotopy && !try_kernel
+        S = sum(generate_homotopy(expr(eq), x))
+        return cache.(unique([one(x); [equivalent(t, x) for t in terms(S)]]))
+    end
+
+    S = 0
     eq = expand(expr(eq))
-    S = 0 #Set{Any}()
+
     for t in terms(eq)
         q = equivalent(t, x)
         f = kernel(q)

src/homotopy.jl

@@ -57,19 +57,36 @@ end
 
 ##############################################################################
 
+Symbolics.@register_symbolic Ei(z)
+Symbolics.@register_symbolic Si(z)
+Symbolics.@register_symbolic Ci(z)
+Symbolics.@register_symbolic Li(z)
+
+Symbolics.derivative(::typeof(Ei), args::NTuple{1, Any}, ::Val{1}) = exp(args[1]) / args[1]
+Symbolics.derivative(::typeof(Si), args::NTuple{1, Any}, ::Val{1}) = sin(args[1]) / args[1]
+Symbolics.derivative(::typeof(Ci), args::NTuple{1, Any}, ::Val{1}) = cos(args[1]) / args[1]
+Symbolics.derivative(::typeof(Li), args::NTuple{1, Any}, ::Val{1}) = 1 / log(args[1])
+
+@syms si(𝑥) ci(𝑥) ei(𝑥) li(𝑥)
+
+##############################################################################
+
 function substitute_x(eq, x, sub)
     eq = substitute(eq, sub)
     substitute(eq, Dict(𝑥 => x))
 end
 
 function generate_homotopy(eq, x)
     eq = eq isa Num ? eq.val : eq
+    x = x isa Num ? x.val : x
+
     q, sub = transform(eq, x)
     S = 0
 
     for i in 1:length(sub)
         μ = u[i]
         h₁, ∂h₁ = apply_partial_int_rules(sub[μ])
+        h₁ = substitute(h₁, Dict(si => Si, ci => Ci, ei => Ei, li => Li))
         h₂ = expand_derivatives(Differential(μ)(q))
 
         h₁ = substitute_x(h₁, x, sub)
@@ -78,8 +95,6 @@ function generate_homotopy(eq, x)
         S += expand((1 + h₁) * (1 + h₂))
     end
 
-    S = simplify(S)
-
     unique([one(x); [equivalent(t, x) for t in terms(S)]])
 end
 
@@ -92,8 +107,8 @@ end
 
 partial_int_rules = [
                      # trigonometric functions
-                     @rule 𝛷(sin(~x)) => (cos(~x), ∂(~x))
-                     @rule 𝛷(cos(~x)) => (sin(~x), ∂(~x))
+                     @rule 𝛷(sin(~x)) => (cos(~x) + si(~x), ∂(~x))
+                     @rule 𝛷(cos(~x)) => (sin(~x) + ci(~x), ∂(~x))
                      @rule 𝛷(tan(~x)) => (log(cos(~x)), ∂(~x))
                      @rule 𝛷(csc(~x)) => (log(csc(~x) + cot(~x)), ∂(~x))
                      @rule 𝛷(sec(~x)) => (log(sec(~x) + tan(~x)), ∂(~x))
@@ -137,7 +152,8 @@ partial_int_rules = [
                      @rule 𝛷(log(~x)) => (~x + ~x * log(~x) +
                                           sum(candidate_pow_minus(~x, -1); init = one(~x)),
                                           ∂(~x))
-                     @rule 𝛷(exp(~x)) => (exp(~x), ∂(~x))
+                     @rule 𝛷(^(log(~x), -1)) => (log(log(~x)) + li(~x), ∂(~x))
+                     @rule 𝛷(exp(~x)) => (exp(~x) + ei(~x), ∂(~x))
                      @rule 𝛷(^(exp(~x), ~k::is_neg)) => (^(exp(-~x), -~k), ∂(~x))
                      # square-root functions
                      @rule 𝛷(^(~x, ~k::is_abs_half)) => (sum(candidate_sqrt(~x, ~k);

src/integral.jl

@@ -33,7 +33,7 @@ Base.signbit(x::SymbolicUtils.Sym{Number, Nothing}) = false
     solved is the solved integral and unsolved is the residual unsolved portion of the input
     err is the numerical error in reaching the solution
 """
-function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 5,
+function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 10,
                    radius = 1.0,
                    show_basis = false, opt = STLSQ(exp.(-10:1:0)), bypass = false,
                    symbolic = true, max_basis = 100, verbose = false, complex_plane = true,
@@ -150,13 +150,14 @@ function integrate_term(eq, x, l; kwargs...)
     end
 
     eq = cache(eq)
-    basis = generate_basis(eq, x; homotopy)
+    basis1 = generate_basis(eq, x, true; homotopy)
+    basis2 = generate_basis(eq, x, false; homotopy)
 
     if show_basis
-        inform(l, "Generating basis (|β| = $(length(basis)))", basis)
+        inform(l, "Generating basis (|β| = $(length(basis)))", basis1)
     end
 
-    if length(basis) > max_basis
+    if length(basis1) > max_basis
         result(l, "|β| = $(length(basis)) is too large")
         return 0, expr(eq), Inf
     end
@@ -170,12 +171,12 @@ function integrate_term(eq, x, l; kwargs...)
     yᵣ = 0
 
     for i in 1:num_steps
-        if length(basis) > max_basis
+        if length(basis1) > max_basis
             break
         end
 
         if symbolic
-            y, ϵ = try_symbolic(Float64, expr(eq), x, expr.(basis), deriv!.(basis, x);
+            y, ϵ = try_symbolic(Float64, expr(eq), x, expr.(basis1), deriv!.(basis1, x);
                                 kwargs...)
 
             if !isequal(y, 0) && accept_solution(eq, x, y, radius) < abstol
@@ -187,6 +188,7 @@ function integrate_term(eq, x, l; kwargs...)
         end
 
         for j in 1:num_trials
+            basis = isodd(j) ? basis1 : basis2
             r = radius #*sqrt(2)^j
             y, ϵ = try_integrate(Float64, eq, x, basis, r; kwargs...)
 
@@ -203,10 +205,11 @@ function integrate_term(eq, x, l; kwargs...)
         inform(l, "Failed numeric")
 
         if i < num_steps
-            basis = expand_basis(basis, x)
+            basis1 = expand_basis(basis1, x)
+            basis2 = expand_basis(basis2, x)
 
             if show_basis
-                inform(l, "Expanding the basis (|β| = $(length(basis)))", basis)
+                inform(l, "Expanding the basis (|β| = $(length(basis)))", basis1)
             elseif verbose
                 inform(l, "Expanding the basis (|β| = $(length(basis)))")
             end

src/special.jl

@@ -0,0 +1,71 @@
+using SpecialFunctions
+
+# Adopted from mpmath (https://mpmath.org/), which is a Python library for real and 
+# complex floating-point arithmetic with arbitrary precision
+
+E1(z) = Complex(expint(z))
+
+function Ei(z)
+    if z isa Real
+        return Complex(-expint(Complex(-z)))
+    end
+
+    v = -expint(-z)
+
+    if imag(z) > 0
+        v += Complex(0, π)
+    else
+        v -= Complex(0, π)
+    end
+
+    return v
+end
+
+function Ci(z)
+    if z isa Real
+        return Complex(cosint(z))
+    end
+
+    v = (Ei(z * im) + Ei(-z * im)) / 2
+
+    if real(z) ≈ 0
+        if imag(z) >= 0
+            v += Complex(0, π / 2)
+        else
+            v -= Complex(0, π / 2)
+        end
+    end
+
+    if real(z) < 0
+        if imag(z) >= 0
+            v += Complex(0, π)
+        else
+            v -= Complex(0, π)
+        end
+    end
+    return v
+end
+
+function Si(z)
+    if z isa Real
+        return Complex(sinint(z))
+    end
+
+    v = (Ei(z * im) - Ei(-z * im)) / 2im
+
+    if real(z) ≈ 0
+        return v
+    end
+
+    if real(z) > 0
+        v -= π / 2
+    else
+        v += π / 2
+    end
+
+    return v
+end
+
+function Li(z)
+    return Ei(log(z)) - Ei(log(2.0))
+end

test/runtests.jl

@@ -186,6 +186,17 @@ basic_integrals = [
     x / (x^4 - 4),
     x^3 / (x^4 + 1),
     1 / (x^4 + 1),
+    # exponential/trigonometric/logarithmic integral functions
+    exp(2x) / x,
+    exp(x + 1) / (x + 1),
+    x * exp(2x^2 + 1) / (2x^2 + 1),
+    sin(3x) / x,
+    sin(x + 1) / (x + 1),
+    cos(5x) / x,
+    x * cos(x^2 - 1) / (x^2 - 1),
+    1 / log(3x - 1),
+    1 / (x * log(log(x))),
+    x / log(x^2),
     # bypass = true
     β,      # turn of bypass = true
     (log(x - 1) + (x - 1)^-1) * log(x),