new expand_basis added

Author: shahriariravanian

Project.toml

@@ -1,7 +1,7 @@
 name = "SymbolicNumericIntegration"
 uuid = "78aadeae-fbc0-11eb-17b6-c7ec0477ba9e"
 authors = ["Shahriar Iravanian <[email protected]>"]
-version = "1.2.0"
+version = "1.2.1"
 
 [deps]
 DataDrivenDiffEq = "2445eb08-9709-466a-b3fc-47e12bd697a2"

src/integral.jl

@@ -86,6 +86,8 @@ function integrate(eq, x = nothing; abstol = 1e-6, num_steps = 2, num_trials = 1
     s, u, ε = integrate_sum(eq, x; bypass, abstol, num_trials, num_steps,
         radius, show_basis, opt, symbolic,
         max_basis, verbose, complex_plane, use_optim)
+        
+    s = beautify(s)
 
     if detailed
         return s, u, ε

src/sparse.jl

@@ -9,7 +9,7 @@ function solve_sparse(eq, x, basis, radius; kwargs...)
     l = find_independent_subset(A; abstol)
     A, basis = A[l, l], basis[l]
 
-    y₁, ε₁ = sparse_fit(A, basis, opt; abstol)
+    y₁, ε₁ = sparse_fit(A, basis; abstol, opt)
     if ε₁ < abstol
         return y₁, ε₁, basis
     end
@@ -30,7 +30,7 @@ function solve_sparse(eq, x, basis, radius; kwargs...)
 
     # moving toward the poles
     modify_basis_matrix!(A, X, eq, x, basis, radius; abstol)
-    y₄, ε₄ = sparse_fit(A, basis, opt; abstol)
+    y₄, ε₄ = sparse_fit(A, basis; abstol, opt)
 
     if ε₄ < abstol || ε₄ < ε₁
         return y₄, ε₄
@@ -49,6 +49,8 @@ end
 
 function init_basis_matrix(eq, x, basis, radius, complex_plane; abstol = 1e-6)
     n = length(basis)
+    eq = value(eq)
+    basis = value.(basis)
 
     # A is an nxn matrix holding the values of the fragments at n random points
     A = zeros(Complex{Float64}, (n, n))
@@ -81,7 +83,7 @@ function init_basis_matrix(eq, x, basis, radius, complex_plane; abstol = 1e-6)
                 end
             end
         catch e
-            println("Error from init_basis_matrix!: ", e)
+              println("Error from init_basis_matrix!: ", e)
         end
         l -= 1
     end
@@ -113,7 +115,7 @@ function DataDrivenSparse.active_set!(idx::BitMatrix, p::SoftThreshold,
     DataDrivenSparse.active_set!(idx, p, abs.(x), λ)
 end
 
-function sparse_fit(A, basis, opt; abstol = 1e-6)
+function sparse_fit(A, basis; abstol = 1e-6, opt = STLSQ(exp.(-10:1:0)))
     n, m = size(A)
 
     try

src/symbolic.jl

@@ -31,7 +31,7 @@ function equiv(y, x)
     if is_add(y)
         return sum(equiv(u, x) for u in args(y))
     elseif is_mul(y)
-        return prod(isdependent(u, x) ? u : 1 for u in args(y))
+        return prod(isdependent(u, x) ? u : 1 for u in args(y))        
     elseif is_div(y)
         return expand_fraction(y, x)
     elseif is_number(y)
@@ -132,21 +132,22 @@ end
     Returns:
         a dict of sym : value pairs
 """
-function subs_symbols(eq, x; include_x = false, radius = 5.0)
+function subs_symbols(eq, x; include_x = false, radius = 5.0, as_complex=true)
     S = Dict()
     for v in get_variables(value(eq))
         if !isequal(v, x)
-            S[v] = randn()
+            S[v] = as_complex ? Complex(randn()) : randn()
         end
     end
 
     if include_x
-        S[x] = randn()
+        S[x] = as_complex ? Complex(randn()) : randn()
     end
 
     return S
 end
 
+
 """
     Splits terms into a part dependent on x and a part constant w.r.t. x    
     For example, `atomize(a*sin(b*x), x)` is `(a, sin(b*x))`)
@@ -209,8 +210,8 @@ end
 function make_eqs(eq, x, basis)
     frags = Dict()
 
-    for d in terms(eq)
-        c, a = atomize(d, x)
+    for term in terms(eq)
+        c, a = atomize(term, x)
         frags[a] = c
     end
 
@@ -219,8 +220,8 @@ function make_eqs(eq, x, basis)
     for (i, b) in enumerate(basis)
         db = expand_fraction(diff(b, x), x)
 
-        for d in terms(db)
-            c, a = atomize(d, x)
+        for term in terms(db)
+            c, a = atomize(term, x)
             frags[a] = get(frags, a, 0) + c * θ[i]
         end
 
@@ -277,7 +278,14 @@ function make_square(eq, x, vars, frags)
     for i in 1:n
         S = subs_symbols(eq, x; include_x = true)
         q = sum(substitute(k, S) * v for (k, v) in frags)
-        push!(eqs, q ~ 0)
+        Q = (q ~ 0)
+        if Q isa Array
+            # Q returns a complex array
+            # a different pathway is needed here!
+            return nothing
+        else
+            push!(eqs, Q)
+        end
     end
 
     return eqs
@@ -303,29 +311,6 @@ function apply_coefs(q, ker, x)
     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.
     
@@ -336,31 +321,49 @@ end
     Returns:
         the integral or nothing if no solution
 """
-function integrate_symbolic(eq, x; abstol = 1e-6, radius = 1.0, verbose = false)
+function integrate_symbolic(eq, x; abstol = 1e-6, radius = 1.0, verbose = false, num_steps=2)
     eq = expand(eq)
+    coef, eq = atomize(eq, x)
 
     if is_holonomic(eq, x)
         basis = blender(eq, x)
     else
         basis = generate_homotopy(eq, x)
     end
+    
+    for k = 1:num_steps    
+        sol = try_symbolic(eq, x, basis; abstol, radius, verbose)
+    
+        if sol != nothing
+            return coef * sol            
+        end        
+       
+        if k < num_steps 
+            basis = expand_basis_symbolic(basis, x) 
+        end
+    end
+    
+    return nothing
+end
 
+
+function try_symbolic(eq, x, basis; abstol = 1e-6, radius = 1.0, verbose = false)
     ker = best_hints(eq, x, basis)
 
     if ker == nothing
         return nothing
     end
 
     ker = [atomize(y, x)[2] for y in ker]
-
     eqs, vars, frags = make_eqs(eq, x, ker)
-
     sol = solve_eqs(eq, x, ker, eqs, vars; abstol, radius, verbose)
 
     if sol == nothing
         try
             eqs = make_square(eq, x, vars, frags)
-            sol = solve_eqs(eq, x, ker, eqs, vars; abstol, radius, verbose)
+            if eqs != nothing
+                sol = solve_eqs(eq, x, ker, eqs, vars; abstol, radius, verbose=false)
+            end
         catch e
             if verbose
                 @warn(e)
@@ -370,3 +373,36 @@ function integrate_symbolic(eq, x; abstol = 1e-6, radius = 1.0, verbose = false)
 
     return sol
 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
+
+
+function expand_basis_symbolic(basis, x)
+    b = sum(basis)
+    basis = split_terms(expand((1+x)*(b + diff(b, x))), x)    
+    
+    return basis
+end
+