[WIP] Add symbolic integral (#240)

* Basic integral with equality check

* Export integral  and add test

* Fix Typo

* Using Domains for integral

* Fix indentation and isequal

* Add derivative of integral

* Add replaceSym function

* Update expand_derivatives for integral and Leibniz rule

* Modify expand_derivatives

* Remove replaceSym and use substitute instead

* Fix typo

* Fix typo, change return type to Term

* Add tests

* Add tests, use istree

* Fix isequal , add tests

* remove ==true

* Add u^2 test

* fix test

* add istree

* fix whitespace

Co-authored-by: Christopher Rackauckas <accounts@chrisrackauckas.com>

Project.toml

@@ -44,7 +44,6 @@ Reexport = "0.2, 1"
 ReferenceTests = "0.9"
 Requires = "1.1"
 RuntimeGeneratedFunctions = "0.4.3, 0.5"
-SafeTestsets = "0"
 SciMLBase = "1.8"
 Setfield = "0.7"
 SpecialFunctions = "0.7, 0.8, 0.9, 0.10, 1.0"

src/Symbolics.jl

@@ -6,6 +6,9 @@ using LinearAlgebra
 
 using Reexport
 
+using DomainSets
+
+import DomainSets: Domain
 @reexport using SymbolicUtils
 
 import SymbolicUtils: Term, Add, Mul, Pow, Sym, symtype,
@@ -57,6 +60,9 @@ export Differential, expand_derivatives
 
 include("diff.jl")
 
+export Integral
+include("integral.jl")
+
 include("array-lib.jl")
 
 include("linear_algebra.jl")

src/diff.jl

@@ -113,6 +113,26 @@ function expand_derivatives(O::Symbolic, simplify=false; occurances=nothing)
                     return expand_derivatives(operation(arg)(inner), simplify)
                 end
             end
+        elseif isa(operation(arg), Integral)
+            if isa(operation(arg).domain, AbstractInterval)
+                domain = operation(arg).domain
+                a, b = DomainSets.endpoints(domain)
+                c = 0
+                inner_function = expand_derivatives(arguments(arg)[1])
+                if istree(value(a))
+                    t1 = SymbolicUtils.substitute(inner_function, Dict(operation(arg).x => value(a)))
+                    t2 = D(a)
+                    c -= t1*t2
+                end
+                if istree(value(b))
+                    t1 = SymbolicUtils.substitute(inner_function, Dict(operation(arg).x => value(b)))
+                    t2 = D(b)
+                    c += t1*t2
+                end
+                inner = expand_derivatives(D(arguments(arg)[1]))
+                c += operation(arg)(inner)
+                return value(c)
+            end
         end
 
         l = length(arguments(arg))
@@ -151,6 +171,8 @@ function expand_derivatives(O::Symbolic, simplify=false; occurances=nothing)
             x = +((!_iszero(c) ? vcat(c, exprs) : exprs)...)
             return simplify ? SymbolicUtils.simplify(x) : x
         end
+    elseif istree(O) && isa(operation(O), Integral)
+        return operation(O)(expand_derivatives(arguments(O)[1]))
     elseif !hasderiv(O)
         return O
     else

src/integral.jl

@@ -0,0 +1,16 @@
+# basic integral struct with upper bound and lower bound.
+struct Integral{X, T <: Domain} <: Function
+    x
+    domain::T
+    Integral(x,domain) = new{typeof(x),typeof(domain)}(Symbolics.value(x), domain)
+end
+
+(I::Integral)(x) = Term{SymbolicUtils.symtype(x)}(I, [x])
+(I::Integral)(x::Num) = Num(I(Symbolics.value(x)))
+SymbolicUtils.promote_symtype(::Integral, x) = x
+
+function Base.show(io::IO, I::Integral)
+    print(io, "Integral(", I.x, ", ", I.domain, ")")
+end
+
+Base.:(==)(I1::Integral, I2::Integral) = convert(Bool, simplify(isequal(I1.x, I2.x) && isequal(I1.domain, I2.domain)))

test/integral.jl

@@ -0,0 +1,35 @@
+using Symbolics, Test
+using DomainSets
+@variables x y
+
+I1 = Integral(x, ClosedInterval(1, 5))
+I2 = Integral(x, ClosedInterval(1, 5))
+@test I1 == I2
+
+@variables v(..) u(..) x y
+D = Differential(x)
+Dxx = Differential(x)^2
+I = Integral(y, ClosedInterval(1, 5))
+eq = D(I(u(x,y))) ~ 0
+eq_test = I(D(u(x,y)))
+@test isequal(expand_derivatives(eq.lhs), Symbolics.value(eq_test))
+
+eq = Dxx((I(u(x,y)))) + I(D(u(x,y))) ~ 0
+eq_test = I(D(u(x,y))) + I(D(D(u(x,y))))
+@test isequal(expand_derivatives(eq.lhs), Symbolics.value(eq_test))
+
+I = Integral(y, ClosedInterval(1, v(x)))
+eq = D((I(u(x,y)))) ~ 0
+eq_test = I(D(u(x,y))) + D(v(x))*u(x, v(x))
+@test isequal(expand_derivatives(eq.lhs), Symbolics.value(eq_test))
+
+I = Integral(y, ClosedInterval(1, v(x)))
+eq = Dxx((I(u(x,y)))) ~ 0
+eq_test = D(I(D(u(x,y))) + D(v(x))*u(x, v(x)))
+eq_test_ = I(D(D(u(x,y)))) + D(D(v(x)))*u(x, v(x)) + 2D(u(x,v(x)))*D(v(x))
+@test isequal(expand_derivatives(eq.lhs), Symbolics.value(eq_test_))
+@test isequal(expand_derivatives(eq.lhs), expand_derivatives(eq_test))
+
+eq = D((I(u(x,y)^2))) ~ 0
+eq_test = I(2D(u(x,y))*u(x,y)) + D(v(x))*u(x, v(x))^2
+@test isequal(expand_derivatives(eq.lhs), Symbolics.value(eq_test))