Error during integration: not of type (UNSIGNED-BYTE 15)

asked 2020-08-08 06:06:42 +0100

lvb884
11 ●1

Hi everyone, I am new to Sage and coding, and I am trying to solve a tricky integral from equation 13 in this physics paper: http://www.nature.com/articles/srep11876 (www.nature.com/articles/srep11876)

Here is the integrand in code form (I would upload a screenshot of the full equation, but I do not have enough karma points):

sage: var ('alpha beta E_o a d_o z_o n u')
(alpha, beta, E_o, a, d_o, z_o, n, u)
sage: j = CC.0
sage: integrand = e**(-j * n * u) / (((a + d_o + 0.5 * z_o * (((16 * pi) / (alpha * beta)) ** (-3))) + 0.5 * z_o * (((16 *
....: pi) / (alpha * beta)) ** (-3)) * sin(u)) ** 3) - 1

The integral is part of a Fourier coefficient, so it is a definite integral with limits -pi and pi. I want to integrate with respect to u and then solve the full equation for beta. All of the other variables are known.

However, when I try to integrate, I get an error:

sage: integrate(integrand, u, -pi, pi)
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-32-84f9192f9549> in <module>()
----> 1 integrate(integrand, u, -pi, pi)

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/misc/functional.py in integral(x, *args, **kwds)
    751     """
    752     if hasattr(x, 'integral'):
--> 753         return x.integral(*args, **kwds)
    754     else:
    755         from sage.symbolic.ring import SR

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.integral (build/cythonized/sage/symbolic/expression.cpp:64474)()
  12434                     R = ring.SR
  12435             return R(integral(f, v, a, b, **kwds))
> 12436         return integral(self, *args, **kwds)
  12437
  12438     integrate = integral

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/integration/integral.py in integrate(expression, v, a, b, algorithm, hold)
    929         return indefinite_integral(expression, v, hold=hold)
    930     else:
--> 931         return definite_integral(expression, v, a, b, hold=hold)
    932
    933

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/function.pyx in sage.symbolic.function.BuiltinFunction.__call__ (build/cythonized/sage/symbolic/function.cpp:12262)()
   1136             res = self._evalf_try_(*args)
   1137             if res is None:
-> 1138                 res = super(BuiltinFunction, self).__call__(
   1139                         *args, coerce=coerce, hold=hold)
   1140

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/function.pyx in sage.symbolic.function.Function.__call__ (build/cythonized/sage/symbolic/function.cpp:6938)()
    595             for i from 0 <= i < len(args):
    596                 vec.push_back((<Expression>args[i])._gobj)
--> 597             res = g_function_evalv(self._serial, vec, hold)
    598         elif self._nargs == 1:
    599             res = g_function_eval1(self._serial,

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/function.pyx in sage.symbolic.function.BuiltinFunction._evalf_or_eval_ (build/cythonized/sage/symbolic/function.cpp:13412)()
   1224         res = self._evalf_try_(*args)
   1225         if res is None:
-> 1226             return self._eval0_(*args)
   1227         else:
   1228             return res

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/integration/integral.py in _eval_(self, f, x, a, b)
    203         for integrator in self.integrators:
    204             try:
--> 205                 A = integrator(*args)
    206             except (NotImplementedError, TypeError, AttributeError):
    207                 pass

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/symbolic/integration/external.py in maxima_integrator(expression, v, a, b)
     44         result = maxima.sr_integral(expression,v)
     45     else:
---> 46         result = maxima.sr_integral(expression, v, a, b)
     47     return result._sage_()
     48

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/interfaces/maxima_lib.py in sr_integral(self, *args)
    788         """
    789         try:
--> 790             return max_to_sr(maxima_eval(([max_integrate],[sr_to_max(SR(a)) for a in args])))
    791         except RuntimeError as error:
    792             s = str(error)

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/libs/ecl.pyx in sage.libs.ecl.EclObject.__call__ (build/cythonized/sage/libs/ecl.c:7794)()
    803         """
    804         lispargs = EclObject(list(args))
--> 805         return ecl_wrap(ecl_safe_apply(self.obj,(<EclObject>lispargs).obj))
    806
    807     def __richcmp__(left, right, int op):

/opt/sagemath-9.1/local/lib/python3.7/site-packages/sage/libs/ecl.pyx in sage.libs.ecl.ecl_safe_apply (build/cythonized/sage/libs/ecl.c:5456)()
    375     if ecl_nvalues > 1:
    376         s = si_coerce_to_base_string(ecl_values(1))
--> 377         raise RuntimeError("ECL says: {}".format(
    378             char_to_str(ecl_base_string_pointer_safe(s))))
    379     else:

RuntimeError: ECL says: 1.8189894035458565e-12 is not of type (UNSIGNED-BYTE 15).

I have no idea where 1.8189894035458565e-12 came from, and I get the same error after assigning values to all of the variables except beta and u.

Has anyone encountered this error before? Is there a better way I should evaluate this?

Comments

A few remarks:

$\LaTeX$ing (i. e. viewing) your code show that your integrand has not the right form.
Do not replace exact constants (1/2, i) by floating-point approximations (0.5,CC.0). This may bias some tests (e. nullity tests) and mislead the algorithms.
To avoid drowning in details, simplify your expression by replacing expression not involving the integration variable by constants.
Explicitely add assumptions : this often allows the integration algorithms to progress...

HTH,

Emmanuel Charpentier ( 2020-08-08 17:03:09 +0100 )edit

add a comment

answered 2020-08-09 00:57:06 +0100

Emmanuel Charpentier
7774 ●6 ●51 ●154

updated 2020-08-09 20:02:26 +0100

EDIT : I checked that none of the integrators available to/via Sage (including Rubi installed in Mathematica) can explicitly compute $\displaystyle\int\frac{e^{iru}}{((A+B\sin u))^3 -1}\,du$ where r is real. However, as already noted, using the hint that $n$ is a integer, solutions exist to compute this integral for any given value of $n$.

Eating my own dogfood :

sage: var("A, B, u", domain="real")
(A, B, u)
sage: var("n", domain="integer")
n
sage: Itg=(cos(n*u)+I*sin(n*u))/((A+B*sin(u))^3-1)

Note that we replace $e^{inu}$ by its de Moivre expansion, to avoid repeating an awkward substitution.

Let's see if this expression has poles in the real segment $(-\pi\ \pi)$:

sage: var("t", domain="real")
t
sage: assume(t>=-1, t<=1)
sage: Itg.denominator().subs(sin(u)==t).solve(t)

So, yes, depending on the values of $A$ and $B$, this expression might have poles (points for which the denominator is null) in the integration path. Integrating there would be a problem... left to the reader as an exercise.

Preliminary trials showed that mots of the integrators available with Sage have problems with this kind of expression :

Maxima (default) will integrate some cases, and fail with nosy questions about some nosy questions about (complicated) relations between A and B, not cured by adding the relevant assumptions.
Sympy can be extremely slow, and may return expressions without direct translation to Sage, thus requiring manual preparation (direct calls to sympy methods), and surgical extraction and conversion of the result.
Giac usually fails to return.
Fricas returns very large Fricas expression, not translatable in Sage. No surgery available here, unless you do all the work in Fricas (not fun).
algorithm="mathematica_free" may or may not work. It often fails for an absence of results in WolframAlpha's answer; this seems to be more or less random (WolframAlpha's availability ?).

However, Mathematica seems to work. The best way seems to embed calls to Mathematica in management code in Sage, returning Sage objects.

A few attempts show that finding primitives for these functions is much easier that finding definite integrals. It seems prudent to first find the primitives, evaluate them and, if no bobbytrap (e. g. poles, bizarre integration results, etc...) is found, find the definite integral by taking appropriate limits.

We may help Mathematica by more or less pre-processing the integrand (e.g. separating the real and imaginary parts) :

sage: %time foo=[{n:v, "C":mathematica.Integrate(cos(v*u).trig_expand()/((A+B*sin
....: (u))^3-1), u), "S":mathematica.Integrate(I*sin(v*u).trig_expand()/((A+B*sin
....: (u))^3-1), u)} for v in (-3..3)]
CPU times: user 82.4 ms, sys: 4.03 ms, total: 86.4 ms
Wall time: 8.7 s
sage: %time bar=[{n:v, "P":mathematica.Integrate((cos(v*u)+I*sin(v*u)).trig_expan
....: d()/((A+B*sin(u))^3-1), u)} for v in (-3..3)]
CPU times: user 47.1 ms, sys: 7.97 ms, total: 55 ms
Wall time: 8.31 s
sage: %time gee=[{n:v, "P":mathematica.Integrate(Itg.subs(n==v).trig_expand(), u)
....: } for v in (-3..3)]
CPU times: user 75.4 ms, sys: 0 ns, total: 75.4 ms
Wall time: 8.22 s

The computing times are close. However, the complexity of the answers are not the same (not illustrated here : some expression don't fit a=on a 19"x19" bedsheet. See for yourself...). Finding the limits is usually slow, and may be problematic in some cases. Finding the limits from separate integrations of real and imaginary parts may be easier...

Substituting the original constants and variables in our "big constants" $A$ and $B$ is left as an exercise to the reader.

So if you can somehow access Mathematica, this problem can be solved by Sage.

EDIT : I have not yet been able to get any of these integrators to create a closed-form formula for an arbitrary $n$ ; all the successes I got where for specific values of $n$.

EDIT : my test case (computing $c_i$ for $i$ integer in $[-3\ 3]$) is stupid, because it doesn't exploit the fact that, all our variables being reals, the values of $c_j$ and $c_{-j}$ should be conjugate...

HTH,

edit flag offensive delete link

add a comment

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Error during integration: not of type (UNSIGNED-BYTE 15)

Comments

1 Answer

Your Answer

Question Tools

Stats

Related questions

Error during integration: not of type (UNSIGNED-BYTE 15) edit

Comments

1 Answer

Your Answer

Question Tools

Stats

Related questions

Error during integration: not of type (UNSIGNED-BYTE 15)