I’m trying to perform indefinite and/or definite integrals using Symbolics.Integral(), but I can’t seem to find documentation on how to use it (though, most of Symbolics.jl methods appear undocumented). Currently, through trial and error, I have:
import Symbolics
import IntervalSets
@variables x
f = x^2 + x + 1
domain = IntervalSets.ClosedInterval( 0.0, 1.0 )
var_domain = Symbolics.VarDomainPairing( x, domain )
I = Symbolics.Integral( var_domain )
fI = I(f)
def_integral = fI.val
Which results in:
julia> @variables x
1-element Vector{Num}:
x
julia> f = x^2 + x + 1
1 + x + x^2
julia> domain = IntervalSets.ClosedInterval( 0.0, 1.0 )
0.0 .. 1.0
julia> var_domain = Symbolics.VarDomainPairing( x, domain )
Symbolics.VarDomainPairing(x, 0.0 .. 1.0)
julia> I = Symbolics.Integral( var_domain )
(::Integral{Symbolics.VarDomainPairing}) (generic function with 2 methods)
julia> fI = I(f)
Integral(x, 0.0 .. 1.0)(1 + x + x^2)
julia> def_integral = fI.val
Integral(x, 0.0 .. 1.0)(1 + x + x^2)
julia> @variables x
1-element Vector{Num}:
x
julia> f = x^2 + x + 1
1 + x + x^2
julia> domain = Symbolics.Num.( [ 0, 1 ] )
2-element Vector{Num}:
0
1
julia> fI = SymbolicNumericIntegration.integrate( f, x )[1]
x + (1//3)*(x^3) + (1//2)*(x^2)
julia> f₀ = Symbolics.substitute( fI, Dict( [ x=>domain[1] ] ) )
0//1
julia> f₁ = Symbolics.substitute( fI, Dict( [ x=>domain[2] ] ) )
11//6
julia> def_integral = f₁ - f₀
11//6
My primary questions can thus be summarized as:
How do I evaluate the results of Symbolics.Integral()?
Does SymbolicNumericIntegration.jl support definite integration?
a. Is there a reason SymbolicNumericIntegration doesn’t support definite integrals?
b. Is it reasonable to create a pull-request to SymbolicNumericIntegration for a definite_integral method using symbolic substitution for evaluating definite integrals?
Ok, so I take it then that symbolic integration isn’t implemented? I had assumed it was since .Integral() appeared via Symbolics.<tab><tab> exploring – which is how I do most of my package capabilities exploring.
A related (and perhaps more uplifting) comment, I’m coming back to Julia now to convert some educational material from Matlab and the performance between Julia and Matlab symbolic calculations is ridiculous. Obviously, Matlab supports a broader range of functionality (e.g. multivariate integration), but it appears Julia has a good foundation. Here’s a direct (manual) translation of Matlab code into Julia. Since the ultimate goal is educational, I’m trying to avoid too many performance tricks that might obfuscate the code for non-CS undergrads.
function LagrangeBasis( degree, basis_idx, variate, domain )
ref_domain = PolynomialReferenceDomain( Lagrange )
variate = AffineMapping1D( domain, ref_domain, variate )
node = LinRange( -1.0, 1.0, degree + 1 )
basis_val = 1.0
for n = 1 : degree + 1
if n != basis_idx
basis_val = basis_val * ( (variate - node[n] ) / ( node[basis_idx] - node[n] ) )
end
end
return basis_val
end
function PolynomialReferenceDomain( basis::Basis )
if basis == Bernstein
return [0.0, 1.0]
elseif basis == Chebyshev
return [-1.0, 1.0]
elseif basis == Lagrange
return [-1.0, 1.0]
elseif basis == LagrangeCheby
return [-1.0, 1.0]
elseif basis == Legendre
return [-1.0, 1.0]
elseif basis == Monomial
return [0.0, 1.0]
end
end
function AffineMapping1D( domain, target_domain, x )
x_target = x .- domain[1]
x_target *= ( ( target_domain[2] - target_domain[1] ) / ( domain[2] - domain[1] ) )
x_target += target_domain[1]
end
@enum Basis begin
Bernstein
Chebyshev
Lagrange
LagrangeCheby
Legendre
Monomial
end
Whereas Matlab’s timeit function on its equivalent code (using Symbolic Toolbox calculations) results in 0.0196 sec execution time – meaning Julia’s Symbolics.jl is roughly 500x-600x faster (in this instance) without the user really worrying about performance.
Yeah I think that’s a good way to put it. Symbolics has a good foundation but needs more features to fill it out. That is happening over time, but integration is one piece that’s notably missing.
Does this also affect MethodOfLines.jl directly which tells it doesn’t support IntegralEquations for discretization? given Symbolics.jl is used for defining
If needed, you can always fall back to SymbolicRegression.jl and learn the symbolic integral with a genetic algorithm, by evaluating the derivatives of candidate expressions and seeing if they match the function. (This is actually an even more powerful strategy than available in Matlab/Mathematica/SymPy!)
Then you can just convert the recovered expression to Symbolics.jl afterwards with node_to_symbolic!
Here’s a demo of integrating the function f(x) = x^2 + 2x + 1 + \sin(x) by numerically evaluating its derivatives inside a genetic algorithm:
Code:
using MLJ # for fit/predict
using SymbolicRegression # for SRRegressor
using Zygote # For `enable_autodiff=true`
using SymbolicUtils
f(x) = x^2 + 2x + 1 + sin(x)
X = reshape(0.0:0.01:10.0, :, 1)
y = f.(X[:, 1])
function derivative_loss(tree, dataset::Dataset{T,L}, options, idx) where {T,L}
# Column-major:
X = idx === nothing ? dataset.X : view(dataset.X, :, idx)
∂y = idx === nothing ? dataset.y : view(dataset.y, idx)
ŷ, ∂ŷ, completed = eval_grad_tree_array(tree, X, options; variable=true)
!completed && return L(Inf)
mse = sum(i -> (∂ŷ[i] - ∂y[i])^2, eachindex(∂y)) / length(∂y)
return mse
end
model = SRRegressor(
binary_operators=[+, -, *, /],
unary_operators=[sin, cos, exp],
loss_function=derivative_loss,
enable_autodiff=true,
batching=true,
batch_size=25,
niterations=100,
)
mach = machine(model, X, y)
fit!(mach)
r = report(mach)
eq = r.equations[r.best_idx]
symbolic_eq = node_to_symbolic(eq, model)
While certainly (very) interesting and (very) cool, it’s a little heavy for the purposes of teaching undergraduates L^2 inner products:
Let’s compute the inner product between x^0 and x^1 on the domain [-1,1]. First we will train a symbolic regression model to integrate the integrand with a genetic algorithm – an NP-hard problem…
Meanwhile the students be like:
Why all this machinery when \int x \ \mathrm{dx} is obviously just \frac{x^2}{2} ?
It’s just much cleaner right now in my Matlab-based materials:
No, the ability to represent integrals is separate from the ability to symbolically solve them. Symbolics has the ability to represent integrals right now, which can be used for example to represent integro-differential equations for NeuralPDE and MOL, but it cannot symbolically solve integrals.
You may want to check out SymbolicNumericIntegration.jl which could in theory have a SymbolicRegression.jl backend to how it’s doing the symbolic regression part. It’s a bit smarter than just a direct symbolic regression though, I think you’d enjoy it.
But indeed it’s the kind of algorithm that should exist in a suite of symbolic integration algorithms, and this probably isn’t the the kind of thing you’d tell people to use first. RUBI type methods should come first.
