Ask Your Question
1

Difference between numerical_integral and scipy.integrate.quad

asked 2021-08-21 12:38:32 +0200

ypcman gravatar image

updated 2021-08-21 18:11:04 +0200

slelievre gravatar image

I plot an integral with quad() function in a Python file.

I wanted to plot the same result with SageMath but the values returned are different and I really don't understand why the two graphics are not the same ...

Here is my Sage code :

import scipy.integrate as sp_int

int1(a, j, k, L) =  piecewise([[[0, 1/4], 1 * cos(pi * (j + 1) * a/L) * cos(pi * (k + 1) * a/L)], [(1/4, 3/4), 2 * cos(pi * (j + 1) * a/L) * cos(pi * (k + 1) * a/L)], [[3/4, 1], 1 * cos(pi * (j + 1) * a/L) * cos(pi * (k + 1) * a/L)]])

def int2(a, j, k, L):
    if a < 1 / 4 or a > 3 / 4:
        return cos(pi * (j + 1) * a/L) * cos(pi * (k + 1) * a/L)
    else:
        return 2* cos(pi * (j + 1) * a/L) * cos(pi * (k + 1) * a/L)

def integral1(j):
    k = 2
    L = 0.5
    return numerical_integral(int1, 0, 1, max_points=10, params=[j, k, L])[0]

def integral2(j):
    k = 2
    L = 0.5
    return sp_int.quad(int2, 0, 1, args=(j, k, L))[0]

i1 = plot(integral1, 0, 10, rgbcolor=(0.8,0, 0),
          legend_label='SageMath/numerical_integral')
i2 = plot(integral2, 0, 10, rgbcolor=(0, 0.8, 0),
          legend_label='Python/Quad')
show(i2 + i1)

Plots: numerical_integral vs scipy.integrate.quad

edit retag flag offensive close merge delete

Comments

int1 is piecewise on L, not what you want

FrédéricC gravatar imageFrédéricC ( 2021-08-21 14:33:33 +0200 )edit

Could you plot the functions int1 and int2 in order to see if they are the same?

tolga gravatar imagetolga ( 2021-08-21 15:02:14 +0200 )edit

Thanks for your comments.

  • Frédéric, certes ... : How to piecewise on 'a' ?
  • Tolga : I will try
ypcman gravatar imageypcman ( 2021-08-21 16:23:53 +0200 )edit

If Sage functions can only be piecewise on the last variable, change the order of the variables.

slelievre gravatar imageslelievre ( 2021-08-21 18:14:18 +0200 )edit

2 Answers

Sort by » oldest newest most voted
4

answered 2021-08-22 07:31:18 +0200

dsejas gravatar image

updated 2021-08-22 07:40:02 +0200

Hello, @ypcman! This is actually a good question, and it took me a long time to figure out what the problem was. There are many things happening here, so bear with me. (For the sake of clarity, I will refer to SageMath programming functions, like print(), piecewise(), etc., as "commands", "instructions" or "methods", and mathematical functions, like $f(x)=x^2$, $g(a,x)=a+x$, etc., as "functions".)

Let us consider a simple piecewise function:

$$f(a,x) = \begin{cases}a+x & \text{if $a\in[0,1]$},\cr a-x&\text{if $a\in(-1,0)$}.\end{cases}$$

One would expect to be able to define this by executing

f(a, x) = piecewise([[[0,1], a+x], [(-1,0), a-x]])

Then, by calling f(1,1), you would expect the value $2$. However, if you proceed this way, you will only get the error message

TypeError: __call__() takes from 3 to 4 positional arguments but 5 were given

In simplified terminology, the instruction above expects a mathematical expression on the right-hand-side of the equal sign in order to assign it to the "notation" on the left-hand-side. This way the instruction can define a function. However, the piecewise() method does not define such an expected expression, but a whole function on its own. Then, Sage gets confused by this and tries to "wrap" the piecewise function into another function, which is then called "f". If my words are not clear, try executing the following:

f1(a, x) = a - x
print(f1)
f2 = piecewise([[[0,1], a+x], [(-1,0), a-x]])
print(f2)
f3(a, x) = piecewise([[[0,1], a+x], [(-1,0), a-x]])
print(f3)

The first print instruction will show you (a, x) |--> a - x, which is the usual mathematical notation for functions indicating that a and x are the (independent) variables, and the way this function will be evaluated is by subtracting x from a. The second print will show you

piecewise(a|-->a + x on [0, 1], a|-->a - x on (-1, 0); a)

Although this is a human-readable representation of f2 (the internal representation is more complicated), you will notice that this is a set of different functions on the variable a, each define in its own domain. The third print will show you an even more complicated output:

(a, x) |--> piecewise(a|-->a + x on [0, 1], a|-->a - x on (-1, 0); a)

This means that this is function defined on the pair (a,x) that returns a piecewise function with independent variable a, and with an additional parameter x. The confusing part is that the first and second as and xs are not recognized as the same variables by Sage. Moreover, if you try to obtained the wrapped piecewise function by evaluating something like f3(1,2), Sage calls the subroutine for evaluating piecewise functions (which takes 3 or 4 arguments) with two additional parameters, $a=1$ and $x=2$, thus making a total of 5 or 6 arguments, which produces the error message above.

On the other hand, notice that piecewise functions can only be defined on one independent variable; the remaining variables are assumed to be parameters. How does Sage determine the independent variable of the piecewise function? Apparently it does this by lexicographical precedence. For example, if a and x are present (as in our examples above), it assumes a to be the variable; if c, d and e are present, then it assumes c is the variable. What do you do if you want to specify the variable? You can use the var argument. For example,

f = piecewise([[[0,1], a+x], [(-1,0), a-x]], var=x)

sets x as the independent variable of g2 and lets a as a parameter.

How do you integrate a piecewise function? You can use the integral() method. Unfortunately, I couldn't use the numerical_integral() method on piecewise functions, but that should not be a problem, since the values computed with the integral() method can be used by Sage to plot. In your particular case, you can define

def integral1(j):
    k = 2
    L = 0.5
    return int1.subs(j=j,k=k,L=L).integral(a,0,1)

The int1.subs(j=j,k=k,L=L) replace the values j, k and L in the definition of int1, and then you integrate with respect to a from $0$ to $1$.

I took the liberty of creating a version of your code with all these considerations I have discussed. You can check it here.

I hope this helps!

edit flag offensive delete link more

Comments

Very nice explanation. I agree.

My comment on the existence of a closed-form solution stands, but is irrelevant to what dsejas explained.

Emmanuel Charpentier gravatar imageEmmanuel Charpentier ( 2021-08-22 08:02:43 +0200 )edit

What a great answer ! Deep but clear explanation. I learned a lot from reading you. Thanks.

ypcman gravatar imageypcman ( 2021-08-22 08:55:33 +0200 )edit

@ypcman - wonderful. To mark your question as solved, consider accepting the answer.

slelievre gravatar imageslelievre ( 2021-08-22 11:38:40 +0200 )edit
0

answered 2021-08-22 06:55:57 +0200

Emmanuel Charpentier gravatar image

(Not an answer, but a note needing more space than one allowed in a comment. I'll update this answer if necessary)

What version of Sage do you use ? At least with Sagemath 9.4.rc1, piecewise functions seem problematic :

sage: var("a, j, k, L")
(a, j, k, L)
sage: int1(a, j, k, L) =  piecewise([[[0, 1/4], 1 * cos(pi * (j + 1) * a/L) * cos(pi * (k + 1) * a/L)], [(1/4, 3/4), 2 * cos(pi * (j + 1) * a/L) * cos(pi * (k + 1) * a/L)], [[3/4, 1], 1 * cos(pi * (j + 1) * a/L) * cos(pi * (k + 1) * a/L)]])

So far so good. but :

sage: int1(a, j, k, L)
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-162-c48732677c56> in <module>
----> 1 int1(a, j, k, L)

/usr/local/sage-9/local/lib/python3.9/site-packages/sage/symbolic/expression.pyx in sage.symbolic.expression._eval_on_operands.new_f (build/cythonized/sage/symbolic/expression.cpp:70643)()
  13102         new_args = list(ex._unpack_operands())
  13103         new_args.extend(args)
> 13104         return f(ex, *new_args, **kwds)
  13105     return new_f
  13106 

TypeError: __call__() takes from 3 to 4 positional arguments but 7 were given

Ditto when calling with numerical arguments. This could be a newly-introduced bug...

BTW, I do not see the point of using numerical integration in this specific case : your function $f_1$ can be described as :

$$ f_0 : \left( a, j, k, L \right) \ {\mapsto} \ \cos\left(\frac{\pi a j}{L} + \frac{\pi a}{L}\right) \cos\left(\frac{\pi a k}{L} + \frac{\pi a}{L}\right) $$

$$ f_1 : \begin{cases} f_0\left(a, j, k, L\right) & a < \frac{1}{4} \\ 2\,f_0\left(a, j, k, L\right) & a \geq \frac{1}{4} \wedge a \leq \frac{3}{4} \\ f_0\left(a, j, k, L\right) & a > \frac{3}{4} \end{cases} $$

Your integral $\displaystyle\int_0^1 f_1\left(a, j, k, l\right)\,\mathrm{d}a$ is simply $\displaystyle{\int_0^\frac{1}{4} f_0\left(a, j, k, l\right)\,\mathrm{d}a + \int_\frac{1}{4}^\frac{3}{4} 2\,f_0\left(a, j, k, l\right)\,\mathrm{d}a + \int_\frac{3}{4}^1 f_0\left(a, j, k, l\right)\,\mathrm{d}a}$.

Since

$$ \displaystyle{\int f_0\left(a, j, k, l\right)\,=\,\frac{{\left(L j - L k\right)} \sin\left(\frac{\pi a j}{L} + \frac{\pi a k}{L} + \frac{2 \, \pi a}{L}\right) - {\left(L j + L k + 2 \, L\right)} \sin\left(-\frac{\pi a j}{L} + \frac{\pi a k}{L}\right)}{2 \, {\left(\pi j^{2} - \pi k^{2} + 2 \, \pi j - 2 \, \pi k\right)}}} $$

your (definite) integral, which is a function of j, k and L, has an explicit (closed form) expression (which my laziness leaves to the reader as an exercise...).

HTH,

edit flag offensive delete link more

Comments

I'm using Sage 9.3

ypcman gravatar imageypcman ( 2021-08-22 16:08:46 +0200 )edit

Your Answer

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

Add Answer

Question Tools

1 follower

Stats

Asked: 2021-08-21 12:38:32 +0200

Seen: 152 times

Last updated: Aug 22 '21