T6: Integration

T6: Integration#

Note

Click the and open this notebook in Colab to enable interactivity.

Note

To save your progress, make a copy of this notebook in Colab File > Save a copy in Drive and you’ll find it in My Drive > Colab Notebooks.

Integration is one of the fundamental concepts of calculus, so it’s important that you can do it analytically and numerically. We will once again be relying on the powerful SciPy library, specifically the scipy.integrate module.

Example 1#

Find the centroid of the region bounded by the curve \( y = \sin x\) and the parabola \(y = x^2\). You are given that \(\sin x = x^2\) has \(0\) and \(x_0 \approx 0.8767\) as the only two solutions.

In the workbook, you used the trapezoid(y, x) function for integrating fixed samples, i.e., y was an array of discrete values. Here, we will demonstrate a slightly more general approach to numerical integration, using the quad(func, a, b) function (short for quadrature) to integrate func between a and b. funcis once again a callable that we will have to define ourselves.

While that is doable for an integral in one dimension, we also need to compute an integral in two dimensions to get the moments and center of mass. For that, we will use the general dblquad(func, a, b, gfun, hfun) function, which integrates:

  • the function func(y, x), where y must be first and x is second;

  • from x between a and b;

  • from y between gfun and hfun, both callable functions that we have to write.

import numpy as np
from scipy.integrate import quad, dblquad

x0 = 0.8767

def ylower(x):
    return x ** 2

def yupper(x):
    return np.sin(x)
    
def func(x):
    return yupper(x) - ylower(x)

A = quad(func, 0, x0)[0]
print(f"A = {A:.4f}\n")


def func_xA(y, x):
    return x

def func_yA(y, x):
    return y

xA = dblquad(func_xA, 0, x0, ylower, yupper)[0]
print(f"xA = {xA:.4f}")
print(f"xcom = {xA/A:.4f}\n")

yA = dblquad(func_yA, 0, x0, ylower, yupper)[0]
print(f"yA = {yA:.4f}")
print(f"ycom = {yA/A:.4f}")

A = 0.1357

xA = 0.0601
xcom = 0.4431

yA = 0.0445
ycom = 0.3277

Example 2#

Evaluate \(\displaystyle\int_{0}^{\infty} e^{-x^2} dx\).

def f(x):
    return np.exp(-x**2)

res = quad(f, 0, np.inf)[0]
print(f"Integral: {res:.4f}")
print(f"Numerical √π/2: {np.sqrt(np.pi)/2:.4f}")

Integral: 0.8862
Numerical √π/2: 0.8862

Example 3#

Find the moment of inertia about the \(z\)-axis of the region bounded by parabolas \(4 + y = \dfrac{1}{2}(1 + x)^2\) and \(4 - y = \dfrac{1}{2} (1 - x)^2\) where the density function is \(\delta(x,y) = \dfrac{xy}{x^2 + y^2}\).

The first thing we need to do is to find the domain of integration, which means the points where the two parabolas intersect need to be determined. Straightforward equation solving gives \((\sqrt{7}, \sqrt{7})\) and \((-\sqrt{7}, -\sqrt{7})\) as the points of intersection.

It is also much easier to integrate \(y\) first, then \(x\), as it gives us:

\[ \int_{-\sqrt{7}}^{\sqrt{7}} \int_{\frac{1}{2}(1+x)^2-4}^{4 - \frac{1}{2}(1-x)^2} f(x,y)\, \mathrm{d}y\, \mathrm{d}x, \]

where \(f(x,y) = xy\).

import numpy as np
from scipy.integrate import dblquad

def ylower(x):
    return 1/2 * (1 + x)**2 - 4

def yupper(x):
    return 4 - 1/2 * (1 - x)**2

def fxy(y, x):
    return x * y

I = dblquad(fxy, -np.sqrt(7), np.sqrt(7), ylower, yupper)[0]
print(f"Integral: {I:.4f}")
print(f"Numerical: {196 / 15 * np.sqrt(7):.4f}")

Integral: 34.5712
Numerical: 34.5712
Contents