# Integração numérica

Também conhecido na literatura como quadratura.

Como calcular

$$
\int_a^b f(\xi)\:d\xi
$$

para uma função qualquer?

 * Na mão: difícil mesmo nos casos mais simples
 * No computador com cálculo simbólico: útil algumas vezes mas pouco eficiente
 * Métodos numéricos - na prática o que é usado.


<!-- TEASER_END -->

# Usar a aproximação!

É bom lembrar que chegamos em várias maneiras para aproximar uma função:

$$
y = f(x) \approx y^\delta = \sum_{i=1}^n a_i\phi_i(x)
$$

Uma primeira possibilidade é integrar termo a termo:

$$
\int_a^bf(\xi)\:d\xi \approx \sum_{i=1}^n a_i\int_a^b\phi_i(\xi)\:d\xi
$$

Como $\phi_i(x)$ foram escolhidos levando em conta simplicidade, esta abordagem pode ser utilizada. 

# Regra do ponto médio:

Dados os pontos $x_1, x_2, \ldots, x_n$, usando interpoladores constantes


```julia
using PyPlot
using CurveFit
using BenchmarkTools
```


```julia
x = 0:0.1:1
xx = 0:0.002:1
f(x) = sin(π*x) + 0.5 + cos(x)
y = f.(x);
yy = f.(xx);
```


```julia
plot(xx, yy, "b-")
plot(x, y, "rs")

nx = length(x)

lty = "k-"

for i in 1:nx-1
    xm = 0.5 * (x[i] + x[i+1])
    ym = f(xm)
    plot([x[i], x[i], x[i+1], x[i+1]], [0, ym, ym, 0], lty )
    plot([xm, xm], [0, ym], "k:")
    
end

```

Analisando a figura acima, pode-se facilmante ver que a integral vale (para pontos igualmente espaçados):

$$
\int_a^b f(\xi)\:d\xi \approx \frac{b-a}{n-1} \sum_{i=1}^{n-1}\left[ f\left(x_{i+\frac{1}{2}}\right) \right]
$$



```julia
function midrule(a, b, f, n)
    dx = (b-a) / n
    x = range(a, b, length=n+1)
    s = 0.0
    for i in 1:n
        s += f(x[i]+dx/2)
    end
    return dx * s
end
```


```julia
## Análise de convergência:
nn = 2:1000; a = 0.0; b = 1.0; 
Ie = 2/π + 0.5 + sin(1); Inn = midrule.(a, b, f, nn); err = abs.(Inn .- Ie);
loglog(nn, err, "bo"); xlabel("n"); ylabel("ε")
nn1 = 10:300; ee1 = 5 .* nn1.^(-2)
plot(nn1, ee1, "r--"); text(50, 0.005, L"$\mathcal{O}(n^{-2})$")
```

Assim, a estimativa do erro vale:

$$
\varepsilon = \mathcal{O}\left(n^{-2}\right)
$$

# Regra do trapézio

Agora, ao invés de usar interpolação consante entre os pontos, pode-se ligar os pontos por retas.


```julia
function trapezoidal(a, b, f, n)
    x = range(a, b, length=n+1)
    s = f(x[1])/2 + f(x[end])/2
    for i in 2:n
        s += f(x[i])
    end
    return (b-a) * s/n
end    
```


```julia
nn2 = nn1
Inn2 = trapezoidal.(a, b, f, nn2);
err2 = abs.(Inn2 .- Ie);
```


```julia
loglog(nn2, err2, "bo"); xlabel("n"); ylabel("ε")
mm = 10:300; ee1 = 5 .* mm.^(-2)
plot(mm, ee1, "r--")
text(50, 0.005, L"$\mathcal{O}(n^{-2})$")
```

### Estimativa de erro da regra do trapézio:

O erro ao integrar uma seção usando a regra do trapézio vale:

$$
\epsilon_i = R_i(h) - \int_{i-1}^i f(x)\:dx = f''(\xi) \frac{h^3}{12}
$$

Erro total:

$$
\epsilon = \sum_{i=1}^n \epsilon_i \approx f''(\xi)\frac{h^2(b-a)}{12} \sim h^2 
$$

# Regra de Simpson

Usa três pontos para aproximar uma parábola:

$$
\int_{i-1}^{i+1} f(x)\:dx = S_i(h) = \frac{h}{3}\left[f(x_{i-1}) + 4f(x_i) + f(x_{i+1})\right]
$$

Assim, no domínio inteiro, 

$$
\int_a^b f(x)\:dx \approx \frac{h}{3}\left[f(x_1) + 4f(x_2) + 2f(x_3) + 4f(x_4) + \ldots + 2f(x_{n-2}) + 4f(x_{n-1}) + f(x_n)\right]
$$

### Estimativa do erro:

$$
\int_{i-1}^{i+1} f(x)\:dx = \frac{h}{3}\left[f(x_{i-1}) + 4f(x_i) + f(x_{i+1})\right] - \frac{h^5}{90}f^{(4)}(\xi)
$$

$$
\epsilon_i = \mathcal{O}(h^5) \quad\Longrightarrow\quad \epsilon = \mathcal{O}(h^4)
$$



```julia
# Detalhe n tem que ser divisível por 2.
function simpsonsrule(a, b, f, n)
    x = range(a, b, length=n+1)
    s = f(x[1]) + f(x[end])
    w, w2 = 4,2
    for i = 2:n
        s += w*f(x[i])
        w, w2 = w2, w
    end
    return s * (b-a)/(3n)
end
```


```julia
nn3 = 2:2:1000
Inn3 = simpsonsrule.(a, b, f, nn3);
err3 = abs.(Inn3 .- Ie);
```


```julia
loglog(nn3, err3, "bo")
grid(ls=":")
xlabel("n")
ylabel("ε")
nn1 = 10:300
ee1 = 5 .* nn1.^(-4)
plot(nn1, ee1, "r--")
text(50, 0.000002, L"$\mathcal{O}(n^{-4})$")
```

# Generalização do processo

Usando as idéias acima com interpoladores polinomiais de ordem crescente pode-se obter regras de ordem superior para a integração. 

A idéia é usar interpoladores de Lagrange para calcular as integrais.Com isso se chega às formulas de Newton-Cotes (<https://pt.wikipedia.org/wiki/F%C3%B3rmulas_de_Newton-Cotes>)

**Cuidado com isso!!!** Lembre-se das oscilações dos polinômios de Lagrange.


# Integração de Romberg

Um processo parecido com a extrapolação de Richardson (veremos isso mais tarde). Suponhamos que queremos estimar

$$
\int_a^b f(\xi) \: d\xi
$$




Uma primeira aproximação é usar a regra do trapézio nas extremidades. O intervalo pode ser dividido por dois sequencialmente, obtendo-se:

\begin{align}
R_{1,1} &= \frac{h_1}{2}\left[f(a) + f(b)\right] \\
R_{2,1} &= \frac{h_2}{2}\left[f(a) + f(b) + 2f(a+h_2)\right]  = \frac{1}{2}\left[ R_{1,1} + h_1 f(a+h_2)\right]\\
\vdots & \\
R_{k,1} &= \frac{1}{2}\left\{ R_{k-1,1} + h_{k-1}\sum_{i=1}^{2^k-2}f\left[a + \left(2i-1\right)h_k\right]\right\}
\end{align}

onde 

$$h_k = \frac{b-a}{2^{k-1}}$$



O erro da integral vale (difícil de demonstrar mas possível):

$$
\int_a^b f(x)\:dx - R_{k,1} = \sum_{i=1}^\infty K_i h_k^{2i} = K_1 h_k^2 + \sum_{i=2}^\infty K_i h_k^{2i}
$$

por outro lado, lembrando que 

$$\epsilon_{k+1} = \frac{\epsilon_{k}}{4}$$

chega-se a:

$$
\int_a^b f(x)\:dx - R_{k+1,1} = \sum_{i=1}^\infty K_i h_{k+1}^{2i} = \sum_{i=1}^\infty \frac{K_i h_{k}^{2i}}{2^{2i}} =  \frac{K_1 h_k^2}{4} + \sum_{i=2}^\infty \frac{K_i h_k^{2i}}{4^i}
$$




substraindo a equação anterior de 4 vezes a equação acima, chega-se à seguinte relação:

$$
R_{k,2} = R_{k,1} + \frac{R_{k,1} - R_{k-1,1}}{3}
$$

Repetindo este processo:

$$
R_{k,j} = R_{k,j-1} + \frac{R_{k,j-1} - R_{k-1,j-1}}{4^{j-1}-1}
$$



```julia
function romberg1(a, b, f, n)
    R = zeros(n, n)
    h = Float64(b-a)
    R[1,1] = h/2 * (f(a) + f(b))
    
    for i in 2:n
        s = 0.0
        for k in 1:(2^(i-2))
            s += f(a + (k-0.5)*h)
        end
        R[i,1] = 0.5 * (R[i-1,1] + h * s)
        for k in 2:i
            R[i,k] = R[i,k-1] + (R[i,k-1] - R[i-1,k-1]) / (4^(k-1) - 1)
        end
        h = h/2
    end
    
    return R
end

```


```julia
romberg1(0, π, sin, 6)
```


```julia
function romberg(a, b, f, n)
    R1 = zeros(n)
    R2 = zeros(n)
    h = Float64(b-a)
    R1[1] = h/2 * (f(a) + f(b))
    for i in 2:n
        s = 0.0
        for k in 1:(2^(i-2))
            s += f(a + (k-0.5)*h)
        end
        R2[1] = 0.5 * (R1[1] + h * s)
        for k in 2:i
            R2[k] = R2[k-1] + (R2[k-1] - R1[k-1]) / (4^(k-1) - 1)
        end
        for k in 1:i
            R1[k] = R2[k]
        end
        h = h/2
    end
    
    return R2[n]
    
end

```


```julia
function countromberg(n)
    #R1 = zeros(n)
    #R2 = zeros(n)
    #h = Float64(b-a)
    #R1[1] = h/2 * (f(a) - f(b))
    countfun = 2
    countfloat = 2
    for i in 2:n
        for k in 1:(2^(i-2))
            countfun += 1
            countfloat += 4
        end
        for k in 2:i
            countfloat += 5
        end
        countfloat += 1
    end
    
    return countfun, countfloat
    
end

```


```julia
nnx = 1:20
nny = countromberg.(nnx)
nfun = [x[1] for x in nny]
nfloat = [x[2] for x in nny];

```


```julia
semilogy(nnx, nfun)
grid()
xlabel("Número de iterações")
ylabel("Número de chamadas de função")
```


```julia
nn4 = nnx
Inn4 = romberg.(a, b, f, nn4);
err4 = abs.(Inn4 .- Ie);
```


```julia
rom = loglog(nfun, err4, "bo", label="Romberg")
sim = plot(nn3, err3, "r-", label="Regra de Simpson")
mid = plot(nn, err, "g-", label="Regra do ponto médio")
grid()
xlabel("Número de Chamadas de função")
ylabel("ε")
legend()
```


```julia
## Benchmark para comparar tempo de execução:
I1a = midrule(a, b, f, 400)
I2a = simpsonsrule(a, b, f, 22)
I3a = romberg(a, b, f, 4)

abs(I1a - Ie), abs(I2a - Ie), abs(I3a - Ie)
```


```julia
@benchmark midrule(a, b, f, 400)
```


```julia
@benchmark simpsonsrule(a, b, f, 22)
```


```julia
@benchmark romberg(a, b, f, 4)
```


```julia
I1b = midrule(a, b, f, 800000)
I2b = simpsonsrule(a, b, f, 1000)
I3b = romberg(a, b, f, 6)

abs(I1b - Ie), abs(I2b - Ie), abs(I3b - Ie)
```


```julia
@benchmark midrule(a, b, f, 800000)
```


```julia
@benchmark simpsonsrule(a, b, f, 1000)
```


```julia
@benchmark romberg(a, b, f, 6)
```

# Problemas

## Problema 1

Faça uma função para calcular a integral usando a regra do ponto médio com pontos cujo espaçamento não é igual



```julia

```

## Problema 2
Faça uma função para calcular a integral usando a regra do trapézio com pontos cujo espaçamento não é igual



```julia

```

## Problema 3

Calcule usando os diversos métodos apresentados a seguinte integral:

$$
\int_{-1}^1 \left(1- | x |\right) \: dx
$$





```julia

```

## Problema 4

A regra do ponto médio apresenta o mesmo erro assintótico que a regra do trapézio. Vamos explorar um pouco isso:

 1. Você consegue criar um método com erro $\mathcal{O}(h)$?
 2. Se você criou o método acima do jeito mais simples você consegue criar um outro método bem semelhante com o mesmo tipo de erro.
 3. Qual a relação entre estes dois métodos e a regra ponto médio
 



```julia

```

e97f98e322af90cbe2b788baf28950ee87d50606
ipynb
Jupyter Notebook
03-integracao-numerica.ipynb
pjabardo/sci-comp
2f590363d0b7edd87d8125494eebac1c346da78d
MIT
03-integracao-numerica.ipynb
pjabardo/sci-comp
2f590363d0b7edd87d8125494eebac1c346da78d
MIT
03-integracao-numerica.ipynb
pjabardo/sci-comp
2f590363d0b7edd87d8125494eebac1c346da78d
MIT
Qwen/Qwen-72B
1. YES
2. YES