Code Session 8: Integração

import numpy as np 
import matplotlib.pyplot as plt

A integração numérica em uma variável pode ser realizada em Python utilizando a função quad do módulo scipy.integrate.

`quad`¶

Esta função calcula a integral definida $\int_a^b f(x) \, dx$ numericamente através de regras de quadratura.

Os argumentos de entrada obrigatórios desta função são:

a função f a ser integrada
o limite inferior a
o limite superior b

Os principais argumentos de saída são:

y: valor numérico da integral
abserr: estimativa do erro absoluto

Como importá-la?

from scipy.integrate import quad

from scipy.integrate import quad

Problema 1¶

O período de um pêndulo simples de comprimento $L$ é $\tau = 4\sqrt{ \dfrac{L}{g} } h(\theta_0)$ , onde $g$ é a aceleração gravitacional, $\theta_0$ representa a amplitude angular, e

h(\theta_0) = \int_0^{\pi/2} \dfrac{1}{ \sqrt{ 1 - \textrm{sen}^2(\theta_0/2)\textrm{sen}^2 \theta } } \, d\theta

(1)

Calcule $h(\tau)$ , para $\tau = 15^{\circ}, \, 30^{\circ}, \, 45^{\circ}$ e compare estes valores com $h(0^{\circ}) = \pi/2$ (a aproximação usada para pequenas amplitudes).

Resolução¶

Em primeiro lugar, fazemos os cálculos diretos da integral para os distintos valores de $\tau$ .

# cálculo direto das integrais caso a caso

theta0 = np.array([0,15,30,45]) # ângulos

vals,errs = [],[] # integrais, erros

for t in theta0:
    f = lambda theta: 1/(np.sqrt(1. - np.sin(t/2)**2 * np.sin(theta)**2))
    v,e = quad(f,0,np.pi/2)
    print(f'Integral h({t}) = {v:g}')
    vals.append(v)
    errs.append(e)

# converte listas para arrays
vals = np.asarray(vals)
errs = np.asarray(errs)

Integral h(0) = 1.5708
Integral h(15) = 2.49203
Integral h(30) = 1.79372
Integral h(45) = 1.67896

Vemos que o valor das integrais é muito sensível. Para realizar uma comparação mais interessante, utilizaremos um cálculo relativo tomando o valor em $h(15)$ como referência.

plt.stem(theta0, (vals - vals[0])/vals[0], use_line_collection=True);
plt.xlabel('$\\theta_0$')
plt.ylabel('erro relativo');

/var/folders/ll/g0vl8b194pbfyp4cwpzz4p740000gn/T/ipykernel_51604/945931442.py:1: MatplotlibDeprecationWarning: The 'use_line_collection' parameter of stem() was deprecated in Matplotlib 3.6 and will be removed two minor releases later. If any parameter follows 'use_line_collection', they should be passed as keyword, not positionally.
  plt.stem(theta0, (vals - vals[0])/vals[0], use_line_collection=True);

Podemos verificar que a mudança de valor da integral entre $15^{\circ}$ e $30^{\circ}$ é da ordem de 60%, enquanto que nos demais casos, ela se limita a 20%.

Problema 2¶

Uma corrente elétrica alternada é descrita por

i(t) = i_0 \left( \textrm{sen}\left( \frac{\pi t}{t_0} \right) - \beta \, \textrm{sen}\left( \frac{2 \pi t}{t_0} \right) \right),

(2)

onde $i_0 = 1 \, A$ , $t_0 = 0.05 \, s$ e $\beta =0.2$ . Calcule a corrente RMS definida por

i_{rms} = \sqrt{ \dfrac{1}{t_0} \int_0^{t_0} i^2(t) \, dt}

(3)

Resolução¶

Neste caso, basta passarmos os valores iniciais e finais para computar a integral.

i0, t0, beta = 1.0, 0.05, 0.2 # parâmetros iniciais

i2 = lambda t: ( i0 * ( np.sin( (np.pi*t) / t0 ) - beta * np.sin( (2*np.pi*t) / t0 ) ) )**2 # função

i_rms = np.sqrt( 1.0/t0 * quad(i2, 0, t0)[0] )

print(f'Corrente RMS = {i_rms:g} A')

Corrente RMS = 0.72111 A

Regra do Trapézio Generalizada¶

A regra do trapézio generalizada (composta) pode ser calculada usando

scipy.integrate.cumtrapz

Vamos utilizar a função $i(t)$ do Problema 2 e estimar sua integral no intervalo $t = [0,5]$ utilizando a regra do trapézio generalizada.

# Visualização
t = np.linspace(0,5)
i = lambda t: i0 * ( np.sin( (np.pi*t) / t0 ) - beta * np.sin( (2*np.pi*t) / t0 ) )  # função

plt.plot(t,i(t),'o');
plt.xlabel('$t$')
plt.ylabel('$i(t)$');

from scipy.integrate import cumtrapz

i = lambda t: i0 * ( np.sin( (np.pi*t) / t0 ) - beta * np.sin( (2*np.pi*t) / t0 ) ) # função

T = cumtrapz(i(t), t)[-1] # pega último valor, já que é cumulativa

print(f'Integral por Trapézio = {T:g}')

Integral por Trapézio = -4.97605e-15

Quadratura Gaussiana (QG)¶

O cálculo de uma integral por quadratura Gaussiana pode ser calculado usando

scipy.integrate.quadrature

Vamos utilizar a função $i(t)$ do Problema 2 e estimar sua integral no intervalo $t = [0,5]$ utilizando a QG para várias ordens (controladas pelo argumento miniter).

from scipy.integrate import quadrature

for ordem in range(1,11):
    I_QG,err_QG = quadrature(i,0,5,miniter=ordem)
    print(f'Integral por QG (ordem {ordem}) = {I_QG:g}')

Integral por QG (ordem 1) = 3.60822e-14
Integral por QG (ordem 2) = 9.71445e-15
Integral por QG (ordem 3) = 4.35416e-14
Integral por QG (ordem 4) = 6.92502e-14
Integral por QG (ordem 5) = 6.57807e-14
Integral por QG (ordem 6) = -5.7801e-14
Integral por QG (ordem 7) = -1.40166e-14
Integral por QG (ordem 8) = -1.60635e-14
Integral por QG (ordem 9) = -7.63278e-16
Integral por QG (ordem 10) = -4.0341e-14

quad¶

Problema 1¶

Resolução¶

Problema 2¶

Resolução¶

Regra do Trapézio Generalizada¶

Quadratura Gaussiana (QG)¶

`quad`¶