special integral functions (Ei/Si/Ci/Li) implemented

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,28 +57,43 @@ 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)
         h₂ = substitute_x(h₂ * ∂h₁^-1, x, sub)
 
         S += expand((1 + h₁) * (1 + h₂))
     end
-    
-    S = simplify(S)
 
     unique([one(x); [equivalent(t, x) for t in terms(S)]])
 end
@@ -91,9 +106,9 @@ function ∂(x)
 end
 
 partial_int_rules = [
-					 # trigonometric functions
-					 @rule 𝛷(sin(~x)) => (cos(~x), ∂(~x))
-                     @rule 𝛷(cos(~x)) => (sin(~x), ∂(~x))
+                     # trigonometric functions
+                     @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))
@@ -119,12 +134,12 @@ partial_int_rules = [
                      @rule 𝛷(^(csch(~x), -1)) => (cosh(~x), ∂(~x))
                      @rule 𝛷(^(sech(~x), -1)) => (sinh(~x), ∂(~x))
                      @rule 𝛷(^(coth(~x), -1)) => (log(cosh(~x)), ∂(~x))
-					 # inverse trigonometric functions
+                     # inverse trigonometric functions
                      @rule 𝛷(asin(~x)) => (~x * asin(~x) + sqrt(1 - ~x * ~x), ∂(~x))
                      @rule 𝛷(acos(~x)) => (~x * acos(~x) + sqrt(1 - ~x * ~x), ∂(~x))
                      @rule 𝛷(atan(~x)) => (~x * atan(~x) + log(~x * ~x + 1), ∂(~x))
-                     @rule 𝛷(acsc(~x)) => (~x * acsc(~x) + atanh(1 - ^(~x,-2)), ∂(~x))
-                     @rule 𝛷(asec(~x)) => (~x * asec(~x) + acosh(~x), ∂(~x))  
+                     @rule 𝛷(acsc(~x)) => (~x * acsc(~x) + atanh(1 - ^(~x, -2)), ∂(~x))
+                     @rule 𝛷(asec(~x)) => (~x * asec(~x) + acosh(~x), ∂(~x))
                      @rule 𝛷(acot(~x)) => (~x * acot(~x) + log(~x * ~x + 1), ∂(~x))
                      # inverse hyperbolic functions
                      @rule 𝛷(asinh(~x)) => (~x * asinh(~x) + sqrt(~x * ~x + 1), ∂(~x))
@@ -134,19 +149,26 @@ partial_int_rules = [
                      @rule 𝛷(asech(~x)) => (asech(~x), ∂(~x))
                      @rule 𝛷(acoth(~x)) => (~x * acot(~x) + log(~x + 1), ∂(~x))
                      # logarithmic and exponential functions
-                     @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)) => (~x + ~x * log(~x) +
+                                          sum(candidate_pow_minus(~x, -1); init = one(~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); init = one(~x)), 1);
+                     @rule 𝛷(^(~x, ~k::is_abs_half)) => (sum(candidate_sqrt(~x, ~k);
+                                                             init = one(~x)), 1);
                      @rule 𝛷(sqrt(~x)) => (sum(candidate_sqrt(~x, 0.5); init = one(~x)), 1);
-                     @rule 𝛷(^(sqrt(~x), -1)) => 𝛷(^(~x, -0.5))                                                      
-					 # rational functions                                                              
-                     @rule 𝛷(^(~x::is_poly, ~k::is_neg)) => (sum(candidate_pow_minus(~x, ~k); init = one(~x)), 1)
+                     @rule 𝛷(^(sqrt(~x), -1)) => 𝛷(^(~x, -0.5))
+                     # rational functions                                                              
+                     @rule 𝛷(^(~x::is_poly, ~k::is_neg)) => (sum(candidate_pow_minus(~x,
+                                                                                     ~k);
+                                                                 init = one(~x)), 1)
                      @rule 𝛷(^(~x, -1)) => (log(~x), ∂(~x))
-                     @rule 𝛷(^(~x, ~k::is_neg_int)) => (sum(^(~x, i) for i in (~k + 1):-1), ∂(~x))
+                     @rule 𝛷(^(~x, ~k::is_neg_int)) => (sum(^(~x, i) for i in (~k + 1):-1),
+                                                        ∂(~x))
                      @rule 𝛷(1 / ~x) => 𝛷(^(~x, -1))
-                     @rule 𝛷(^(~x, ~k)) => (^(~x, ~k + 1), ∂(~x))                     
+                     @rule 𝛷(^(~x, ~k)) => (^(~x, ~k + 1), ∂(~x))
                      @rule 𝛷(1) => (𝑥, 1)
                      @rule 𝛷(~x) => ((~x + ^(~x, 2)), ∂(~x))]
 

src/integral.jl

@@ -33,11 +33,11 @@ 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,
-                   homotopy = true, use_optim=false)
+                   homotopy = true, use_optim = false)
     eq = expand(eq)
     eq = apply_div_rule(eq)
 
@@ -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
@@ -236,13 +239,12 @@ end
 """
 function try_integrate(T, eq, x, basis, radius; kwargs...)
     args = Dict(kwargs)
-	use_optim = args[:use_optim]
+    use_optim = args[:use_optim]
     basis = basis[2:end]    # remove 1 from the beginning
-	
-	if use_optim
-		return solve_optim(T, eq, x, basis, radius; kwargs...)
-	else
-		return solve_sparse(T, eq, x, basis, radius; kwargs...)
-	end
-end
 
+    if use_optim
+        return solve_optim(T, eq, x, basis, radius; kwargs...)
+    else
+        return solve_sparse(T, eq, x, basis, radius; kwargs...)
+    end
+end