Aproximações de Números e Funções por Sequências e Séries

Os objetivos desta lição são

  • Identificar como lidar com sequências e séries em programação;

  • Entender que elementos anteriores da sequência são substituídos;

  • Tentar gastar o mínimo de elementos e variáveis (mantendo a clareza);

  • Tomar cuidado com overflow e underflow;

  • Relembrar/treinar a programação.

Sequências

Uma sequência em matemática é um conjunto de valores ordenados:

\[ (a_1, a_2, a_3, \dots, a_n, \dots) \]

O índice \(1\) usado no primeiro elemento é arbitrário, e serve para normalizar essa ordem. Em alguns casos é mais interessante começar com o índice \(0\).

Matematicamente, costumamos estar interessados na convergência dessa sequência. No entanto, podemos utilizar uma sequência por outros motivos. Computacionalmente, estamos interessados em como obter os elementos de uma sequência, quando necessário.

É importante ressaltar aqui a importância da frase quando necessário, pois como veremos a seguir, não é do nosso interesse calcular todos esses valores ao mesmo tempo, e quase sempre também não seria possível.

Exemplo 1

Dado um objeto com aceleração constante de \(10 m/s^2\) saindo do repouso, queremos calcular sua posição nos instantes \(t = 0, 0.1, 0.2, 0.3, \dots\). Note que temos uma sequência \(t_n = (n-1)/10\), e queremos uma sequência \(x_n\) das posições.

Neste caso, porém, as leis de movimento nos dizem que

\[ x(t) = x_0 + v_0 t + \frac{1}{2}at^2, \]

que em nosso caso se reduz a

\[ x(t) = 5t^2. \]

Em outras palavras, a sequência \(x_n\) pode ser calculada diretamente dado \(t_n\), e \(t_n\) também tem uma fórmula geral dado \(n\),

\[ x_n = 5\bigg(\frac{n-1}{10}\bigg)^2 = \frac{(n-1)^2}{20}. \]

Desse modo, não é necessário guardar nada em especial, pois \(x_n\) pode ser calculado diretamente a partir de \(n\).

Exemplo 2

Dado um número \(x \in \mathbb{N}, x > 1\), queremos calcular a sequência de Collatz desse número. Em particular, queremos calcular qual o primeiro elemento que é \(1\) (se houver).

Aqui é importante ressaltar a conjectura de Collatz: Dado um número \(a_1 \in \mathbb{N}\), e a fórmula recursiva

\[ a_{n+1} = \left\{\begin{array}{ll} a_n/2, & \text{se } a_n \text{ é par}; \\ 3a_n + 1, & \text{se } a_n \text{ é impar}, \end{array}\right. \]

sempre existe um elemento desta sequência tal que \(a_n = 1\). A partir desse ponto a sequência perde a importância. Por exemplo, a sequência de Collatz de 3 é

\[ 3 \rightarrow 10 \rightarrow 5 \rightarrow 16 \rightarrow 8 \rightarrow 4 \rightarrow 2 \rightarrow 1 \rightarrow 4 \rightarrow 2 \rightarrow 1 \cdots \]

Nessa sequência, é necessário saber quem é o ponto atual, e é impossível saber um elemento arbitrário da sequência em casos gerais. Dessa maneira é preciso calcular elemento a elemento.

Lembre-se que para verificar a paridade de um número, utilizamos %.

15 % 2
1
14 % 2
0
15 % 2 == 0 # É par ?
false
14 % 2 == 0 # É par ?
true

Lembre-se também que para dividir um número e mantê-lo inteiro, usamos div.

div(15, 2)
7
div(14, 2)
7
14/2 # 7.0 e 7 são o mesmo valor, mas são diferentes em computação
7.0

Um erro comum de alunos inexperientes nesse caso é tentar numerar todos os elementos.

a1 = 3
if a1 % 2 == 0
    a2 = div(a1, 2)
else
    a2 = 3 * a1 + 1
end
10
if a2 % 2 == 0
    a3 = div(a2, 2)
else
    a3 = 3 * a2 + 1
end
5

Essa estratégia logo fica inviável. No entanto, note que precisamos apenas do valor atual da sequência, de maneira que podemos usar apenas uma mesma variável a que irá guardar sempre o valor mais atual.

a = 300

if a % 2 == 0
    a = div(a, 2)
else
    a = 3 * a + 1
end
a

if a % 2 == 0
    a = div(a, 2)
else
    a = 3 * a + 1
end
a

if a % 2 == 0
    a = div(a, 2)
else
    a = 3 * a + 1
end
a
226

Rodando mais de uma vez esse bloco, fazemos uma nova iteração da sequência. Para deixar mais claro, vamos guardar também uma variável n que diz qual o índice da sequência.

a = 3
n = 1

if a % 2 == 0
    a = div(a, 2)
else
    a = 3 * a + 1
end
n = n + 1
println("a = $a, n = $n")

if a % 2 == 0
    a = div(a, 2)
else
    a = 3 * a + 1
end
n = n + 1
println("a = $a, n = $n")

if a % 2 == 0
    a = div(a, 2)
else
    a = 3 * a + 1
end
n = n + 1
println("a = $a, n = $n")
a = 10, n = 2
a = 5, n = 3
a = 16, n = 4

Note que n também pode ser usado para responder a pergunta original: quando a for 1, teremos o n marcando o índice correspondente.

Isso nos leva a um conceito importante. Note que sequências contém infinitos elementos. No entanto, computacionalmente, não podemos fazer contas infinitas. Isso constitui um loop infinito. Códigos com loops infinitos são considerados "errados" (no sentido de não fazerem o que se espera), então devemos evitá-los. Portanto, o fato de percebermos que para resolver o problema, devemos fazer a conta do Bloco 2 até que n seja 1 é de extrema importância, pois isso evitará um loop infinito (dado que a conjectura de Collatz esteja correta).

Com isso, podemos definir um termo formal para este conjunto de código com um objetivo específico:

Def. (Algoritmo): Um algoritmo é uma coleção de instruções para realizar alguma tarefa específica. Segundo Knuth (The Art of Computer Programming, v.1), ele deve satisfazer as seguintes condições, parafraseadas aqui:

  • Finitude: O algoritmo deve acabar em tempo finito;

  • Bem definido: As intruções devem ser claras e sem ambiguidade;

  • Entrada: O algoritmo tem zero ou mais entradas, que são valores determinados antes do algoritmo começar. Essas entradas são especificados a partir de conjuntos de objetos;

  • Saída: O algoritmo tem uma ou mais saídas, que são quantidades relacionadas com as entradas;

  • Eficácia/Computabilidade: As operações feitas no algoritmo devem ser suficientemente básicas para que a princípio possam ser executadas por uma pessoa num espaço finito e tempo finito com papel e caneta.

Podemos fazer uma prévia do nosso algoritmo:

  1. Entrada: \(x\)

  2. Inicialização: \(a \leftarrow x, n \leftarrow 1\).

  3. Enquanto \(a > 1\)

    1. Se \(a\) é par, faça \(a \leftarrow a/2\)

    2. Caso contrário, faça \(a \leftarrow 3a + 1\).

    3. Incremente \(n\).

  4. Saída: \(n\).

Nosso algoritmo para encontrar o primeiro elemento está quase feito. A entrada, saída, e todos os passos estão claros. Todos os passos são triviamente feitos com papel e caneta. No entanto, temos alguns pequenos problemas:

  • E se a conjectura for falsa?

Se a conjectura é verdadeira, o algoritmo para quando encontra 1, que pela conjectura, existe. Mas se não existe, pode ser que tenhamos um loop infinito. Dessa maneira, iremos colocar uma condição de parada para a falha do algoritmo. Não existe uma regra para a condição de parada. Deve-se tomar cuidado, no entanto, para não excluir a possibilidade do algoritmo funcionar em muitos casos. Uma condição de parada possível neste caso é que n seja muito grande. Se, por exemplo, fizermos 1 milhão de iterações e não encontrarmos a_n = 1, então talvez ele não exista. Note que essa condição exclui algumas soluções de aparecerem, por exemplo

\[ a_1 = 2^{1.000.000}. \]

É fácil ver que \(a_n = 2^{1.000.001 - n}\), de modo que para n = 1.000.001 teremos \(a_n = 1\).

Muitas vezes a condições de parada de falha é um limitante físico para deixar o algoritmo tratável. É muito fácil escolher um número x tal que nosso algoritmo leve horas para convergir. Não é do nosso interesse esperar tanto para um algoritmo recreativo. Poderíamos colocar uma condição de tempo então. Qual condição escolher é um assunto complicado, que depende de muitos outros fatores além da disciplina. Em especial, a teoria matemática arredor do problema deve ser considerada. Imagine que a conjectura de Collatz é falsa, e que x é um valor finito tal que a sequência gerada com a_1 = x nunca decresça à 1. Como é possível verificar que x é um contra-exemplo para a conjectura?

  • O maior inteiro.

O computador é uma máquina física, com limitações impostas para que seja possível fazer contas determinísticas. Dessa maneira, ele segue regras para definir números inteiros e "reais" (as aspas serão explicadas depois), e essas regras limitam o maior inteiro que pode ser representado por um tipo específico de dado. Mais importante, esse valor é independente da linguagem. A saber, o tipo de inteiro básico do Julia é o Int64, de 64 bits, e cujo maior valor é \(2^{63}-1 \approx 10^{18}\). Coisas estranhas acontecem se fizermos qualquer conta com esse valor, de modo que é imperativo não ultrapassá-lo. Vamos colocar uma nova condição de parada no nosso código: se \(a > 10^{17}\), também paramos.

  1. Entrada: \(x \in \mathbb{N}\)

  2. Inicialização: \(a \leftarrow x, n \leftarrow 1, a\_\text{bound} = 10^17, n\_\text{bound} = 10^6\)

  3. Se \(a \leq 0\) faça \(n \leftarrow 0\) termine o algoritmo

  4. Enquanto \(a > 1\)

    1. Se \(a\) é par, faça \(a \leftarrow a/2\)

    2. Caso contrário, faça \(a \leftarrow 3a + 1\).

    3. Incremente \(n\).

    4. Se \(a > a\_\text{bound}\) faça \(n \leftarrow -1\) e termine o algoritmo

    5. Se \(n > n\_\text{bound}\) faça \(n \leftarrow -2\) e termine o algoritmo

  5. Saída: \(n\) Se \(n > 0\), então \(a_n = 1\): Saída bem sucedida, Se \(n = 0\), então \(x \leq 0\): Saída de erro na entrada, Se \(n = -1\), então \(a\) ficou muito grande: Saída de valor muito grande, Se \(n = -2\), então \(n\) ficou muito grande: Saída de muitas iterações.

Keyword argument

function collatz(a::Int; a_bound = 10^17, n_bound = 1_000_000) # O mesmo que começar com x e atribuir x à a.
    n = 1
    if a <= 0
        return 0
    end
    while a > 1
        print("$a ")
        if a % 2 == 0
            a = div(a, 2)
        else
            a = 3a + 1
        end
        n += 1
        if a > a_bound
            return -1
        elseif n > n_bound
            return -2
        end
    end
    return n
end
collatz (generic function with 1 method)
collatz(151)
151 454 227 682 341 1024 512 256 128 64 32 16 8 4 2 
16
collatz(2^30, n_bound=10)
1073741824 536870912 268435456 134217728 67108864 33554432 16777216 8388608 4194304 2097152 
-2
collatz(1)
1
collatz(10^17+2)
100000000000000002 50000000000000001 
-1

Fibonacci

Vamos continuar nosso exemplo com mais uma sequência interessante, a de Fibonacci:

\[ F_1 = F_2 = 1 \qquad F_{n+1} = F_n + F_{n-1}. \]

A sequência de Fibonacci envolve dois termos, e cresce infinitamente. Existem vários motivos para se trabalhar com a série de Fibonacci, mas vamos utilizá-la recreativamente. Vamos criar um simples algoritmo para calcular o n-ésimo termo da série de Fibonacci. Assuma, por enquanto, que todos os termos serão bem representados por inteiros de 64 bits.

Entrada: \(n\) inteiro

  1. Se \(n \leq 0\) retorne \(0\)

  2. Se \(n = 1\) ou \(n = 2\) retorne \(1\)

  3. Inicialização: \(F_1 \leftarrow 1, F_2 \leftarrow 1, k \leftarrow 2\)

  4. Enquanto \(k < n\)

    1. \(F_{\text{novo}} \leftarrow F_1 + F_2\) - Isso calcula o valor novo

    2. \(F_1 \leftarrow F_2\) - O valor \(F_2\) passa a ser o mais antigo dos dois

    3. \(F_2 \leftarrow F_{\text{novo}}\) - O valor \(F_{\text{novo}}\) passa a ser o anterior

    4. Incremente \(k\)

  5. Retorne \(F_2\). Saída: 0 se a entrada estava incorreta, \(F_n\) caso contrário.

function fibon(n::Int)
    if n ≤ 0 #\le
        return 0
    elseif n == 1 || n == 2
        return 1
    end
    F1, F2, k = 1, 1, 2
    while k < n
        F3 = F1 + F2
        F1, F2 = F2, F3
        #F1, F2 = F2, F1 + F2
        k += 1
    end
    return F2
end

for n = 1:20
    println("F_$n = $(fibon(n))")
end
F_1 = 1
F_2 = 1
F_3 = 2
F_4 = 3
F_5 = 5
F_6 = 8
F_7 = 13
F_8 = 21
F_9 = 34
F_10 = 55
F_11 = 89
F_12 = 144
F_13 = 233
F_14 = 377
F_15 = 610
F_16 = 987
F_17 = 1597
F_18 = 2584
F_19 = 4181
F_20 = 6765

É bastante simple calcular um elemento da sequência de Fibonacci, e é interessante que podemos explorar outra maneira de fazer isso.

Recursão

function fat_por_recursao(n::Int)
    if n < 0
        error("BLAH")
    end
    if n == 0
        return 1
    end
    return n * fat_por_recursao(n - 1)
end

fat_por_recursao(5)
120
fat_por_recursao(21)
-4249290049419214848

operador ternário

PROP ? Resultado se SIM : Resultado se NAO

5 > 0 ? "ok" : "oops"
"ok"
fat_uma_linha(n::Int) = n == 0 ? 1 : n * fat_uma_linha(n - 1)
fat_uma_linha (generic function with 1 method)
fat_uma_linha(5)
120

Recursão é uma das partes mais importantes no desenvolvimento de códigos. Muitas estratégias computacionais e matemáticas envolvem o chamado "dividir e conquistar" ou "reduzir para um caso conhecido". Em particular, demonstrações por indução finita fazem uma coisa parecida, onde você supõe saber resolver para um caso e consegue resolver um caso de nível superior reduzindo-o ao caso conhecido.

Computacionalmente, podemos pensar em recursão como a divisão de um problema em outros, cada um destes menor que o primeiro. Para cada problema menor, repetimos o processo. Como não podemos ter um loop infinito, é necessário terminar essa redução de alguma maneira. Essa maneira envolve, simplesmente, saber resolver alguns casos de problema menor (que garantam a convergência).

No caso de Fibonacci, o cálculo do n-ésimo termo é feito calculando os termos n-1 e n-2. Os casos n-1 e n-2 são, de fato, menores. Além disso, repetindo este processo chegaremos a dois casos: n = 1 ou n = 2. Desse modo, temos

  • casos pequenos que conseguimos resolver trivialmente;

  • uma maneira de reduzir um problema em outros menores;

  • garantia matemática que essa redução leva aos problemas menores que sabemos resolver.

Abaixo temos um pseudo-código.

Rotina \(\text{FIBO}(n)\)

Entrada: \(n\)

  1. Se \(n = 1\) ou \(n = 2\) retorne \(1\)

  2. Senão, retorna \(\text{FIBO}(n-1) + \text{FIBO}(n-2)\)

function FIBO(n)
    if n <= 2
        return min(1, max(n, 0)) # Pequeno truque
    else
        return FIBO(n-1) + FIBO(n-2)
    end
end

FIBO(10)

for n = 1:20
    println("F_$n = $(FIBO(n))")
end
F_1 = 1
F_2 = 1
F_3 = 2
F_4 = 3
F_5 = 5
F_6 = 8
F_7 = 13
F_8 = 21
F_9 = 34
F_10 = 55
F_11 = 89
F_12 = 144
F_13 = 233
F_14 = 377
F_15 = 610
F_16 = 987
F_17 = 1597
F_18 = 2584
F_19 = 4181
F_20 = 6765

Muitas vezes, o código de recursão será mais curto. Em particular, em Julia (e C) temos o chamado operador ternário, que nos permite fazer a pergunta if ... else ... end em uma única linha. Desse modo, temos uma implementação curtíssima de Fibonacci:

F(n) = n <= 2 ? min(1, max(n, 0)) : F(n-1) + F(n-2)

for n = 1:20
    println("F_$n = $(F(n))")
end
F_1 = 1
F_2 = 1
F_3 = 2
F_4 = 3
F_5 = 5
F_6 = 8
F_7 = 13
F_8 = 21
F_9 = 34
F_10 = 55
F_11 = 89
F_12 = 144
F_13 = 233
F_14 = 377
F_15 = 610
F_16 = 987
F_17 = 1597
F_18 = 2584
F_19 = 4181
F_20 = 6765

No entanto, é importante ter cuidado ao utilizar recursão. O custo computacional fica escondido nas chamadas de função, e às vezes podemos deixar um código muito mais complicado do que deveria. No caso de Fibonacci, por exemplo, ao calcular F(n), pedimos o cálculo de F(n-1) e F(n-2), mas o cálculo de F(n-1) irá pedir o cálculo de F(n-2) e F(n-3), ou seja, F(n-2) será calculado duas vezes!

function FIBO_comprint(n)
    println("Chamada: F$n")
    if n <= 2
        return min(1, max(n, 0)) # Pequeno truque
    else
        return FIBO_comprint(n-1) + FIBO_comprint(n-2)
    end
end

FIBO_comprint(5)
Chamada: F5
Chamada: F4
Chamada: F3
Chamada: F2
Chamada: F1
Chamada: F2
Chamada: F3
Chamada: F2
Chamada: F1
5

Em outras palavras, Fibonacci com recursão acaba sendo muito mais caro que Fibonacci usando for ou while. Tome cuidado com essas armadilhas.

Leitura adicional: Existe um conceito chamada "Avaliação preguiçosa", que consiste na avaliação dos valores somente quando são necessários. Em Julia, o pacote Lazy.jl implementa esse conceito. Se a avaliação de F(n) for feita de maneira preguiçosa, não haverá o custo adicional do cálculo de \(F_n\) repetidos.

Mais Fibonacci

Vamos voltar nossa atenção ao problema inicial de Fibonacci: calcular o n-ésimo termo. Nossas últimas discussões foram todas no sentido de como calcular os elementos de Fibonacci termo a termo ou recursivamente. Porém, matematicamente, o n-ésimo termo de Fibonacci é bem definido!

\[ F_n = \frac{\phi^n - \psi^n}{\sqrt{5}}, \]

onde

\[\phi = \frac{1 + \sqrt{5}}{2}\]

e

\[\psi = \frac{1 - \sqrt{5}}{2} = 1 - \phi. \]

Isso quer dizer que podemos calcular Fibonacci com uma fórmula direta.

phi = (1 + sqrt(5))/2
psi = 1 - phi
for n = 1:20
    Fn = round(Int, (phi^n - psi^n)/sqrt(5))
    println("F_$n = $Fn")
end
F_1 = 1
F_2 = 1
F_3 = 2
F_4 = 3
F_5 = 5
F_6 = 8
F_7 = 13
F_8 = 21
F_9 = 34
F_10 = 55
F_11 = 89
F_12 = 144
F_13 = 233
F_14 = 377
F_15 = 610
F_16 = 987
F_17 = 1597
F_18 = 2584
F_19 = 4181
F_20 = 6765

Perceba que o problema aqui agora é outro. Os valores não estão mais "corretos". Na verdade, quase todos contém um pequeno erro.

Quão pequeno?

phi = (1 + sqrt(5))/2
psi = 1 - phi
F1 = F2 = 1
Fnovo = 1
for n = 3:93
    Fn = (phi^n - psi^n)/sqrt(5)
    Fnovo = F1 + F2
    F2 = F1
    F1 = Fnovo
    println("F_$n = $Fnovo, Diferença F_$n = $(Fn - Fnovo)")
end
F_3 = 2, Diferença F_3 = 0.0
F_4 = 3, Diferença F_4 = 0.0
F_5 = 5, Diferença F_5 = 8.881784197001252e-16
F_6 = 8, Diferença F_6 = 1.7763568394002505e-15
F_7 = 13, Diferença F_7 = 1.7763568394002505e-15
F_8 = 21, Diferença F_8 = 3.552713678800501e-15
F_9 = 34, Diferença F_9 = 7.105427357601002e-15
F_10 = 55, Diferença F_10 = 1.4210854715202004e-14
F_11 = 89, Diferença F_11 = 2.842170943040401e-14
F_12 = 144, Diferença F_12 = 5.684341886080802e-14
F_13 = 233, Diferença F_13 = 5.684341886080802e-14
F_14 = 377, Diferença F_14 = 1.7053025658242404e-13
F_15 = 610, Diferença F_15 = 3.410605131648481e-13
F_16 = 987, Diferença F_16 = 4.547473508864641e-13
F_17 = 1597, Diferença F_17 = 9.094947017729282e-13
F_18 = 2584, Diferença F_18 = 1.8189894035458565e-12
F_19 = 4181, Diferença F_19 = 2.7284841053187847e-12
F_20 = 6765, Diferença F_20 = 4.547473508864641e-12
F_21 = 10946, Diferença F_21 = 7.275957614183426e-12
F_22 = 17711, Diferença F_22 = 1.0913936421275139e-11
F_23 = 28657, Diferença F_23 = 2.1827872842550278e-11
F_24 = 46368, Diferença F_24 = 3.637978807091713e-11
F_25 = 75025, Diferença F_25 = 5.820766091346741e-11
F_26 = 121393, Diferença F_26 = 8.731149137020111e-11
F_27 = 196418, Diferença F_27 = 1.7462298274040222e-10
F_28 = 317811, Diferença F_28 = 2.9103830456733704e-10
F_29 = 514229, Diferença F_29 = 4.656612873077393e-10
F_30 = 832040, Diferença F_30 = 8.149072527885437e-10
F_31 = 1346269, Diferença F_31 = 1.1641532182693481e-9
F_32 = 2178309, Diferença F_32 = 2.3283064365386963e-9
F_33 = 3524578, Diferença F_33 = 3.725290298461914e-9
F_34 = 5702887, Diferença F_34 = 6.51925802230835e-9
F_35 = 9227465, Diferença F_35 = 1.1175870895385742e-8
F_36 = 14930352, Diferença F_36 = 1.862645149230957e-8
F_37 = 24157817, Diferença F_37 = 2.9802322387695312e-8
F_38 = 39088169, Diferença F_38 = 4.470348358154297e-8
F_39 = 63245986, Diferença F_39 = 6.705522537231445e-8
F_40 = 102334155, Diferença F_40 = 1.341104507446289e-7
F_41 = 165580141, Diferença F_41 = 2.086162567138672e-7
F_42 = 267914296, Diferença F_42 = 3.5762786865234375e-7
F_43 = 433494437, Diferença F_43 = 5.960464477539062e-7
F_44 = 701408733, Diferença F_44 = 1.0728836059570312e-6
F_45 = 1134903170, Diferença F_45 = 1.6689300537109375e-6
F_46 = 1836311903, Diferença F_46 = 2.6226043701171875e-6
F_47 = 2971215073, Diferença F_47 = 4.76837158203125e-6
F_48 = 4807526976, Diferença F_48 = 7.62939453125e-6
F_49 = 7778742049, Diferença F_49 = 1.33514404296875e-5
F_50 = 12586269025, Diferença F_50 = 1.9073486328125e-5
F_51 = 20365011074, Diferença F_51 = 3.4332275390625e-5
F_52 = 32951280099, Diferença F_52 = 5.340576171875e-5
F_53 = 53316291173, Diferença F_53 = 9.1552734375e-5
F_54 = 86267571272, Diferença F_54 = 0.000152587890625
F_55 = 139583862445, Diferença F_55 = 0.000274658203125
F_56 = 225851433717, Diferença F_56 = 0.00042724609375
F_57 = 365435296162, Diferença F_57 = 0.00067138671875
F_58 = 591286729879, Diferença F_58 = 0.0010986328125
F_59 = 956722026041, Diferença F_59 = 0.0018310546875
F_60 = 1548008755920, Diferença F_60 = 0.0029296875
F_61 = 2504730781961, Diferença F_61 = 0.0048828125
F_62 = 4052739537881, Diferença F_62 = 0.00830078125
F_63 = 6557470319842, Diferença F_63 = 0.013671875
F_64 = 10610209857723, Diferença F_64 = 0.021484375
F_65 = 17167680177565, Diferença F_65 = 0.037109375
F_66 = 27777890035288, Diferença F_66 = 0.0625
F_67 = 44945570212853, Diferença F_67 = 0.09375
F_68 = 72723460248141, Diferença F_68 = 0.171875
F_69 = 117669030460994, Diferença F_69 = 0.28125
F_70 = 190392490709135, Diferença F_70 = 0.4375
F_71 = 308061521170129, Diferença F_71 = 0.6875
F_72 = 498454011879264, Diferença F_72 = 1.1875
F_73 = 806515533049393, Diferença F_73 = 2.0
F_74 = 1304969544928657, Diferença F_74 = 3.0
F_75 = 2111485077978050, Diferença F_75 = 5.25
F_76 = 3416454622906707, Diferença F_76 = 8.5
F_77 = 5527939700884757, Diferença F_77 = 14.0
F_78 = 8944394323791464, Diferença F_78 = 24.0
F_79 = 14472334024676221, Diferença F_79 = 40.0
F_80 = 23416728348467685, Diferença F_80 = 64.0
F_81 = 37889062373143906, Diferença F_81 = 104.0
F_82 = 61305790721611591, Diferença F_82 = 160.0
F_83 = 99194853094755497, Diferença F_83 = 272.0
F_84 = 160500643816367088, Diferença F_84 = 448.0
F_85 = 259695496911122585, Diferença F_85 = 736.0
F_86 = 420196140727489673, Diferença F_86 = 1216.0
F_87 = 679891637638612258, Diferença F_87 = 2048.0
F_88 = 1100087778366101931, Diferença F_88 = 3200.0
F_89 = 1779979416004714189, Diferença F_89 = 5120.0
F_90 = 2880067194370816120, Diferença F_90 = 8192.0
F_91 = 4660046610375530309, Diferença F_91 = 14336.0
F_92 = 7540113804746346429, Diferença F_92 = 22528.0
F_93 = -6246583658587674878, Diferença F_93 = 1.844674407370959e19

Impressionantemente, esse erro não é tão pequeno assim. Ele cresce com o tamanho de n, de modo que essa aproximação é inviável para valores grandes de n.

A fórmula para o n-ésimo número de Fibonacci é

\[ F_n = \frac{\phi^n - \psi^n}{\sqrt{5}}. \]

Acontece que \(\phi > 1\) e \(0 < \psi < 1\), então a distância entre \(\phi^n\) e \(\psi^n\) cresce bastante com \(n\). Além disso, ambos são irracionais, então seus valores já são aproximados. Neste caso, é mais simples calcular os números de Fibonacci utilizando for ou while com inteiros. Claro que isso também vai depender do motivo para se calcular esses números.

Aproximações de Séries

Vamos para um assunto levemente diferente. Uma série é uma soma de infinitos termos numa ordem dada.

\[ S = \sum_{k = 1}^{\infty} a_k. \]

A questão aqui é, em geral, se essa soma existe ou não. Para tanto, define-se uma sequência \((s_1,s_2,\dots,s_n,\dots)\) dada por

\[ s_n = \sum_{k = 1}^n a_k. \]

Se essa sequência convergence, então essa série converge.

Computacionalmente, em geral, estamos mais preocupados em calcular essa soma. Um resultado básico de série diz que para que a série convirja, é necessário que \(a_k \rightarrow 0\). Como já vimos, isso irá nos causar alguns problemas, dado que

\[ s_{n+1} = s_n + a_{n+1}, \]

e \(a_{n+1}\) vai eventualmente ser muito pequeno.

Matematicamente, muitas vezes temos um objetivo, como calcular \(\pi\), que pode ser obtido pelo cálculo de uma série. Devemos nos preocupar em como fazer isso de uma maneira computacional eficiente que não perca muita precisão. Frequentemente, teremos limites computacionais para essa eficiência, e aí devemos voltar à matemática para conseguir alguma maneira melhor de resolver o problema inicial.

Expansão de Taylor

Como deve ter sido visto no curso de Cálculo, uma função continuamente diferenciável até ordem \(n\) admite uma aproximação polinomial em torno de um ponto \(a\) de seu domínio dada por

\[ P_n(x) = f(a) + f'(a)(x-a) + \frac{1}{2}f''(a)(x-a)^2 + \frac{1}{3!}f'''(a)(x-a)^3 + \dots + \frac{1}{n!}f^{(n)}(a)(x-a)^n. \]

Esse polinômio é chamado polinômio de Taylor de ordem \(n\) em torno do ponto \(a\), e existem alguns teoremas indicando o quão boa é essa aproximação.

Teorema: Se \(f\) é continuamente diferenciável até ordem \(n\) no ponto \(a\), então

\[ f(x) = P_n(x) + r_n(x), \]

onde

\[ \lim_{x \rightarrow a} \frac{ r_n(x) }{ |x - a|^n } = 0. \]

Teorema: Se \(f\) é continuamente diferenciável até ordem \(n+1\) no intervalo fechado de \(a\) à \(x\), então

\[ f(x) = P_n(x) + \int_a^x \frac{f^{(n+1)}(t)}{n!}(x-t)^n \text{d} t. \]

Teorema: Se \(f\) é continuamente diferenciável até ordem \(n+1\) num intervalo aberto contento \(a\) e \(f^{(n)}\) é contínua no intervalo fechado de \(a\) à \(x\), então

\[ f(x) = P_n(x) + \frac{1}{(n+1)!} f^{(n+1)}(\xi)(x - a)^{n + 1}, \]

onde \(\xi\) é um número real entre \(a\) e \(x\).

Podemos usar o conceito da expansão de Taylor para calcular numericamente alguns valores de funções não polinomiais. O caso mais comum é o da função \(\exp(x) = e^x\), onde \(e\) é o número de Euler.

MathConstants.e
ℯ = 2.7182818284590...

Aviso: A função exponencial já costuma estar implementada em baixo nível, no entanto faremos este estudo por questões didáticas.

A função \(e^x\) tem a seguinte expansão de Taylor em torno do ponto \(0\).

\[ e^x = 1 + x + \frac{x^2}{2} + \frac{x^3}{3!} + \frac{x^4}{4!} + \dots + \frac{x^n}{n!} + \dots. \]

Então é bastante simples definir um algoritmo que calcule a aproximação da exponecial dado \(x\) e \(n\).

Entrada: x e n ≧ 0

  1. Inicialização: \(E \leftarrow 1.0\)

  2. Para k de 1 à n

    1. \(E \leftarrow E + x^k/k!\)

  3. Retorne E Saída: E ≈ eˣ com n termos da expansão de Taylor

E = 1.0
E += 1.0   # k = 1
E += 0.5   # k = 2
E += 1/6   # k = 3
E += 1/24  # k = 4
E += 1/120 # k = 5
2.7166666666666663
MathConstants.e - E
0.0016151617923787498
function exponencial(x, n)
    #Exercício. Usa a função factorial(n) para calcular n!
    E = 1.0
    for k = 1:n
        E += x^k / factorial(k)
    end
    return E
end

exponencial(1.0, 20)
2.7182818284590455
MathConstants.e - exponencial(1.0, 20)
-4.440892098500626e-16
abs( exp(2) - exponencial(2.0, 20) ) / exp(2)
6.250497655727703e-15
abs( exp(3) - exponencial(3.0, 20) ) / exp(3)
1.1790767393200277e-11
abs( exp(10) - exponencial(10.0, 20) ) / exp(10)
0.0015882606618580852
abs( exp(3) - exponencial(3.0, 21) ) / exp(3)
OverflowError: 21 is too large to look up in the table; consider using `factorial(big(21))` instead

Vamos testar nossa função contra exp, calculando o erro dessa aproximação. No entanto, note que a função exponencial cresce rapidamente, e como vimos, a precisão de um valor é relativo ao valor. Sendo assim, vamos utilizar o erro relativo na nossa comparação.

Definição: O erro cometido ao se aproximar \(x\) por \(x'\) é \(x - x'\).

Definição: O erro absoluto cometido ao se aproximar \(x\) por \(x'\) é \(|x - x'|\).

Definição: O erro relativo cometido ao se aproximar \(x\) por \(x'\) é \(\dfrac{|x - x'|}{|x|}\).

using Plots
gr(size=(400,300))

N = collect(1:20)
Erro = [abs(exponencial(1.0, n) - exp(1.0))/exp(1.0) for n in N]
scatter(N, Erro, yaxis=:log)

Note que após \(n = 17\), o erro é sempre a precisão da máquina. Isso quer dizer que os termos \(a_n\) com \(n > 17\) não afetam a soma, computacionalmente.

N = collect(1:20)
Erro = [abs(exponencial(5.0, n) - exp(5.0))/exp(5.0) for n in N]
scatter(N, Erro, yaxis=:log)
png(joinpath(@OUTPUT, "plt-exp2"))

Veja que é erro é muito maior para \(x = 5\) do que para \(x = 1\). Isso acontece porque os termos além de \(n > 20\) ainda são importantes para a aproximação de \(e^5\). A primeira vista, uma simples solução seria utilizar \(n > 20\), mas veja o que acontece com nossa implementação.

Se você se lembra bem, isso acontece porque para \(21!\) é maior que o máximo dos inteiros de 64 bits. Isso quer dizer que para calcular \(e^x\) corretamente, devemos dar mais atenção aos detalhes do algoritmo. Note também que \(x^n\) pode acabar ficando muito grande para \(x\) e \(n\) grandes e causar overflow também.

Revisitando nossa soma:

\[ s_n = 1 + x + \frac{x^2}{2} + \dots + \frac{x^n}{n!}, \]

de modo que

\[ s_n = s_{n-1} + \frac{x^n}{n!}. \]

O termo à direita da soma não pode ser calculado diretamente pois cada termo da fração pode "explodir" (termo que utilizarei para dizer causar overflow). No entanto, o resultado da fração é bem comportado, então podemos tentar chegar nesse valor de outra maneira.

Veja que

\[ \frac{x^n}{n!} = \frac{x\times x\times x\times\dots\times x}{1\times2\times3\times\dots\times n} = \frac{x}{1}\times\frac{x}{2}\times\frac{x}{3}\times\dots\times\frac{x}{n}, \]

então podemos fazer esse cálculo seguindo essa ordem, de maneira que o produto todo ficará equilibrado.

No entanto, é possível deixar esse produto mais eficiente. Note que se chamarmos \(t_n = \dfrac{x^n}{n!}\), então

\[ s_n = s_{n-1} + t_n, \]

e

\[ t_n = t_{n-1}\times\frac{x}{n}. \]

Em outras palavras, a sequência que é somada também pode ser calculada utilizando o termo anterior. Cada iteração fará então apenas um produto, uma divisão, e uma soma.

Gráfico do termo

x = 10.0
n = 60
E = [1.0]
t = 1.0
for k = 1:n
    t = t * x / k
    push!(E, t)
end
scatter(E, yaxis=:log)

Entrada: \(x\) e \(n \geq 0\)

  1. Inicialização: \(E \leftarrow 1.0, t \leftarrow 1.0\)

  2. Para \(k\) de \(1\) à \(n\)

    1. \(t \leftarrow t * x / k\)

    2. \(E \leftarrow E + t\)

  3. Retorne \(E\) Saída: \(E \approx e^x\) com n termos da expansão de Taylor

function exponencial2(x, n)
    #Implemente
    E = 1.0
    t = 1.0
    for k = 1:n
        t *= x / k
        E += t
    end
    return E
end

exponencial2(1.0, 20)
2.7182818284590455
(exponencial2(10.0, 60) - exp(10.0)) / exp(10)
-1.6516398231937218e-16
using Plots
N = collect(1:50)
Erro = [abs(exponencial2(1.0, n) - exp(1.0))/exp(1.0) for n in N]
Erro[findall(Erro .== 0.0)] .= eps()
scatter(N, abs.(Erro), yaxis=:log)
N = collect(1:50)
Erro = [abs(exponencial2(10.0, n) - exp(10.0))/exp(10.0) for n in N]
Erro[findall(Erro .== 0.0)] .= eps()
scatter(N, Erro, yaxis=:log)
ylims!(eps()/2, 1.0)

Perceba agora que a função exp não recebe um \(n\) específico. A quantidade de termos calculados é específico para o valor \(x\). Para \(x = 0\), \(n = 0\) basta, para \(x = 1\), \(n = 17\) basta, para \(x = 10\), \(n = 45\) foi necessário.

Em particular, como estamos fazendo uma atualização do tipo \(s_n = s_{n-1} + t_n\) e \(t_n\) tende a 0, podemos verificar se obtivemos um erro de arredondamento na soma, de modo que os termos seguintes também não acrescentarão nada.

Entrada: \(x\)

  1. Inicialização: \(E \leftarrow 1.0, t \leftarrow x, k \leftarrow 1\)

  2. Enquanto \(E + t \neq E\)

    1. \(E \leftarrow E + t\)

    2. Incremente \(k\)

    3. \(t \leftarrow t * x / k\)

  3. Retorne \(E\) Saída: \(E \approx e^x\)

function exponencial3(x)
    #Implemente
    E = 1.0
    t = x
    k = 1
    while E + t != E
        E += t
        k += 1
        t *= x / k
    end
    return E
end

plot(x->exp(x) - exponencial3(x), 0, 10.0)
plot(x->abs(exp(x) - exponencial3(x))/exp(x), 0, 10.0)

Por fim, vejamos o que acontece com \(x < 0\).

plot(x->abs(exp(x) - exponencial3(x))/exp(x), -10.0, 10.0)

O erro cresce novamente. Isso acontece porque o termo \(t_n\) fica alternando de sinal, e fica pequeno antes de contribuir o suficiente para fazer a diferença necessária. Uma maneira de remediar esse problema é usar a relação \(e^x = \dfrac{1}{e^{-x}}\).

function exponencial4(x)
    if x < 0
        return 1.0/exponencial4(-x)
    end
    #Implemente
    return exponencial3(x)
end

plot(x->abs(exp(x) - exponencial4(x))/exp(x), -10.0, 10.0)

Cálculo de \(\pi\)

Um tópico bastante interessante de computação matemática é o cálculo do valor de \(\pi\). O interesse no cálculo de \(\pi\) é antigo, e existem várias maneiras de fazê-lo. Além disso, a busca pelo valor de \(\pi\) com a maior quantidade de dígitos corretos é uma disputa matemática bastante acirrada.

Das maneiras de calcular \(\pi\), uma das mais interessantes, em minha opinião, é através da identidade

\[ \arctan 1 = \frac{\pi}{4} \qquad \Rightarrow \qquad \pi = 4 \arctan 1. \]

Agora, utilizamos a derivada de \(\arctan\):

\[ \frac{\text{d}}{\text{d}x}\arctan x = \frac{1}{1+x^2}, \]

e a expansão da fração para \(|x| < 1\):

\[ \frac{1}{1+x^2} = 1 - x^2 + x^4 - x^6 + \dots + (-1)^n x^{2n} + \dots. \]

Integrando, temos

\[ \arctan x = x - \frac{x^3}{3} + \frac{x^5}{5} - \frac{x^7}{7} + \dots + (-1)^n\frac{x^{2n+1}}{2n+1} + \dots\]

É possível mostrar que essa série converge para \(x = 1\) também, de modo que

\[ \frac{\pi}{4} = \arctan 1 = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \dots \frac{(-1)^n}{2n+1} + \dots. \]

Utilizando conceitos parecidos com o da expansão da exponencial, podemos calcular \(\pi\) sem precisar passar um valor de \(n\), porém, essa convergência é muito lenta. No exemplo abaixo, colocaremos uma condição de parada para \(n > 1.000.000.000\).

function pi_atan()
    S = 1.0
    n = 1
    σ = -1
    t = σ/(2n+1)
    while S + t != S
        S += t
        σ = -σ
        n += 1
        t = σ/(2n+1)
        if n > 1_000_000_000
            break
        end
    end
    return 4S
end

pi_atan()
3.1415926545880506
@time pi_atan() - pi
  1.182689 seconds
9.982574766809194e-10

O valor é, de fato, calculado até uma aproximação razoável, mas demora muito e leva muitas iterações. Por sorte, existem dezenas de outras maneiras de se calcular \(\pi\).

Uma dessas maneiras é a série

\[ \pi = \sqrt{12}\bigg(1 - \frac{1}{3\times 3} + \frac{1}{5\times 3^2} - \frac{1}{7\times 3^3} + \frac{1}{9\times 3^4} - \dots\bigg) = \sqrt{12}\sum_{k=0}^{\infty}\frac{(-1)^k}{(2k+1)3^k}. \] 
function pi_madhava()
    S = 1.0
    third = 1.0/3.0
    σ = -1
    pow_third = third
    t = σ * third * third
    n = 1
    while S + t != S
        S += t
        n += 1
        pow_third *= third
        σ = -σ
        t = σ * pow_third / (2n + 1)
    end
    println("n = $n") # Usar para mostrar quantas iterações
    return sqrt(12) * S
end

pi_madhava() - pi
n = 31
8.881784197001252e-16

Note que com 31 iterações já chegamos num valor de \(\pi\) decente.

\(\pi\) com BigFloat

Como dito anteriormente, existe uma busca por \(\pi\) com diversas casas decimais corretas. Você já deve ter percebido que o Julia traz uma implementação própria de \(\pi\), chamada através de pi. Essa aproximação pode ser vista com BigFloat também.

BigFloat(pi)
3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823537875937519577818577805321712268066130019278766111959092164201989380952572010654858632788659361533818279682303019520353018529689957736225994138912497217752834791315155748572424541506959508295331168617278558890750983817546374649393192550604009277016711390098488240128583616035637076601047101819429555961989467678374494482553797747268471040475346462080466842590694912933136770289891521047521620569660240580381501935112533824300355876402474964732639141992726042699227967823547816360093417216412199245863150302861829745557067498385054945885869269956909272107975093029553211653449872027559602364806654991198818347977535663698074265425278625518184175746728909777727938000816470600161452491921732172147723501414419735685481613611573525521334757418494684385233239073941433345477624168625189835694855620992192221842725502542568876717904946016534668049886272327917860857843838279679766814541009538837863609506800642251252051173929848960841284886269456042419652850222106611863067442786220391949450471237137869609563643719172874677646575739624138908658326459958133904780275900994657640789512694683983525957098258226205224894077267194782684826014769909026401363944374553050682034962524517493996514314298091906592509372216964615157098583874105978859597729754989301617539284681382686838689427741559918559252459539594310499725246808459872736446958486538367362226260991246080512438843904512441365497627807977156914359977001296160894416948685558484063534220722258284886481584560285060168427394522674676788952521385225499546667278239864565961163548862305774564980355936345681743241125150760694794510965960940252288797108931456691368672287489405601015033086179286809208747609178249385890097149096759852613655497818931297848216829989487226588048575640142704775551323796414515237462343645428584447952658678210511413547357395231134271661021359695362314429524849371871101457654035902799344037420073105785390621983874478084784896833214457138687519435064302184531910484810053706146806749192781911979399520614196634287544406437451237181921799983910159195618146751426912397489409071864942319615679452080951465502252316038819301420937621378559566389377870830390697920773467221825625996615014215030680384477345492026054146659252014974428507325186660021324340881907104863317346496514539057962685610055081066587969981635747363840525714591028970641401109712062804390397595156771577004203378699360072305587631763594218731251471205329281918261861258673215791984148488291644706095752706957220917567116722910981690915280173506712748583222871835209353965725121083579151369882091444210067510334671103141267111369908658516398315019701651511685171437657618351556508849099898599823873455283316355076479185358932261854896321329330898570642046752590709154814165498594616371802709819943099244889575712828905923233260972997120844335732654893823911932597463667305836041428138830320382490375898524374417029132765618093773444030707469211201913020330380197621101100449293215160842444859637669838952286847831235526582131449576857262433441893039686426243410773226978028073189154411010446823252716201052652272111660396665573092547110557853763466820653109896526918620564769312570586356620185581007293606598764861179104533488503461136576867532494416680396265797877185560845529654126654085306143444318586769751456614068007002378776591344017127494704205622305389945613140711270004078547332699390814546646458807972708266830634328587856983052358089330657574067954571637752542021149557615814002501262285941302164715509792592309907965473761255176567513575178296664547791745011299614890304639947132962107340437518957359614589019389713111790429782856475032031986915140287080859904801094121472213179476477726224142548545403321571853061422881375850430633217518297986622371721591607716692547487389866549494501146540628433663937900397692656721463853067360965712091807638327166416274888800786925602902284721040317211860820419000422966171196377921337575114959501566049631862947265473642523081770367515906735023507283540567040386743513622224771589150495309844489333096340878076932599397805419341447377441842631298608099888687413260472156951623965864573021631598193195167353812974167729478672422924654366800980676928238280689964004824354037014163149658979409243237896907069779422362508221688957383798623001593776471651228935786015881617557829735233446042815126272037343146531977774160319906655418763979293344195215413418994854447345673831624993419131814809277771038638773431772075456545322077709212019051660962804909263601975988281613323166636528619326686336062735676303544776280350450777235547105859548702790814356240145171806246436267945612753181340783303362542327839449753824372058353114771199260638133467768796959703098339130771098704085913374641442822772634659470474587847787201927715280731767907707157213444730605700733492436931138350493163128404251219256517980694113528013147013047816437885185290928545201165839341965621349143415956258658655705526904965209858033850722426482939728584783163057777560688876446248246857926039535277348030480290058760758251047470916439613626760449256274204208320856611906254543372131535958450687724602901618766795240616342522577195429162991930645537799140373404328752628889639958794757291746426357455254079091451357111369410911939325191076020825202618798531887705842972591677813149699009019211697173727847684726860849003377024242916513005005168323364350389517029893922334517220138128069650117844087451960121228599371623130171144484640903890644954440061986907548516026327505298349187407866808818338510228334508504860825039302133219715518430635455007668282949304137765527939751754613953984683393638304746119966538581538420568533862186725233402830871123282789212507712629463229563989898935821167456270102183564622013496715188190973038119800497340723961036854066431939509790190699639552453005450580685501956730229219139339185680344903982059551002263535361920419947455385938102343955449597783779023742161727111723643435439478221818528624085140066604433258885698670543154706965747458550332323342107301545940516553790686627333799585115625784322988273723198987571415957811196358330059408730681216028764962867446047746491599505497374256269010490377819868359381465741268049256487985561453723478673303904688383436346553794986419270563872931748723320837601123029911367938627089438799362016295154133714248928307220126901475466847653576164773794675200490757155527819653621323926406160136358155907422020203187277605277219005561484255518792530343513984425322341576233610642506390497500865627109535919465897514131034822769306247435363256916078154781811528436679570611086153315044521274739245449454236828860613408414863776700961207151249140430272538607648236341433462351897576645216413767969031495019108575984423919862916421939949072362346468441173940326591840443780513338945257423995082965912285085558215725031071257012668302402929525220118726767562204154205161841634847565169998116141010029960783869092916030288400269104140792886215078424516709087000699282120660418371806535567252532567532861291042487761825829765157959847035622262934860034158722980534989650226291748788202734209222245339856264766914905562842503912757710284027998066365825488926488025456610172967026640765590429099456815065265305371829412703369313785178609040708667114965583434347693385781711386455873678123014587687126603489139095620099393610310291616152881384379099042317473363948045759314931405297634757481193567091101377517210080315590248530906692037671922033229094334676851422144773793937517034436619910403375111735471918550464490263655128162288244625759163330391072253837421821408835086573917715096828874782656995995744906617583441375223970968340800535598491754173818839994469748676265516582765848358845314277568790029095170283529716344562129640435231176006651012412006597558512761785838292041974844236080071930457618932349229279650198751872127267507981255470958904556357921221033346697499235630254947802490114195212382815309114079073860251522742995818072471625916685451333123948049470791191532673430282441860414263639548000448002670496248201792896476697583183271314251702969234889627668440323260927524960357996469256504936818360900323809293459588970695365349406034021665443755890045632882250545255640564482465151875471196218443965825337543885690941130315095261793780029741207665147939425902989695946995565761218656196733786236256125216320862869222103274889218654364802296780705765615144632046927906821207388377814233562823608963208068222468012248261177185896381409183903673672220888321513755600372798394004152970028783076670944474560134556417254370906979396122571429894671543578468788614445812314593571984922528471605049221242470141214780573455105008019086996033027634787081081754501193071412233908663938339529425786905076431006383519834389341596131854347546495569781038293097164651438407007073604112373599843452251610507027056235266012764848308407611830130527932054274628654036036745328651057065874882256981579367897669742205750596834408697350201410206723585020072452256326513410559240190274216248439140371

Podemos modificar nossa implementação para buscar \(\pi\) com BigFloat também.

function pi_madhava_BF()
    S = BigFloat(1.0)
    third = BigFloat(1.0)/BigFloat(3.0)
    σ = -1
    pow_third = third
    t = σ * third * third
    n = 1
    while S + t != S
        S += t
        n += 1
        pow_third *= third
        σ = -σ
        t = σ * pow_third / (2n + 1)
        if n > 1_000_000
            break
        end
    end
    println("n = $n") # Usar para mostrar quantas iterações
    return sqrt(BigFloat(12.0)) * S
end

setprecision(2^10)
pi_madhava_BF() - pi
n = 641
0.0
setprecision(2^11)
round(Int, -log10(eps(BigFloat)))
616
for n = 10:14
    setprecision(2^n)

    E = max(abs(pi_madhava_BF() - BigFloat(pi)), eps(BigFloat))
    casas = round(Int, -log10(E))
    println("Com $(2^n) bits obtive $casas casas")
end
n = 641
Com 1024 bits obtive 308 casas
n = 1286
Com 2048 bits obtive 616 casas
n = 2578
Com 4096 bits obtive 1231 casas
n = 5161
Com 8192 bits obtive 2465 casas
n = 10329
Com 16384 bits obtive 4930 casas

setprecision(2^15)
@time pi_madhava_BF();
n = 20666
  1.050171 seconds (206.68 k allocations: 420.997 MiB, 2.97% gc time)

function logaritmo2()
    ln2 = 1.0 / 2.0
    k = 2
    r = 1 / 4
    while ln2 + r /k != ln2
        ln2 += r / k
        k += 1
        r = r / 2
        if k > 1_000_000
            break
        end
    end
    println("k = $k")
    return ln2
end

logaritmo2() - log(2)
k = 49
-2.220446049250313e-16

Exercícios

Faça os exercícios do capítulo 1 do livro Cálculo Numérico de Ruggiero e Lopes.

  1. Calcule a expansão de Taylor das seguintes funções, em torno do ponto a dado:

  • \(f(x) = x^2 - 5x + 6\), \(a = 2\);

  • \(g(x) = x^3 - 3x^2 + 3x - 1\), \(a = 1\);

  • \(h(x) = e^x\), \(a = 0\);

  • \(z(x) = \ln (x)\), \(a = 1\);

  1. Considere as seguintes séries e sequências convergentes:

  • \(2 = 1 + \frac{1}{2} + \frac{1}{4} + \frac{1}{8} + \dots + \frac{1}{2^n} + \dots\).

  • \(a_{n+1} = \frac{a_n}{2} + \frac{1}{a_n}\), com \(a_1 = 1\).

  • \(b_{n+1} = \sqrt{b_n + 1}\), com \(b_1 = 1\).

  • \(\phi_n = \dfrac{F_n}{F_{n-1}}\), onde \(F_n\) é o n-ésimo termo de Fibonacci.

  • Alguma série para \(\pi\) da internet.

Para cada uma delas, faça os seguintes itens

  • Implemente uma função que recebe um número \(n\) e calcula a soma parcial \(s_n\) se o item for uma série, ou n-ésimo elemento se o item for uma sequência. Faça sua função ser o mais econômica possível e evite overflows.

  • Faça um gráfico do erro pelo número de termos \(n\).

  • Dado \(\varepsilon > 0\), é possível determinar n de modo que o erro seja menor que \(\varepsilon\)? Veja quantos casos você consegue fazer isso, ou encontrar algum limitante desse tipo.