# Solução dos exercícios propostos na aula 4

Ao fim da aula 4, foram propostos alguns exercícios envolvendo interpolação. A seguir estão os exercícios. Em seguida, vem uma breve recapitulação da aula. A solução está logo após.

## Exercícios

### Problema 1

Interpole a função de Runge com $-1 \le x \le 1$:
$$
f(x) = \frac{1}{1 + 25x^2}
$$

 1. Use 11 pontos uniformemente distribuídos
 2. Aumente o número de pontos
 3. Tente usar os pontos $x_k = \cos\left(\frac{k\pi}{N}\right)$ para $k = 0\ldots N$.
 4. Brinque com o número de pontos

### Problema 2

Procure na Net o método de diferenças divididas de Newton a interpole a função anterior nos mesmos pontos. Este método é simplesmente um jeito inteligente de resolver a matriz apresentada lá em cima.

### Problema 3

Use a biblioteca Interpolations.jl e Dierckx.jl para fazer as interpolações. Compare a interpolação linear com os splines.



<!-- TEASER_END -->

# Interpolação

Dado um conjunto de n pontos $(x_i, y_i)$, qual o poliômio que passa por todos?

$$
y_i = a_0 + a_1 x_i + a_2 x_i^2 + \ldots a_n x_i^n \qquad i=1, \ldots, m
$$

## Vandermonde

Com n+1 pontos distintos, se o polinômio for de grau n, pode-se montar o seguinte sistema linear:

$$
\begin{bmatrix}
1 & x_0 & x_0^2 & \cdots & x_0^n \\
1 & x_1 & x_1^2 & \cdots & x_1^n \\
\vdots & \vdots & \vdots & \ddots & \vdots\\
1 & x_n & x_n^2 & \cdots & x_n^n\\
\end{bmatrix}\cdot
\left\{ \begin{matrix} a_0 \\ a_1 \\ \vdots \\ a_n\\ \end{matrix}\right\}
= \left\{\begin{matrix} y_0 \\ y_1 \\ \vdots \\ y_n\\\end{matrix}\right\}
$$

Mas esta é uma operação cara, $\mathcal{O}(n^3)$!

## Outra possibilidade

$$
y = f(x) \approx a_0 + a_1(x-x_0) + a_2(x-x_0)(x-x_1) + \cdots + a_n(x-x_0)(x-x_1)\cdots(x - x_{n-1})
$$

Com isso chegamos ao seguinte sistema linear triangular:
$$
\begin{bmatrix}
1 & 0 & 0 &  0 &\cdots & 0\\
1 & (x_1-x_0) & 0 & 0 &  \cdots & 0\\
1 & (x_2 - x_0) &  (x_2 - x_0)(x_2 - x_1) & 0 & \cdots & 0\\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots\\
1 & (x_n - x_0) &  (x_n - x_0)(x_n - x_1) & (x_n - x_0)(x_n - x_1)(x_n-x_2) & \cdots &(x_n-x_0)(x_n-x_1)\cdots(x_n-x_{n-1})\\
\end{bmatrix}\cdot\left\{\begin{matrix} a_0\\ a_1 \\ a_2 \\ \vdots \\ a_n\end{matrix}\right\} = 
\left\{\begin{matrix} y_0\\ y_1 \\ y_2 \\ \vdots \\ y_n\end{matrix}\right\}
$$

Resolver este sistema é muito mais barato: $\mathcal{O}(n^2)$

## Interpolação de Lagrange:

$$
y(x) = \sum_{i=1}^n y_i h_i(x)
$$

onde $h_i(x)$ é o interpolador de Lagrange:

$$
h_k(x) = \prod_{i=1\ldots n,}^n \frac{x - x_i}{x_k - x_i} \qquad i\ne k
$$

Propriedade:
$$
h_i(x_i) = \delta_{ij} \quad \text{onde} \quad \delta_{ij} = \left\{\begin{matrix}1, \: i=j \\ 0, i\ne j\\ \end{matrix}\right.
$$

### Vamos organizar a interpolação de Lagrange


```julia
struct Lagrange
    x::Vector{Float64}
    y::Vector{Float64}
    Lagrange(x, y) = new(copy(x), copy(y))
end
Base.Broadcast.broadcastable(lgr::Lagrange) = Ref(lgr)

function lagrange(k, z, x)
 h = 1.0
    n = length(z)
    for i = 1:(k-1)
        h *= (x - z[i]) / (z[k] - z[i])
    end
    
    for i = (k+1):n
        h *= (x - z[i]) / (z[k] - z[i])
    end
    return h
end

function interp(lgr::Lagrange, x)
    
    y = lgr.y[1] * lagrange(1, lgr.x, x)
    
    for i = 2:length(lgr.x)
        y += lgr.y[i] * lagrange(i, lgr.x, x)
    end
    
    return y
end

    
(lgr::Lagrange)(x) = interp(lgr, x)    
```

## Interpolação Linear




```julia
struct LinearInterp
    x::Vector{Float64}
    y::Vector{Float64}
    LinearInterp(x, y) = new(copy(x), copy(y))
end
Base.Broadcast.broadcastable(lin::LinearInterp) = Ref(lin)

function interp1(lin::LinearInterp, x)
    
    if x < lin.x[1] || x > lin.x[end]
        error("Fora do Range")
    end
    
    index = 2
    n = length(lin.x)
    for i = 2:n
        if lin.x[i] >= x
            index = i
            break
        end
    end
    i1 = index-1
    return lin.y[i1] + (lin.y[index] - lin.y[i1]) * (x - lin.x[i1]) / (lin.x[index] - lin.x[i1])
    
end
    
function interp2(lin::LinearInterp, x)
    index = searchsortedfirst(lin.x, x)
    if index == 1
        index = 2
    end
    i1 = index-1
    return lin.y[i1] + (lin.y[index] - lin.y[i1]) * (x - lin.x[i1]) / (lin.x[index] - lin.x[i1])
    
end

(lin::LinearInterp)(x) = interp1(lin, x)

```

# Pacotes

 * [Interpolations](https://github.com/JuliaMath/Interpolations.jl)
 * [Dierckx](https://github.com/kbarbary/Dierckx.jl)
 * [GridInterpolations](https://github.com/sisl/GridInterpolations.jl)


```julia
using PyPlot
using BenchmarkTools
using LinearAlgebra
```

# Problema 1


```julia
f(x) = 1.0 / (1.0 + 25x^2)

xx = -1:0.001:1
yy = f.(xx);
```

## Pontos igualmente espaçados


```julia
function ploteqspacing(N, dx=0.001)
    xp = range(-1.0, 1.0, length=N)
    yp = f.(xp)
    lgr = Lagrange(xp, yp)
    x1 = -1:dx:1
    y1 = lgr.(x1)
    y0 = f.(x1)
    plot(xp, yp, "ro")
    plot(x1, y0, "r-")
    plot(x1, y1)    
end

```


```julia
ploteqspacing(11)

```


```julia
ploteqspacing(5)
```


```julia
ploteqspacing(8)
```


```julia
ploteqspacing(15)
```


```julia
ploteqspacing(16)
```

## Distribuição do cosseno



```julia
function plotcosspacing(N, dx=0.001)
    xp = [cos(k*π/N) for k in 0:N]
    yp = f.(xp)
    lgr = Lagrange(xp, yp)
    x1 = -1:dx:1
    y1 = lgr.(x1)
    y0 = f.(x1)
    plot(xp, yp, "ro")
    plot(x1, y0, "r-")
    plot(x1, y1)    
end
```


```julia
plotcosspacing(10)
```


```julia
plotcosspacing(4)
```


```julia
plotcosspacing(8)
```


```julia
plotcosspacing(16)
```


```julia
plotcosspacing(20)
```

# Problema 2

Diferenças divididas


```julia
function divdiffmat(x)
    n = length(x)
    
    A = zeros(n,n)
    for i = 1:n
        A[i,1] = 1.0
        for k in 2:i
            A[i, k] = A[i,k-1] * (x[i] - x[k-1])
        end
    end
    
    return A
    
end

```


```julia
struct DividedDiff
    x::Vector{Float64}
    a::Vector{Float64}
end
function divideddiff(x, y)
    A = LowerTriangular(divdiffmat(x))
    a = A\y
    return DividedDiff(copy(x), a)
end
Base.Broadcast.broadcastable(ddif::DividedDiff) = Ref(ddif)


function interp(ddif, x)
    xx = ddif.x
    a = ddif.a
    
    y = a[end] 
    
    for i in (lastindex(a)-1):-1:1
        y = a[i] + (x-xx[i])*y
    end
    return y
end
(ddif::DividedDiff)(x) = interp(ddif, x)
```


```julia
fun2(x) = 1.0 + 2.0*x + 3.0*x*(x-1.0) + 4.0*x*(x-1.0)*(x-2.0) + 5.0*x*(x-1.0)*(x-2.0)*(x-3.0)
x1 = [0.0, 1.0, 2.0, 3.0, 4.0]
y1 = fun2.(x1);
```


```julia
diff = divideddiff(x1, y1)
```


```julia
x2 = 0:0.01:4
y2 = fun2.(x2)
y2b = diff.(x2);
```


```julia
plot(x1, y1, "rs")
plot(x2, y2, "r-")
plot(x2, y2b, "b:")

```


```julia
x = -1:0.4:1
y = sin.(π*x);
```


```julia
diff = divideddiff(x, y)
```


```julia
y1 = diff.(x);
```


```julia
xx = -1:0.01:1
yy = sin.(π*xx)
yy1 = diff.(xx);
```


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

## Algorítmo de Newton

Um jeito interessante de calcular a interpolação usando diferenças divididas. Ao se resolver o sistema linear, pode-se escrever a solução como:

$$
\begin{align}
a_0 &= f(x_0)\\
a_1 & = \frac{f(x_1) - f(y_0)}{x_1 - x_0}\\
a_2 &= \frac{\frac{f(x_2)-f(x_0)}{x_2-x_0} - \frac{f(x_1) - f(x_0)}{x_1-x_0} }{x_2 - x_1}\\
\vdots &=  \vdots\\
\end{align}
$$

o k-ésimo coeficiente vale:
$$
a_k = \mathcal{F}\left(x_0, x_1, \ldots, x_k\right)
$$

com 
$$
\mathcal{F}\left(x_0, x_1, \ldots, x_k\right) = \frac{\mathcal{F}\left(x_0, x_1, \ldots, x_{k-1}\right) -\mathcal{F}\left(x_1, x_1, \ldots, x_k\right)}{x_0-x_k}
$$


```julia
function newton_divdiff(x, y)
    
    n = length(x)
    F = zeros(n,n)
    
    for i in 1:n
        F[i,1] = y[i]
    end
    
    for i = 2:n
        for k in 2:i
            F[i,k] = (F[i,k-1] - F[i-1,k-1]) / (x[i] - x[i-k+1])
        end
    end
    a = [F[i,i] for i in 1:n]
    
    return DividedDiff(copy(x), a)
    
end

```


```julia
fun2(x) = 1.0 + 2.0*x + 3.0*x*(x-1.0) + 4.0*x*(x-1.0)*(x-2.0) + 5.0*x*(x-1.0)*(x-2.0)*(x-3.0)
x1 = [0.0, 1.0, 2.0, 3.0, 4.0]
y1 = fun2.(x1);
diff2 = newton_divdiff(x1, y1)

```

## Uma abordagem que aloca menos memória

Mesmo algorítmo anterior mas não aloca a matriz inteira.


```julia
function newton_divdiff2(x, y)
    
    n = length(x)
    F = zeros(n)
    F0 = zeros(n)
    F1 = zeros(n)
    
    F0[1] = y[1]
    F1[1] = y[1]
    F[1] = y[1]
    
    for i = 2:n
        F1[1] = y[i]
        for k in 2:i
            F1[k] = (F1[k-1] - F0[k-1]) / (x[i] - x[i-k+1])
        end
        F[i] = F1[i]
        for k in 1:i
            F0[k] = F1[k]
        end
    end
    
    return DividedDiff(copy(x), F)
    
end

```


```julia
fun2(x) = 1.0 + 2.0*x + 3.0*x*(x-1.0) + 4.0*x*(x-1.0)*(x-2.0) + 5.0*x*(x-1.0)*(x-2.0)*(x-3.0)
x1 = [0.0, 1.0, 2.0, 3.0, 4.0]
y1 = fun2.(x1);
diff2 = newton_divdiff2(x1, y1)

```

# Problema 3

Usando a biblioteca Interpolations


```julia
using Interpolations
```


```julia
x = -1:0.2:1
y = f.(x);
```

# Interpolação com malha (*grid*)


```julia
itp1 = interpolate((x,),  y, Gridded(Constant()));
itp2 = interpolate((x,),  y, Gridded(Linear()));
```


```julia
xx = -1:0.001:1
yy = f.(xx)
yy1 = itp1.(xx);
yy2 = itp2.(xx);
plot(x, y, "rs")
plot(xx, yy, "r-")
plot(xx, yy1, "b-")
plot(xx, yy2, "b--")

```


```julia
itp3 = CubicSplineInterpolation((x,), y; bc=Line(OnGrid()));

yy3 = itp3.(xx)

plot(x, y, "rs")
plot(xx, yy, "r-")
plot(xx, yy3, "b:")
```


```julia
function compare(h, x)
    xp = -1:h:1
    yp = f.(xp)
    
    
    itp1 =  interpolate((xp,),  yp, Gridded(Constant()));
    itp2 =  interpolate((xp,),  yp, Gridded(Linear()));
    itp3 = CubicSplineInterpolation((xp,), yp, bc=Flat(OnGrid()));
    fx = f(x)
    return itp1(x)-fx, itp2(x)-fx, itp3(x)-fx
end
```


```julia
hstep = [1.0, 0.5, 0.4, 0.2, 0.1, 0.05, 0.02, 0.01, 0.001]
nn = 2.0 ./ hstep
erra = [compare(h, -1+h/4) for h in hstep]
errb = [compare(h, h/4) for h in hstep];
errc = [compare(h, -0.9905) for h in hstep];
errd = [compare(h, 0.0005) for h in hstep];
```


```julia
ea1 = [abs(e[1]) for e in erra]
ea2 = [abs(e[2]) for e in erra]
ea3 = [abs(e[3]) for e in erra]

eb1 = [abs(e[1]) for e in errb]
eb2 = [abs(e[2]) for e in errb]
eb3 = [abs(e[3]) for e in errb];

ec1 = [abs(e[1]) for e in errc]
ec2 = [abs(e[2]) for e in errc]
ec3 = [abs(e[3]) for e in errc];

ed1 = [abs(e[1]) for e in errd]
ed2 = [abs(e[2]) for e in errd]
ed3 = [abs(e[3]) for e in errd];

```


```julia
loglog(nn, ea1, "ro-")
loglog(nn, ea2, "bo-")
loglog(nn, ea3, "go-")


```


```julia
loglog(nn, eb1, "rs-")
loglog(nn, eb2, "bo-")
loglog(nn, eb3, "g^-")

```


```julia
loglog(nn, ec1, "rs-")
loglog(nn, ec2, "bo-")
loglog(nn, ec3, "g^-")

```


```julia
loglog(nn, ed1, "rs-")
loglog(nn, ed2, "bo-")
loglog(nn, ed3, "g^-")

```

93c2667e01f66fe7d7b525b8c6fbf4cd008bdd74
ipynb
Jupyter Notebook
04-solucao.ipynb
iptscicomp/aulas
41a5ac19dcf4e16118f0e5ecf9b389353e8de347
MIT
04-solucao.ipynb
iptscicomp/aulas
41a5ac19dcf4e16118f0e5ecf9b389353e8de347
MIT
04-solucao.ipynb
iptscicomp/aulas
41a5ac19dcf4e16118f0e5ecf9b389353e8de347
MIT
2019-03-15T04:10:15.000Z
2019-03-15T04:10:15.000Z
Qwen/Qwen-72B
1. YES
2. YES