Difference between numerical_integral and scipy.integrate.quad

savecancel

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!