Well, I’m back at it again – needing to integrate Symbolic expressions as part of some educational collateral. I’m wondering what the most robust approach to integrating complicated expressions might be. Consider the following:
import Symbolics
import SymbolicNumericIntegration
using SymbolicNumericIntegration: integrate
import QuadGK
Symbolics.@variables x
m = [ x^n for n in range( 0, 4 ) ]
f = sin( π * x )^2 + cos( x ) - 1
integrand = m * f
julia> display( integrand )
5-element Vector{Symbolics.Num}:
-1 + cos(x) + sin(πx)^2
x*(-1 + cos(x) + sin(πx)^2)
(x^2)*(-1 + cos(x) + sin(πx)^2)
(x^3)*(-1 + cos(x) + sin(πx)^2)
(x^4)*(-1 + cos(x) + sin(πx)^2)
Using integrate from SymbolicsNumericIntegration.jl, which now supports definite integration, sometimes (often) returns the wrong answer:
julia> integrate( integrand[end], (x, 0, 1); detailed=false )
-0.06692331486014069
julia> integrate( integrand[end], (x, 0, 1); detailed=false )
-0.06692331486014069
julia> integrate( integrand[end], (x, 0, 1); detailed=false )
-0.009884419990295636
Cranking down the absolute tolerance doesn’t always help:
## abstol = 1e-12
julia> integrate( integrand[end], (x, 0, 1); detailed=false, abstol=1e-12 )
-0.06692331486014069
julia> integrate( integrand[end], (x, 0, 1); detailed=false, abstol=1e-12 )
-0.009884419990287524
julia> integrate( integrand[end], (x, 0, 1); detailed=false, abstol=1e-12 )
-0.009884419990296808
## abstol = 1e-20
julia> integrate( integrand[end], (x, 0, 1); detailed=false, abstol=1e-20 )
[ Info: Integration is partially successful. Pass `detailed = true` to `integrate` for details
This is especially annoying when broadcasting, as (seemingly) randomly some terms don’t integrate and thus return nothing:
var_tuple = fill( ( x, 0, 1 ), ( 5 ) )
julia> integrate.( integrand, var_tuple; detailed=false )
[ Info: Integration is partially successful. Pass `detailed = true` to `integrate` for details
[ Info: Integration is partially successful. Pass `detailed = true` to `integrate` for details
5-element Vector{Union{Nothing, Float64}}:
0.3414709848078965
0.13177329067603627
nothing
nothing
-0.06692331486014069
julia> integrate.( integrand, var_tuple; detailed=false )
[ Info: Integration is partially successful. Pass `detailed = true` to `integrate` for details
[ Info: Integration is partially successful. Pass `detailed = true` to `integrate` for details
5-element Vector{Union{Nothing, Float64}}:
0.3414709848078965
0.13177329067603627
nothing
-0.07826184164390249
nothing
I’ve generally found trying to tweak kwargs to be an exercise in futility and expression-dependent.
I’ve also tried using QuadGK.jl, resigning myself to using a purely numerical approach, but Symbolics don’t appear compatible with QuadGK?
julia> QuadGK.quadgk( integrand[end], 0, 1 )
ERROR: Sym (x^4)*(-1 + cos(x) + sin(πx)^2) is not callable. Use @syms (x^4)*(-1 + cos(x) + sin(πx)^2)(var1, var2,...) to create it as a callable.
It’s not clear to me how to create a Symbolics callable object in this case (where it’s desired to still be able to do things like transpose and * in other parts of the code):
Symbolics.@syms integrand_5( x )
julia> integrand_5(x) = x^2
ERROR: cannot define function integrand_5; it already has a value
julia> integrand_5
x^2
julia> integrand_5(2)
ERROR: Sym x^2 is not callable. Use @syms x^2(var1, var2,...) to create it as a callable.
Examples of transpose and * not working:
m_1(x) = x^0
m_2(x) = x^1
m_3(x) = x^2
m_4(x) = x^3
m_5(x) = x^4
M = [ m_1, m_2, m_3, m_4, m_5 ]
julia> transpose( M )
1×5 transpose(::Vector{Function}) with eltype Union{}:
Error showing value of type LinearAlgebra.Transpose{Union{}, Vector{Function}}:
ERROR: MethodError: no method matching transpose(::typeof(m_1))
julia> M * f
ERROR: MethodError: no method matching *(::typeof(m_1), ::Symbolics.Num)
julia> reshape( M, 5, 1 ) * reshape( M, 1, 5 )
ERROR: MethodError: no method matching zero(::Type{Function})
Side note #1: the exact integral of integrand[end] can be found in SymPy / Matlab / Mathematica / etc., as:
-\frac{1}{10} + \frac{3}{4\pi^4} - \frac{1}{2\pi^2} + 13\sin(1) - 20\cos(1) \approx -0.009884419990295
Side note #2: the Symbolics.jl documentation claims that it supports “Symbolic Linear Algebra” however this rarely seems to work for me:
julia> var_tuple = fill( ( x, 0, 1 ), ( 5, 5 ) )
julia> A = integrate.( m * transpose( m ), var_tuple; detailed=false )
julia> var_tuple = fill( ( x, 0, 1 ), ( 5 ) )
julia> b = integrate.( m, var_tuple; detailed=false )
julia> A \ b
ERROR: MethodError: no method matching Float64(::Symbolics.Num)
I guess I expected Julia to be a bit more “composable” ala Matlab, isn’t that what multiple dispatch is supposed to provide?
Also, if you use matrix inversion the precision is a bit disappointing:
julia> inv( A ) * b
5-element Vector{Any}:
1.0000000000005969
-5.4569682106375694e-12
-1.8189894035458565e-12
3.637978807091713e-12
-3.637978807091713e-12
I thought the benefit of Symbolic.Num was exactness (to the extent possible)? That’s at least a major benefit of other symbolic engines. Consider Matlab:
>> x = sym( "x", "real" )
>> m(x) = [x^0; x^1; x^2; x^3; x^4]
>> A = int( m * transpose( m ), [0, 1] )
>> b = int( m, [0, 1] )
>> sol = inv( A ) * b
sol =
1
0
0
0
0
>> isAlways( all( sol == sym( [1; 0; 0; 0; 0] ) ) )
ans =
logical
1
