quad

In [7]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

from scipy.integrate   import quad       ### THIS IS NEW

Integration

  • To numerically integrate a function: $I=\int_{x_{lo}}^{x_{hi}}f(x)dx$
    1. Define the function
    2. Set the integration bounds
    3. Call quad(f,xlo,xhi)
      1. quad returns two values, the result and the error.

Example

$$\int_0^13x^2+1dx = \left[x^3+x\right]_0^1 = 2$$
In [8]:
def f(x) :
    return 3.0*x*x + 1.0

I, err = quad(f, 0, 1)    

print("I     = ", I)
print("error = ", err)
I     =  2.0
error =  2.220446049250313e-14

There are a number of other parameters to fine tune the output and operation. Do "help(quad)" for details.

Exercise

Solve for $I$: $$I = \int_0^{\pi/2}\cos{x}dx.$$

In [ ]:
 

Example

Suppose we want to integrate the function

$$f(x) = x^2 + p,$$

where $p$ is a parameter that we might want to change. That is, $f$ is really $$f(x,p) = x^2 + p,$$ but we are only integrating with respect to $x$. Like so:

$$I = \int_{x_{lo}}^{x_{hi}}f(x,p)dx$$
  • This happens all the time. For example, finding the time to fill a tank with water requires integration. But the problem depends on the tank diameter, a possible exit pipe diameter, a valve coefficient, etc.

Try solving the above problem from $x_{lo}=0$ to $x_{hi}=2$ using quad for $p=1$, $p=2$, and $p=3$.

Approach 1

  • This is possible as we have done before by leaving out the $p$ dependence and communicating $p$ as a global variable.
In [13]:
def f(x):
    return x**2 + p

p = 1; I1 = quad(f, 0, 1)[0]  # quad returns tuple (I, err), so quad()[0] is just the I part
p = 2; I2 = quad(f, 0, 1)[0]
p = 3; I3 = quad(f, 0, 1)[0]

print(I1, I2, I3)
1.3333333333333333 2.3333333333333335 3.3333333333333335

Approach 2

  • But it would be good to have our function depend on $p$ as an argument and be able to pass $p$ in properly.
  • The problem is that quad only knows to call our function as $f(x)$.
  • If we want quad to call our function with other parameters, then we need to tell quad about them.
  • Use the args argument of the quad function.
    • We set args to be a tuple of all of our extra function arguments besides $x$.
In [14]:
def f(x, p):
    return x**2 + p

I1 = quad(f, 0, 1, args=(1,))[0]     # args=(1,) sends p=1 to f(x,p) when quad calls f
I2 = quad(f, 0, 1, args=(2,))[0]     # args=(2,) sends p=2 to f(x,p) when quad calls f
I3 = quad(f, 0, 1, args=(3,))[0]     # args=(3,) sends p=3 to f(x,p) when quad calls f

print(I1, I2, I3)
1.3333333333333333 2.3333333333333335 3.3333333333333335

Now, go back and reread the bullets under Approach 2

Note

  • Many scipy tools use the args argument to pass extra function arguments along to our function.

Exercise

Solve the following: $$I = \int_0^5f(x,a,b)dx,$$ where $$f(x,a,b) = \sin(x^a)+b,$$ and $a$ and $b$ are parameters. Define function $f(x,a,b)$, the solve using quad and passing $a=2$, and $b=3$ as extra arguments.

In [ ]: