Way to symbolically integrate function fit with table data?

General Usage
1

Is there a way to symbolically integrate a function fit with table data?

…and not have it take more than O(seconds)

// note: already using SymPy.jl, so can’t change that

2

breaking example from github,

  1. lookup table
table_T = [
     0.00
     1.00
     2.50
     3.00
     5.00
     6.50
     7.00
     7.50
     9.00
    10.00
    15.00
    20.00
    25.00
]

table_sigma_v = [
  0.0000e-00
  5.4835e-21
  7.6096e-19
  1.7128e-18
  1.2891e-17
  3.1494e-17
  3.9797e-17
  4.9111e-17
  8.2598e-17
  1.0886e-16
  2.6542e-16
  4.2434e-16
  5.5940e-16
]
  1. functions that should interpolate and integrate table data
using SymPy

@syms T_0

function foo(rho)
  SymPy.interpolate(table_T, table_sigma_v, T_0 / ( 1 - rho ) ^ 2 )  * rho ^ 3 + rho
end

@syms rho_0

function bar()
  b = 3
  c = 1.3

  SymPy.integrate(foo(rho_0), (rho_0, 0.0, 1.0)) * 1e2 * b / c
end
  1. calling bar() returns this error:
ERROR: PyError (ccall(@pysym(:PyObject_Call), PyPtr, (PyPtr, PyPtr, PyPtr), o, arg, C_NULL)) <class 'sympy.polys.polyerrors.PolynomialDivisionFailed'>
PolynomialDivisionFailed()
  File "/Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/Conda/deps/usr/lib/python2.7/site-packages/sympy/integrals/integrals.py", line 1280, in integrate
    risch=risch, manual=manual)

...

  File "/Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/Conda/deps/usr/lib/python2.7/site-packages/sympy/polys/densearith.py", line 1096, in dup_prem
    raise PolynomialDivisionFailed(f, g, K)

Stacktrace:
 [1] pyerr_check at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/PyCall/src/exception.jl:56 [inlined]
 [2] pyerr_check at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/PyCall/src/exception.jl:61 [inlined]
 [3] macro expansion at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/PyCall/src/exception.jl:81 [inlined]
 [4] #_pycall#66(::Array{Any,1}, ::Function, ::PyCall.PyObject, ::SymPy.Sym, ::Vararg{Any,N} where N) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/PyCall/src/PyCall.jl:620
 [5] _pycall(::PyCall.PyObject, ::SymPy.Sym, ::Vararg{Any,N} where N) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/PyCall/src/PyCall.jl:608
 [6] #pycall#70(::Array{Any,1}, ::Function, ::PyCall.PyObject, ::Type{PyCall.PyAny}, ::SymPy.Sym, ::Vararg{Any,N} where N) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/PyCall/src/PyCall.jl:642
 [7] pycall(::PyCall.PyObject, ::Type{PyCall.PyAny}, ::SymPy.Sym, ::Vararg{Any,N} where N) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/PyCall/src/PyCall.jl:642
 [8] #call_sympy_fun#538(::Array{Any,1}, ::Function, ::PyCall.PyObject, ::SymPy.Sym, ::Vararg{Any,N} where N) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/SymPy/src/SymPy.jl:252
 [9] call_sympy_fun(::PyCall.PyObject, ::SymPy.Sym, ::Vararg{Any,N} where N) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/SymPy/src/SymPy.jl:252
 [10] #sympy_meth#539(::Array{Any,1}, ::Function, ::String, ::SymPy.Sym, ::Vararg{Any,N} where N) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/SymPy/src/SymPy.jl:257
 [11] sympy_meth(::String, ::SymPy.Sym, ::Vararg{Any,N} where N) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/SymPy/src/SymPy.jl:256
 [12] #integrate#326(::Array{Any,1}, ::Function, ::SymPy.Sym, ::Tuple{SymPy.Sym,Float64,Float64}, ::Vararg{Tuple{SymPy.Sym,Float64,Float64},N} where N) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/SymPy/src/SymPy.jl:193
 [13] integrate(::SymPy.Sym, ::Tuple{SymPy.Sym,Float64,Float64}) at /Applications/JuliaPro-0.5.1.1.app/Contents/Resources/pkgs-0.5.1.1/lib/v0.6/SymPy/src/SymPy.jl:193
 [14] bar() at ./REPL[7]:5
3

It is unclear what you are trying to do. If the function is interpolated, integration will be numerical.

There are various methods, depending on the kind of interpolation you want to use. For univariate functions, they should be extremely fast. Eg see the example for ApproxFun.jl.

2 Likes
4

As mentioned above, I’m trying to symbolically integrate a function defined by table data

The emphasis is on symbolic (using SymPy.jl). This needs to get piped into a numerical solver at the end that finds T_0

// Some example code highlighting the problem is in the first reply

5

I don’t know anything about SymPy, but you should definitely rescale your data before doing anything with it.
And the error you’re getting seems to indicate something going wrong inside of SymPy (PolynomialDivisionFailed), so you could maybe try a very simple function first?

6

You are starting with numerical data and finishing with a numerical root finder, so there doesn’t seem to be a good reason to use a true symbolic package in the middle.

As has already been stated, just using a suitable polynomial interpolant should be enough?