---
title: "Análise Numérica - Trabalho Prático 3"
author:
- Diogo Cordeiro
- Hugo Sales
- Pedro Costa
geometry: margin=2cm
output: pdf_document
---

### Motivação

Pretende-se compreender o funcionamento conceptual bem como os desafios da implementação de dois métodos de integração
numéricos: `Regra de Simpson` e `Regra dos Trapézios.`

### 1

<!-- $$ \frac{d^4}{dx^4}_{x \in [0;2]}sin(sin(sin(sin(x)))) $$ -->

Através da análise do gráfico abaixo, verificamos que o majorante, em valor absoluto, da 4ª derivada
da função é menor que 12, este valor é usado para majorar a formula do erro para o calculo de n.

![Quarta derivada da função enunciada](majorante_erro_1.png){ width=8cm }


    /* Função que implementa o método de Simpson */
    double simpson (Function f, double a, double b, int n)
    {
        // Intervalo de passo
        double h = (b - a)/n;
        // Valor de f nos pontos de indice par
        double evens = summation(f, 2, n - 2, 2, h, a);
        // Valor de f nos pontos de indice impar
        double odds = summation(f, 1, n - 1, 2, h, a);
        // Aplicação da método de Simpson
        value = (h/3)*(f(a) + f(b) + 2 * evens + 4 * odds);
    }

\pagebreak

    /* Função para calcular o sumatório de F entre os pontos de indice init
       e stop saltando step pontos, usando h como intervalo de passo */
    double summation (Function f, int init, int stop, int step, double h, double a)
    {
        // Acumulador
        double total = 0;
        for (int i = init; i <= stop; i += step) {
            // Adicionamos o valor de f correspondente ao x de indice i
            total += f(a + i * h);
        }
        return total;
    }

    /* Função para calcular o número de pontos necessários para o calculo
       do integral com erro menor que error */
    int calculateN (double A, double B, double error) {
        int n = ceil( (B - A) / pow((15.0 * error) / 2.0, 1.0 / 4) );
        return n + (n % 2);
    }

    void main ()
    {
        Function f(x) = sin(sin(sin(sin(x))));
        // com 7 casas decimais correctas
        print(simpson(f, 0, 2, calculateN(0, 2, pow(10, -7))));
        // com 12 casas decimais correctas
        print(simpson(f, 0, 2, calculateN(0, 2, pow(10, -12))));
    }

#### Output

Erro           | Resultado
---------------|-------------------
$$ 10^{-7} $$  | 1.0548418906594816
$$ 10^{-12} $$ | 1.0548418772492483

### 2)

    // Valor exacto do integral calculado com o WolframAlpha, arredondado
    // com 15 algarismos significativos, um a mais do que o erro majorado
    // máximo para o caso de 2^20 pontos.
    double I = 1.05484187724912;

    /* Função para calcular o integral recorrendo ao método do Trapézio */
    double trapezio (Function f, double a, double b) {
        // Acumulador
        double summation = 0;

        // Intervalo de passo com n inicial = 2
        double h = (b - a)/2;

        // Valor constante
        double fa_fb = f(a) + f(b);

        // Para cada expoente de 1 a 20 (com passo 1)
        for (int k = 1; k <= 20; ++k) {
            // Guarda sumatorio dos pontos anteriores
            double partial_sum = 0;
            // Número de intervalos
            // Com o left shift fazemos a potência de 2^k
            int n = 1 << k;
            // Evitamos recalcular pontos da função previamente computados
            // guardando o sumatório destes em partial_sum e adicionamos a cada
            // iteração os pontos novos, sendo estes de indice impar
            for (int i = 1; i < n; i += 2) {
                // Adicionamos o valor de f correspondente ao x de indice i
                partial_sum += f(a + i * h);
            }

            summation += partial_sum;

            // Aplicação da formula do Trapézio
            double value = (h/2)*fa_fb + h * summation;

            print(k + "\t| " + value + "\t| " + (I - value));

            // Dividimos o intervalo por 2
            h /= 2;
        }
    }

    void main () {
        Function f(x) = sin(sin(sin(sin(x))));
        print(trapezio(f, 0, 2);
    }

#### Output

k   | $$ I_{n_{k}} $$     | $$ \mid I - I_{n_{k}} \mid $$
----|---------------------|--------------------------
1   | 0.9533749638740736  | $$ 1.1 \cdot 10^{-1}  $$
2   | 1.0308378382617962  | $$ 2.5 \cdot 10^{-2}  $$
3   | 1.0489039934457873  | $$ 6.0 \cdot 10^{-3}  $$
4   | 1.053360809734676   | $$ 1.5 \cdot 10^{-3}  $$
5   | 1.0544718169560368  | $$ 3.8 \cdot 10^{-4}  $$
6   | 1.054749374997165   | $$ 9.3 \cdot 10^{-5}  $$
7   | 1.0548187524860693  | $$ 2.4 \cdot 10^{-5}  $$
8   | 1.0548360961083287  | $$ 5.8 \cdot 10^{-6}  $$
9   | 1.054840431967042   | $$ 1.5 \cdot 10^{-6}  $$
10  | 1.0548415159287925  | $$ 3.6 \cdot 10^{-7}  $$
11  | 1.0548417869190467  | $$ 9.1 \cdot 10^{-8}  $$
12  | 1.0548418546666     | $$ 2.3 \cdot 10^{-8}  $$
13  | 1.054841871603487   | $$ 5.7 \cdot 10^{-9}  $$
14  | 1.0548418758377078  | $$ 1.5 \cdot 10^{-9}  $$
15  | 1.0548418768962615  | $$ 3.5 \cdot 10^{-10} $$
16  | 1.0548418771608998  | $$ 8.9 \cdot 10^{-11} $$
17  | 1.0548418772270582  | $$ 2.3 \cdot 10^{-11} $$
18  | 1.0548418772436017  | $$ 5.6 \cdot 10^{-12} $$
19  | 1.0548418772477444  | $$ 1.4 \cdot 10^{-12} $$
20  | 1.0548418772487873  | $$ 3.3 \cdot 10^{-13} $$