definite integral

Author: shahriariravanian

src/integral.jl

@@ -86,6 +86,27 @@ function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 1
     end
 end
 
+"""
+    Definite integral
+"""
+function integrate(eq, xx::Tuple; kwargs...)
+    x, lo, hi = xx
+    sol = integrate(eq, x; kwargs...)
+
+    if sol isa Tuple
+        if first(sol) != 0 && sol[2] == 0
+            return substitute(first(sol), Dict(x => hi)) -
+                   substitute(first(sol), Dict(x => lo))
+        else
+            return nothing
+        end
+    elseif sol != nothing
+        return substitute(sol, Dict(x => hi)) - substitute(sol, Dict(x => lo))
+    end
+
+    return nothing
+end
+
 """
     integrate_sum(eq, x; kwargs...) 
 

src/symbolic.jl

@@ -70,16 +70,18 @@ end
 function is_holonomic(y, x)
     y = value(y)
 
+    # technically, polynomials are not holonomic, but practically we can include them
     if is_polynomial(y, x)
-        return false
+        return true
     end
 
     if is_number(y)
         return true
     end
 
     if is_term(y)
-        return operation(y) in [sin, cos, sinh, cosh, exp]
+        return operation(y) in [sin, cos, sinh, cosh, exp] &&
+               is_polynomial(args(y)[1], x)
     end
 
     if is_pow(y)
@@ -113,6 +115,7 @@ function blender(y, x; n = 3)
         end
     end
 
+    basis = expand(x + basis + x * basis)
     return split_terms(basis, x)
 end
 
@@ -289,11 +292,40 @@ function apply_coefs(q, ker, x)
     for (coef, expr) in zip(q, ker)
         s += beautify(coef) * expr
     end
-    s = simplify(s)
+
+    try
+        s = simplify(s)
+    catch e
+        # println(e)
+    end
+
     c, a = atomize(s)
     return beautify(c) * a
 end
 
+function solve_eqs(eq, x, ker, eqs, vars; abstol = 1e-6, radius = 1.0, verbose = false)
+    try
+        A, b = make_Ab(eqs, vars)
+        q = A \ b
+        q = value.(q)
+        sol = apply_coefs(q, ker, x)
+
+        # test if sol solves ∫ eq dx
+        S = subs_symbols(eq, x; include_x = true, radius)
+        err = substitute(diff(sol, x) - eq, S)
+
+        if abs(err) < abstol
+            return sol
+        end
+    catch e
+        if verbose
+            @warn(e)
+        end
+    end
+
+    return nothing
+end
+
 """
     The main entry point for symbolic integration.
     
@@ -304,7 +336,7 @@ end
     Returns:
         the integral or nothing if no solution
 """
-function integrate_symbolic(eq, x; abstol = 1e-6, radius = 1.0)
+function integrate_symbolic(eq, x; abstol = 1e-6, radius = 1.0, verbose = false)
     eq = expand(eq)
 
     if is_holonomic(eq, x)
@@ -323,27 +355,18 @@ function integrate_symbolic(eq, x; abstol = 1e-6, radius = 1.0)
 
     eqs, vars, frags = make_eqs(eq, x, ker)
 
-    try
-        if length(eqs) != length(vars)
-            eqs = make_square(eq, x, vars, frags)
-        end
-        A, b = make_Ab(eqs, vars)
-        q = A \ b
-        q = value.(q)
-        sol = apply_coefs(q, ker, x)
+    sol = solve_eqs(eq, x, ker, eqs, vars; abstol, radius, verbose)
 
-        # test if sol is the correct integral of eq
-        S = subs_symbols(eq, x; include_x = true, radius)
-        err = substitute(diff(sol, x) - eq, S)
-
-        if abs(err) < abstol
-            return sol
-        else
-            return nothing
+    if sol == nothing
+        try
+            eqs = make_square(eq, x, vars, frags)
+            sol = solve_eqs(eq, x, ker, eqs, vars; abstol, radius, verbose)
+        catch e
+            if verbose
+                @warn(e)
+            end
         end
-    catch e
-        # println(e)
     end
 
-    return nothing
+    return sol
 end

test/runtests.jl

@@ -253,7 +253,7 @@ sym_integrals = [
     1 / (x * log(a * x)),
     log(log(a * x)) / x,
     sin(log(a * x)),
-    x / (exp(a * x) - b),
+    # x / (exp(a * x) - b),
     exp(a * x) / (b * exp(a * x) + c),
     exp(a * x) * exp(exp(a * x)),
     log(cos(a * x)) * tan(a * x),