A Equação da Torre Eiffel

No mês de maio eu completei cinco anos de casado! Quando eu casei, aproveitei a lua-de-mel em Paris para falar sobre a ciência da Torre Eiffel, então nada mais apropriado que aproveitar o aniversário para falar sobre outro mistério da torre: a equação que descreve o seu formato!

   

A torre mais alta do mundo

O desafio proposto pelo Gustave Eiffel era grandioso: construir a torre mais alta do mundo. Na época, o recorde era o Monumento de Washington, com 169 metros. O Eiffel pretendia erguer sua torre com 324 metros, quase o dobro do recorde anterior.

Para entender o desafio, vale a pena comparar a Torre Eiffel com outras torres da antiguidade. O desenho abaixo está em escala, e mostra a Torre Eiffel comparada com a Grande Pirâmide de Giza no Egito, e com o Campanário de São Marcos em Veneza:

Imagem: Wikipedia (CC-SA)

O primeiro desafio ao fazer uma torre é garantir que ela consiga suportar seu próprio peso. Os egípcios resolveram isso com força bruta: usaram blocos de pedra maciços e fizeram a base bem maior que o topo. O resultado final certamente funcionou, a pirâmide está aí de pé 4500 anos depois. Mas ela é um desperdício de material, com tecnologia moderna você não precisa de tanta matéria-prima para chegar nessa altura.

O Campanário de São Marcos é o extremo oposto. Ele foi construído em 1514, e nessa época a tecnologia já permitia criar uma torre com a largura da base igual à do topo. O problema é que a torre caiu! Em 1902 ela não resistiu ao próprio peso e colapsou (nenhuma vida humana foi perdida, mas o gato do zelador morreu). A torre que você visita hoje em Veneza é uma reconstrução.

Como fazer, então, uma torre alta que suporte o próprio peso? O Eiffel sabia a resposta porque tinha experiência no assunto. Antes de criar os ícones pelo qual é mais lembrado, ele construía pontes. Existem pontes do Eiffel por todo o mundo, desde a França e a Romênia, até o Vietnã e o Chile.

Para entender a tecnologia do Eiffel, vale a pena comparar com uma ponte antiga: a Pont du Gard, um aqueduto romano construído dois mil anos atrás. Ele cobre uma distância de 275 metros a uma altura de 48 metros. Para isso, ela precisa de cinco apoios no chão e três níveis de arcos para aguentar seu próprio peso:

Imagem: Wikipedia (CC-SA)

Agora, compare com uma das mais belas pontes construídas pelo Eiffel: a Ponte Dona Maria, na Cidade do Porto, em Portugal. Ela cobre 353m a 60m de altura, em ambos os casos maior que a ponte romana. E faz isso usando só dois apoios no chão! O segredo do Eiffel para conseguir esse feito é a magia da treliça. Ao usar uma treliça, você consegue suportar mais peso com menos material.

Imagem: Wikipedia (CC-SA)

Todo mundo sabia que o Eiffel era expert no assunto, e que ele saberia construir uma torre que suportasse o próprio peso. Mas, ainda assim, na época ele ouviu muitas críticas. Uma coisa é você construir uma ponte na horizontal, outra é construir uma torre na vertical. Estruturas na vertical tem um problema adicional: será que o Eiffel sabia como fazer a torre resistir ao vento?

O método geométrico do Eiffel

Quando a torre ficou pronta, o Eiffel calou os críticos. Quem já visitou a torre sabe que no terceiro piso venta. Venta muito mesmo! Porém, a torre não se abala. Fica lá, imóvel, mesmo com os ventos mais fortes.

E qual foi o truque do Eiffel para deixar a torre tão estável com o vento? Até pouco tempo atrás, ninguém sabia. As plantas da torre Eiffel são conhecidas, você pode ver todas as especificações em um livrão publicado pela Taschen (todas mesmo, tem até o diâmetro e comprimento de cada parafuso). Mas as plantas não falam o truque da estabilidade, tudo que dizem é que o Eiffel determinou o formato da torre de forma a minimizar a ação do vento, e fez isso usando "métodos geométricos baseados em observações empíricas".

  

O formato da torre foi um mistério por mais de um século, até que Patrick Weidman, professor da Universidade de Colorado, matou o problema, e da maneira mais simples possível. Ele foi atrás dos diários do Eiffel e achou lá a dica que resolve o problema! O trecho do diário dizia:

"Suponha que as faces de uma treliça simples formam uma parede resistente às forças cisalhantes do vento, cujas componentes horizontais são P', P'', P''', P''''. Sabemos que, para calcular as forças atuantes no corte pelo plano MN, precisamos determinar a resultante P de todas as forças exteriores atuando acima da secção, e decompor essa resultante nas três forças passando pelas peças cortadas. Se a forma do sistema for tal que, para cada corte horizontal MN, a extensão das peças externas da treliça se intersectam na força exterior P, então as forças na barra interna serão zero e podemos excluir esse membro."

O diário tinha esse desenho acompanhando:

Lendo com cuidado dá para entender direitinho a motivação do formato da torre. Vamos recriar o raciocínio em câmera lenta. Nós começamos supondo que o vento é uma carga horizontal atuando sobre toda a torre (que ele representou como as forças P', P'', etc). O truque agora é pegar um corte horizontal arbitrário (que ele chamou de MN), e tirar toda a parte de cima da torre. As forças na estrutura não mudam se substituirmos a parte de cima por um ponto material de mesma massa, localizado no centróide do que foi tirado.

    

Vamos então ver qual é o efeito do vento. Eu tenho uma força P horizontal, que é a resultante do vento. Tenho ainda três forças F1, F2, F3, que são as forças atuando sobre as vigas da treliça. Quando uma treliça está em equilíbrio, as forças sobre ela são sempre de tração ou compressão, ou seja, estão na mesma direção de suas vigas respectivas.

A torre está em equilíbrio estático, logo a soma das forças é zero. Além disso, ela também não está rotacionando; logo a soma dos torques também é zero. O torque é o produto da força pela distância; mas não é qualquer distância, é a distância na direção perpendicular à força.

Agora o truque do Eiffel é claro. Ele construiu as laterais da torre de modo que as tangentes às vigas mais externas sempre passam pelo centróide. Quando você faz isso, a direção das forças P, F1 e F3 apontam para o centróide, então o torque delas é zero. Só a força F2 está aplicando torque no centróide!

Podemos calcular a soma dos torques agora. A distância associada à força F2 é a distância d, em verde na figura. Note que ela é estritamente positiva.
$$\begin{align*} P\times 0 +F_1\times 0+F_2\times d+F_3\times 0&=0\\ F_2\times d&=0\\ F_2&=0 \end{align*}$$
Ahá! Se as tangentes externas apontam para o centróide, então a força na viga interna é zero! Se o vento horizontal fosse a única força atuando sobre a torre, então nós poderíamos até tirar a viga interna que não faria diferença!

É claro que na prática não podemos tirar a viga. Ela não ajuda a resistir ao vento, mas ajuda a sustentar o peso. Por outro lado, essa viga pode ser mais fraca que as outras. Com a idéia do Eiffel, só as vigas externas precisam ser fortes. Nessa foto que eu tirei do pilar oeste isso é bem visível: as vigas externas são largas e feitas de material rígido, enquanto que as internas são finas e feitas de mini-treliças, o que deixa a estrutura como um todo mais leve. E quanto mais leve, mais alto você pode subir.

A equação da Torre Eiffel

Agora que nós sabemos o truque do Eiffel, podemos calcular qual é o lugar geométrico dos pontos cujas tangentes passam pelo centróide, isso vai nos dar a equação que descreve o fomato da torre. Aparentemente o Eiffel nunca precisou fazer isso, já que ele usou "um método geométrico baseado em observações empíricas" (para mim isso é um jeito de dizer que ele calculou uma solução numérica ao invés de resolver a equação de fato).

Mas ele poderia ter resolvido a equação se quisesse. Vamos fazer as contas na caixa azul para conferir, pule se você tem medo de equações diferenciais:

Começaremos supondo que o formato da torre pode ser descrito por uma função f(x). Para simplificar as contas (e, acredite, simplifica mesmo), nós vamos rotacionar o gráfico de modo a colocar o eixo x na vertical, orientado para baixo, e o f(x) na horizontal. Vamos ainda supor que f(x) é real, contínua e pode ser diferenciada muitas vezes.
No gráfico acima, f(x) está desenhada em azul. Por um ponto qualquer x eu traço as retas tangentes (em vermelho), elas vão se encontrar em um ponto xC, que nós estamos supondo que sempre será o centróide da parte superior da torre (ou seja, do trecho que vai do ponto escolhido x, até o topo da torre xT). A posição do centróide é dado diretamente pela definição: $$x_C=\frac{\int_{x_T}^{x}t f(t)\,dt}{\int_{x_T}^{x} f(t)\,dt}$$ A tangente podemos calcular de duas maneiras. Nós sabemos que a tangente que passa por x é dada diretamente por f'(x). Por outro lado, ela também é a tangente do ângulo formado pelo eixo x e pela reta vermelha, que você calcula direto como cateto oposto pelo adjacente: $$f'(x)=\frac{\Delta y}{\Delta x}=\frac{f(x) -0}{x-x_C}=\frac{f(x)}{x-x_C}$$ Podemos isolar xC nessa última equação: $$\begin{align*} x-x_C&=\frac{f(x)}{f'(x)}\\ x_C&=\frac{x f'(x)-f(x)}{f'(x)} \end{align*}$$ E depois igualar as duas equações: $$\begin{align*} \frac{x f'(x)-f(x)}{f'(x)} &= \frac{\int_{x_T}^{x}t f(t)\, dt}{\int_{x_T}^{x} f(t)\, dt}\\ \left(x f'(x)-f(x)\right)\int_{x_T}^{x} f(t)\, dt &= f'(x)\int_{x_T}^{x}t f(t)\, dt\\ \left(f'(x)\int_{x_T}^{x}x f(t)\, dt \right)-\left(f(x)\int_{x_T}^{x} f(t)\, dt\right) &= f'(x)\int_{x_T}^{x}t f(t)\, dt\\ f(x)\int_{x_T}^{x} f(t)\, dt &= f'(x)\int_{x_T}^{x}(x-t) f(t)\, dt \end{align*}$$ Na última linha, veja que o x passou para dentro da integral, você pode fazer isso porque ele é constante em relação à variável de integração, que é t. Antes de prosseguir, vamos dar nomes às integrais, chamando-as de g(x) e h(x). Note que essas integrais são funções de x. $$\begin{align*} g(x)&=\int_{x_T}^{x} f(t)\, dt \\ h(x)&= \int_{x_T}^{x}(x-t) f(t)\, dt \end{align*}$$ Vamos substituir as integrais pelas funções que nomeamos, e depois derivar dos dois lados. Quando você deriva dos dois lados, você está comendo uma constante, então no final podem aparecer soluções espúrias. Depois de acabar tudo, temos que pegar as soluções e voltar para o começo, conferindo se elas são realmente válidas. $$\begin{align*} f(x)g(x)&=f'(x)h(x)\\ f'(x)g(x)+f(x)g'(x) &= f''(x)h(x)+f'(x)h'(x) \end{align*}$$ Como calcular g'(x) e h'(x)? Vamos começar por g'(x). Note, primeiro, que você pode usar o teorema fundamental do cálculo para escrever g(x) em função da primitiva de f(x), que iremos chamar de F(x): $$g(x)=\int_{x_T}^x f(t)\,dt=F(x) -F(x_T) $$ Já que F(xT) é uma constante, sua derivada é zero. F(x) é uma primitiva, e a derivada da primitiva de uma função é a própria função (desde que ela seja contínua): $$\begin{align*} g(x)&=F(x) -F(x_T)\\ g'(x)&=F'(x)-0=f(x) \end{align*} $$ Então g'(x)=f(x), essa foi fácil. Para calcular h'(x) é mais complicado. Lembremos a regra de integração por partes para integrais definidas: $$\begin{align*} \int_c^d a(t)b'(t)\,dt&=[a(t)b(t)]_c^d-\int_c^d a'(t)b(t)\,dt \end{align*} $$ Vamos escolher quem são a(t) e b(t): $$\begin{align*} a(t)=x-t&\implies a'(t)=-1\\ b(t)=F(t)&\implies b'(t)=f(t)\\ \end{align*} $$ Agora podemos quebrar h(x) em duas partes: $$\begin{align*} h(x)&=\int_{x_T}^x(x-t)f(t)\,dt\\ &=[(x-t)F(t)]_{x_T}^x-\int_{x_T}^x-F(t)\,dt\\ &=(x-x)F(x)-(x-x_T)F(x_T)+\int_{x_T}^xF(t)\,dt\\ &=-(x-x_T)F(x_T)+\int_{x_T}^xF(t)\,dt\\ &=-(x-x_T)F(x_T)+P(x)-P(x_T)\\ \end{align*} $$ No último passo, por falta de notação melhor, eu chamei de P(x) a primitiva de F(x). Agora podemos derivar h(x). Note que, se P(x) é a primitiva de F(x), então F(x) é a derivada de P(x). Além disso, como P(xT) é constante, sua derivada é zero: $$\begin{align*} h(x)&=-(x-x_T)F(x_T)+P(x)-P(x_T)\\ h'(x)&=-F(x_T)+F(x)-0\\ h'(x)&=\int_{x_T}^xf(t)\,dt\\ h'(x)&=g(x)\\ \end{align*} $$ Concluímos então que h'(x) é o mesmo g(x) que tínhamos visto antes! Antes de continuar, vale a pena resumir o que nós temos até agora: $$\begin{align*} f(x)g(x)&=f'(x)h(x)\\ f'(x)g(x)+f(x)g'(x) &= f''(x)h(x)+f'(x)h'(x)\\ f(x)&=g'(x)\\ f'(x)&=g''(x)\\ f''(x)&=g'''(x)\\ h'(x)&=g(x) \end{align*}$$ Note que dá para reescrever quase tudo usando só g(x) e suas derivadas. Fica faltando h(x), mas isso não é problema porque temos duas equações usando h(x), então é só isolar. Da primeira equação temos: $$\begin{align*} f(x)g(x)&=f'(x)h(x)\\ h(x)&=\frac{f(x)g(x)}{f'(x)}\\ h(x)&=\frac{g'(x)g(x)}{g''(x)}\\ \end{align*}$$ Da segunda equação, temos: $$\begin{align*} f'(x)g(x)+f(x)g'(x)&=f''(x)h(x)+f'(x)h'(x)\\ f''(x)h(x)&=f'(x)g(x)+f(x)g'(x)-f'(x)h'(x)\\ h(x)&=\frac{f'(x)g(x)+f(x)g'(x)-f'(x)h'(x)}{f''(x)}\\ h(x)&=\frac{g''(x)g(x)+g'(x)g'(x)-g''(x)g(x)}{g'''(x)}\\ h(x)&=\frac{g'(x)^2}{g'''(x)}\\ \end{align*}$$ Igualando as duas expressões para h(x): $$\begin{align*} \frac{g'(x)g(x)}{g''(x)} &= \frac{g'(x)^2}{g'''(x)} \\ \frac{g'(x)}{g(x)} &= \frac{g'''(x)}{g''(x)} \end{align*}$$ Tudo que falta agora é resolver a equação para g(x) e finalmente temos o formato da torre. Mas essa equação é uma equação diferencial não-linear de terceira ordem. E isso é os que os matemáticos usam para assustar os filhos: "se você não comer, vai aparecer o bicho-papão, o homem do saco e a equação diferencial não-linear de terceira ordem para te pegar!" Felizmente, essa equação pode ser resolvida se você tiver a visão além do alcance. O que acontece quando você tenta diferenciar log g(x)? $$\begin{align*} \frac{d}{dx}\log(g(x))=\frac{1}{g(x)}\times g'(x)=\frac{g'(x)}{g(x)} \end{align*}$$ Ahá! Essa é exatamente a forma da equação. Fazendo o análogo no outro lado, podemos reescrever a equação original: $$\begin{align*} \frac{g'(x)}{g(x)} &= \frac{g'''(x)}{g''(x)}\\ \frac{d}{dx}\log g(x) &= \frac{d}{dx}\log g''(x) \end{align*}$$ Agora nós podemos integrar dos dois lados (lembrando que vai aparecer uma constante), e depois tirar a exponencial dos dois lados: $$\begin{align*} \frac{d}{dx}\log g(x) &= \frac{d}{dx}\log g''(x)\\ \log g(x) &= \log g''(x) + C\\ g(x) &= k g''(x) \\ \end{align*}$$ Agora ficou fácil né? A equação medonha virou uma equação diferencial linear de segunda ordem, essa é bem-comportada, dá a patinha e finge de morta. A solução geral é da forma abaixo: $$g(x)=A e^{B x}+C e^{-B x}$$ Como f(x)=g'(x), então f(x) também tem a forma acima, só os coeficientes que serão diferentes. A natureza da função depende dos coeficientes; por exemplo, se A=0.5, C=0.5 e B=i, então f(x)=cos(x). Mas nem todas as combinações de coeficientes são soluções válidas! Para determinar os coeficientes corretos, precisamos voltar o f(x) lá na equação original: $$f(x)\int_{x_T}^{x} f(t)\, dt = f'(x)\int_{x_T}^{x}(x-t) f(t)\, dt $$ Daqui em diante é só fazer as contas. Se você substituir tudo, vai notar que o único jeito de conseguir uma solução real válida para qualquer x é fazendo C=0 e xT=-∞; nesse caso A pode ser um real qualquer, e B pode ser um real positivo qualquer. Estritamente falando, a torre teria que ser infinitamente alta para satisfazer a essa equação, mas na prática você pode aproximar "infinitamente alta" por "muito, muito alta". Concluindo, a solução final para o formato da torre Eiffel é: $$f(x)=A e^{B x} $$ Ou seja, a torre tem um perfil exponencial!

Matemática e Engenharia

Agora que sabemos qual o formato teórico da torre, podemos usar os dados reais da torre para estimar os coeficientes da exponencial. Eu rodei um best fit e cheguei em A=3.34323 e B=0.0102592:



A aproximação é boa, mas... parece que não encaixa direitinho? É meio frustrante fazer esse monte de contas, e no fim a curva real não ser tão parecida assim com a equação deduzida.

Mas tem um motivo para isso, e o motivo é que o Eiffel era engenheiro, não matemático. Quando você plota uma exponencial em um gráfico semi-log, o resultado tem que ser uma linha. Mas olha só o que acontece com os dados reais da torre: eles formam duas linhas!

 

Coisas que você aprende na sua primeira aula de engenharia: se você faz o projeto de um elevador, e determina que o cabo do elevador precisa aguentar no mínimo 10 toneladas, qual o cabo que você coloca no elevador? Um matemático diria que é um cabo de 10 toneladas, mas o engenheiro vai falar que o mínimo são 14 toneladas.

Isso é o resultado do fator de segurança. Por mais que o seu modelo matemático seja correto, sempre tem alguma coisa que você não levou em conta (pode ser que o aço que você usou tenha impurezas, ou que o solo não é tão firme quanto você achava, e assim por diante). Por isso, é costume sempre multiplicar o valor final por um fator de segurança para levar em contas esses imprevistos.

O Eiffel sabia disso, e no diário ele explica que imaginava que a base da torre estaria sujeita a torques maiores que o topo, e por isso usou um fator de segurança diferente na base e no topo. Quando você faz o best fit com duas exponenciais ao invés de uma só, o gráfico fica bem melhor!

 

Agora finalmente podemos dar o mistério por resolvido: o formato da Torre Eiffel é descrito por duas exponenciais, escolhidas para minimizar o material necessário para resistir à carga do vento.

A Torre Eiffel é o literalmente o maior monumento à ciência construído, da próxima vez que passar por Paris aproveite para fazer uma reverência aos cientistas do passado :)

Matemática x Tanques de Guerra

Quando eu era criança pequena lá na década de 80, eu tinha um brinquedo de tanque de guerra. Ele chamava Panzer, e era fabricado pela brinquedos Mimo. Para a época era bem divertido: tinha controle remoto, e, além de andar, ainda disparava mísseis!

Mas eu não percebia a ironia do brinquedo. O Panzer não era só um brinquedo, ele existiu de verdade. Na vida real, o Panzer era uma máquina nazista de destruição. A série Panzer de tanques alemães foi responsável pela morte direta de incontáveis aliados. Ver um Panzer surgindo no horizonte era como ver a morte vindo em sua direção. Era um mimo esse Panzer.

Por outro lado, o Panzer era forte, mas não era invencível. Os aliados tinham armas capazes de deter o Panzer, o problema era descobrir quantas dessas armas eles deveriam levar para o front. Poucas semanas antes do Dia D, os alemães começaram a usar o novo modelo Panzer V, e os aliados não tinham como estimar quantas armas levar sem saber quantos Panzer V os inimigos fabricaram.

Foi então que os aliados usaram a arma definitiva contra o Panzer: a Matemática! Usando um truque muito esperto, eles estimaram que os alemães estavam fabricando 270 Panzers por mês. Depois que a guerra acabou, eles foram nos arquivos alemães conferir os números, e na verdade estavam fazendo 276 Panzers por mês. A estimativa foi muito boa!

A Quantidade de Tanques

O truque dos aliados se baseava em um erro dos alemães e em uma hipótese estatística. O erro dos alemães era que os tanques tinham número de série sequencial. Quando saíam da fabrica, eram numerados como 1, 2, 3, e assim por diante. Ao fazer isso, eles deixaram uma brecha para os hackers aliados.

Suponha que você capturou um Panzer alemão, e o serial dele era 14. Você consegue estimar quantos tanques eles produziram baseado só nesse número? No caso geral não, mas você pode fazer uma hipótese simplificadora, que é supor que o tanque que você capturou é aleatório (ou seja, os alemães estavam distribuindo igualmente os tanques por toda a Europa).

Com essa hipótese a intuição é clara. Vamos supor, hipoteticamente, que os alemães fizeram 200 tanques. Qual chance de você capturar um tanque de serial menor ou igual a 14 nessa situação? É 14/200 né? E se os alemães fizeram só 20 tanques? Aí a chance é 14/20, dez vezes maior que 14/200. Com a hipótese de tanque aleatório, é mais provável que os alemães tenham feito 20 tanques do que 200 tanques.

Você pode formalizar esse argumento usando o teorema de Bayes, e calcular o valor esperado do número de tanques. Vamos fazer as contas na caixa azul:

Suponha que conseguimos capturar k tanques. Dentre esses k tanques, o maior serial encontrado tem o valor m. Com essas informações, vamos calcular qual é o valor esperado do número de tanques produzidos, que vamos chamar de n. Podemos abrir esse valor esperado pela definição:
$$E[n\;|\;m,k] = \sum_{n} n \; p(n\; |\; m,k)$$
O problema então é achar a probabilidade condicional de n, dados m e k. Nós ainda não sabemos quanto vale essa probabilidade, mas podemos abri-la usando o teorema de Bayes:
$$p(n\;|\;m,k)=\frac{p(m,k\;|\;n)\;p(n)}{p(m,k)}$$
Essas três probabilidades nós conseguimos calcular! A primeira delas, p(m,k|n), sai por combinatória. Dado um vetor de tamanho n, de quantas maneiras conseguimos escolher k elementos, de modo que o maior deles seja m?

Um dos elementos precisa necessariamente ser o m, então só precisamos calcular a posição dos outros k-1 restantes. E eles precisam ser menores que m, então só podem estar nas m-1 posições menores que m. Por isso, essa quantidade vale binomial(m-1, k-1). Note, entretanto, que isso só vale quando n>=m, caso contrário a probabilidade é zero (não tem como o máximo ser m se o vetor tem um tamanho menor que m).

Para calcular a probabilidade desejada, temos que achar de quantas maneiras podemos escolher k elementos dentre n, sem considerar a restrição do maior deles ser m. Mas isso é o próprio binomial(n, k). Então, a probabilidade é:
$$p(m,k\;|\;n) = \frac{\displaystyle{m-1\choose k-1}} {\displaystyle{n\choose k}}[n\ge m]$$
Vamos agora para a segunda, p(n). Eu não tenho nenhuma informação sobre quantos tanques foram produzidos, então não sei qual é a probabilidade de ter n tanques. Mas eu posso chutar que a distribuição é uniforme. A chance de ter 15 tanques, ou de ter 500 tanques, a priori, deve ser a mesma. Por isso, vamos fazer p(n) ser uma constante independente de n:
$$p(n)=c$$
A última, p(m,k), é a probabilidade do máximo ser m em k escolhas, independente do valor de n. Para isso, podemos usar a lei da probabilidade total: 
$$p(m,k) = \sum_n p(m,k\;|\;n) \;p(n)$$ 
Essas duas já sabemos calcular! É só substituir:
 $$\begin{align*} p(m,k) &= \sum_n p(m,k\;|\;n) \;p(n) \\ &= \sum_n \frac{{m-1\choose k-1}} {{n\choose k}}[n\ge m]\;c\\ &= \sum_{n\ge m} c {m-1\choose k-1} {n\choose k}^{-1}\\ &= c {m-1 \choose k-1} \sum_{n\ge m} {n\choose k}^{-1}\\ \end{align*} $$
A somatória tem uma forma fechada, que você pode achar com o algoritmo de Gosper:
$$\sum_{n=m}^{\infty}{n\choose k}^{-1}=\frac{m}{k-1}{m\choose k}^{-1}$$ 
Substituindo:
$$\begin{align*} p(m,k) &= c {m-1 \choose k-1} \sum_{n\ge m} {n\choose k}^{-1}\\ &= c {m-1 \choose k-1} \times \frac{m}{k-1}{m\choose k}^{-1} \\ &= \frac{cm}{k-1} {m-1 \choose k-1} {m\choose k}^{-1} \\ \end{align*} $$
Agora podemos voltar e substituir as três probalidades na fórmula de Bayes:
$$ \begin{align*} p(n\;|\;m,k)&=\frac{p(m,k\;|\;n)\;p(n)}{p(m,k)} \\ &= \frac{\displaystyle{m-1\choose k-1}{n\choose k}^{-1}[n\ge m]\times c} {\displaystyle\frac{cm}{k-1}{m-1\choose k-1}{m\choose k}^{-1}}\\ &= \frac{k-1}{m}{m\choose k}{n\choose k}^{-1}[n\ge m] \end{align*}$$
Por fim, podemos substituir na fórmula do valor esperado:
 $$\begin{align*} E[n\;|\;m,k] &= \sum_{n} n \; p(n\; |\; m,k)\\ &= \sum_n n \frac{k-1}{m}{m\choose k}{n\choose k}^{-1}[n\ge m]\\ &= \frac{k-1}{m}{m\choose k}\sum_{n\ge m} n {n\choose k}^{-1}\\ \end{align*}$$
Essa somatória também tem uma forma fechada por Gosper:
$$ \sum_{n=m}^{\infty} n{n\choose k}^{-1}=\frac{m(m-1)}{k-2}{m\choose k}^{-1}$$
Chegamos então na substituição final:
$$\begin{align*} E[n\;|\;m,k] &= \frac{k-1}{m}{m\choose k}\sum_{n\ge m} n {n\choose k}^{-1}\\ &= \frac{k-1}{m}{m\choose k}\times\frac{m(m-1)}{k-2}{m\choose k}^{-1} \\ &= \frac{(k-1)(m-1)}{k-2} \\ \end{align*}$$

A Estimativa na Prática

A fórmula final é super simples, mas ela funciona mesmo na prática? Podemos testá-la fazendo uma simulação numérica. Para isso, eu fiz dez mil simulações. Em cada uma eu sorteio o número de tanques fabricados, e escolho aleatoriamente 3 deles. Aí eu calculo a estimativa pela fórmula, e normalizo o erro. O código está no github, e o histograma resultante é o abaixo. A normalização é tal que o valor 0 é uma estimativa totalmente correta:

A estimativa não ficou muito boa não! Ao invés de ter uma gaussiana em torno do zero, o estimador tem um bias muito forte para os negativos. Felizmente, esse problema só acontece porque capturamos apenas 3 tanques. Se, ao invés disso, tivéssemos capturado 30 tanques, então o histograma seria bem melhor. Compare o histograma de 30 tanques sobreposto ao de 3 tanques:

Bem melhor né? Agora quase todas as estimativas estão bem próximas do zero. Esse método é bastante sensível ao número de tanques capturados. Os aliados sabiam disso, por isso que a estimativa deles foi tão boa (estimaram 270 tanques/mês, quando o valor real era 276). Sabe quantos tanques eles capturaram para fazer essa estimativa? Só dois tanques!

Claro que tem um truque. A estimativa com dois tanques é muito ruim, mas eles notaram que os alemães colocavam serial em todas as peças do tanque. Em especial, cada tanque tinha 48 rodas, cada uma com um serial único. Por isso, eles conseguiram usar a fórmula com k=96, o que deu a precisão alta que queriam. E no fim a matemática ganhou a guerra :)

O Objetivo Final da Matemática

Quem acompanhou as contas na caixa azul do post sobre o Pokèmon Twitch deve ter notado que eu pulei um pedaço. Em certo ponto da demonstração, eu magicamente tirei do chapéu uma fórmula nada óbvia:

E como foi que eu cheguei nessa fórmula? Simples, eu calculei no Wolfram Alpha:

"Ricbit, você trapaceou! Que roubalheira! Um matemático de verdade jamais faria isso!"

Mas será que é trapaça mesmo? Nesse post, eu vou mostrar como se resolve essa conta sem usar o Wolfram Alpha; e, se você acha que é trapaça, talvez acabe mudando de idéia :)

A progressão hipergeométrica

Para resolver essa somatória sem usar o Wolfram Alpha, nós precisamos relembrar as progressões que aprendemos no colégio. A primeira delas foi a progressão aritmética, onde a diferença entre dois números consecutivos é constante:

   

A segunda delas é a progressão geométrica, onde a razão entre dois números consecutivos é constante:

   

Uma progressão extremamente útil que não ensinam no colégio é a progressão hipergeométrica. Ela parece com a progressão geométrica, porque é definida a partir da razão entre dois elementos consecutivos. Mas agora, ao invés dessa razão ser constante, nós permitimos que a razão seja uma função racional (que é só um jeito chique de dizer que ela é o quociente entre dois polinômios):

 

Naquela somatória que queremos resolver, os termos sendo somados formam uma progressão hipergeométrica. Confira:

   

O numerador é um polinômio em q, o denominador é um polinômio em q. Então confere com a nossa definição de progressão hipergeométrica.

A soma da progressão hipergeométrica

Certamente vocês lembram que as progressões do colégio tinham uma fórmula para a soma de seus termos. A soma dos n primeiros termos da progressão aritmética é:

A soma dos n primeiros termos da progressão geométrica é:

 

A soma dos n primeiros termos da progressão hipergeométrica é... oops! Não existe uma fórmula geral para a soma de termos hipergeométricos.

Antigamente, achar uma fórmula para uma soma hipergeométrica era tão raro que elas até ganhavam nomes especiais. Por exemplo, a igualdade abaixo tem um nome, é o Binômio de Newton:

 

Outras somas hipergeométricas não tem fórmula, mas são tão úteis, que a própria soma ganha um nome e passa a valer como se fosse uma fórmula:

Por muitos séculos esse foi o estado da arte. Alguém acha uma fórmula aqui, alguém nomeia uma soma ali. Mas isso mudou no meio do século 20.

O algoritmo de Gosper

O primeiro progresso na direção de um método geral de resolução de somas hipergeométricas foi feito pela Mary Celine Fasenmyer em 1945. Depois de se formar em matemática, ela se juntou à ordem religiosa Sisters of Mercy (sim, a mesma que inspirou o nome da banda), e lá ela criou o Método da Irmã Celine, um algoritmo que consegue resolver automaticamente alguns tipos de somas hipergeométicas.

Em 1978 o método foi aprimorado por Bill Gosper e se tornou o algoritmo de Gosper (muito embora eu desconfie que, se a Irmã Celine fosse homem, o algoritmo seria de Gosper-Celine). O método de Gosper se baseia em um teorema surpreendente: A solução de uma soma hipergeométrica é outro hipergeométrico, se e somente se ela puder ser escrita na forma abaixo, onde R[j] é uma função racional.

Essa somatória parece esquisita. Se j é o índice da somatória, então o que ele está fazendo no lado direito da equação? E cadê os limites da somatória?

A resposta é que essa é uma somatória indefinida. Funciona do mesmo jeito que no cálculo: lá você tem integrais definidas e integrais indefinidas, aqui nós temos somatórias definidas (com limites), e somatórias indefinidas (sem limites). Para uma somatória indefinida, sempre vale a propriedade abaixo:

 

Essas duas equações são tudo que nós precisamos para usar o algoritmo de Gosper. Vamos à caixa azul calcular a nossa somatória original:

Começamos substituindo uma equação na outra:

Podemos então dividir tudo por h[j]:

 

Opa, a razão entre os h nós já calculamos lá em cima para a nossa somatória!

 

Falta achar R[j]. Eu sei que ele é um quociente de polinômios, mas não sei qual o grau desses polinômios. Eu vou chutar que o grau é 1, se der errado depois eu volto e aumento o grau. Pois bem, se o grau é 1, então R[j] é da forma abaixo:

 

Agora eu vou substituir. Aperte o cinto que lá vem turbulência:

   

O resultado é um polinômio gigante em j que é identicamente nulo. Se ele é identicamente nulo, então cada um dos seus coeficientes precisa ser zero. Com isso podemos montar um sistema de equações:

 

Agora é "só" resolver o sistema! Vou poupá-los de páginas e páginas, o resultado é o seguinte:

   

Se o sistema não tivesse solução, eu teria que voltar lá em cima e aumentar o grau do polinômio (se as contas foram titânicas com grau 1, imagine com grau 2). Mas o processo não é infinito, existem métodos que permitem achar um limite superior para o grau do polinômio. Se você atingir esse limite sem achar nenhuma solução, então você provou que sua somatória não pode ser escrita como um hipergeométrico.

Voltando à nossa equação original, agora nós sabemos quanto vale a somatória indefinida:

   

Nós queremos a soma de 1 até m, então precisamos converter a soma indefinida em uma soma definida, usando o teorema fundamental do cálculo discreto:

 

Aplicando o teorema, temos:

   

Pronto, esse é o resultado final. Eu acho muito importante que a pessoa resolva o algoritmo de Gosper uma vez na vida, para nunca mais cometer a sandice de repetir o processo. O método é mecânico, você só precisa seguir as regras e não errar as contas. O computador faz isso melhor que você. Use o Wolfram Alpha e economize grafite.

O objetivo final da Matemática

No final do século 19, o Bertrand Russell e o Alfred Whitehead tiveram uma idéia. Ele bolaram um plano para axiomatizar toda a matemática, de uma maneira completa e livre de contradições. Se o plano funcionasse, então toda a Matemática seria reduzida à manipulação de símbolos, que é um processo puramente mecânico e automatizável. O plano começou muito bem, a ponto do Whitehead fazer a seguinte declaração:

"A Civilização avança estendendo o número de operações importantes que podemos fazer sem ter que pensar."

Baseado nessa citação, o Concrete Mathematics do Knuth vai além e conclui:

"O objetivo final da matemática é eliminar a necessidade de todo pensamento inteligente."

É verdade, o objetivo sempre foi esse. Para os babilônios, resolver uma equação de segundo grau era dificílimo, só os grandes mestres conseguiam. Hoje em dia, nós transformamos a equação de segundo grau em uma fórmula que é ensinada para crianças de 14 anos. Você não precisa ser inteligente para usar a fórmula, é só substituir os números e fazer as contas.

Imagine se fosse possível fazer isso com toda a matemática, ao invés de só com equações de segundo grau! Você poderia usar computadores para fazer as contas, no lugar de crianças de 14 anos. A Matemática seria um problema resolvido, e aí você pode se concentrar no que interessa, que são as aplicações. Um físico acabou de descobrir uma propriedade quântica nova, será que ela permite criar um aparelho de teletransporte? Eu monto a equação, jogo no computador, e tcharam! Está aqui a resposta!

Toda a saga por trás do plano do Russell e do Whitehead é contada em quadrinhos no fabuloso Logicomix. Infelizmente, o final da história é que o plano falhou. E falhou espetacularmente. Em 1931, o Gödel provou que nenhum sistema axiomático é capaz de provar todas as sentenças verdadeiras que é capaz de descrever, acabando de uma vez por todas com o sonhado plano da axiomatização completa.

E o século 20 desceu ladeira abaixo. O Church provou que, mesmo se a axiomatização completa tivesse tido sucesso, seria impossível criar um método que prove um teorema qualquer a partir dos axiomas. O Turing provou que é impossível criar um método que analise um programa de computador e diga se ele termina ou se entra em loop infinito. O Matiyasevich provou que é impossível criar um método que resolva todas as equações diofantinas. O Richardson provou que é impossível criar um método que prove todas as identidades trigonométricas. Foi o fim da matemática automatizada... ou não?

Claro que não! Todos esses teoremas provam a impossibilidade do caso geral. Mas você sempre pode automatizar um caso específico. Não existe uma fórmula geral para a soma hipergeométrica, mas o algoritmo de Gosper prova que pelo menos um caso especial pode ser automatizado: o caso da soma hipergeométrica cujo resultado também é um termo hipergeométrico. E esse caso é imensamente útil na prática, especialmente se você estuda computação ou combinatória.

E no fim das contas, usar o Wolfram Alpha para resolver uma somatória é trapacear ou não? Por dentro, tudo que o Wolfram Alpha faz é usar o algoritmo de Gosper. Se você acredita que usar um método computacional é roubar, então por consistência precisa abrir de mão de todos os outros métodos computacionais. Não pode mais resolver a equação de segundo grau usando a fórmula de Bhaskara, não pode mais calcular o máximo divisor comum com o algoritmo de Euclides, não pode nem calcular quanto é 463 vezes 367 sem antes decorar a tabuada do 367.

Para mim isso não faz sentido. E por isso eu uso o Wolfram Alpha sem medo de ser feliz :)

O Paradoxo do Pokémon Twitch

A última moda da internet é o Pokémon Twitch. Como essas modas evaporam tão rapidamente quanto surgem, vale a pena explicar do que se trata:

O twitch é um site para livestreaming de videogames. Você conecta seu videogame no site, e as pessoas do mundo todo podem te assistir enquanto você joga. Ele tem também uma janela de chat, então as pessoas podem comentar enquanto assistem.

Mas algum gênio aprimorou a idéia original. Ele ligou o twitch em um emulador de Game Boy, e configurou o emulador para pegar o texto do chat e interpretar como se fosse o joystick. Com isso, qualquer um que tenha uma conta lá pode controlar o personagem do jogo (no caso, a versão original de Pokémon). E melhor ainda, como todos estão controlando o mesmo jogo, então o que ele fez foi criar uma espécie de crowdgaming, onde a sabedoria das massas escolhe o melhor caminho do personagem.

Quer dizer, na teoria. Na prática, a internet será a internet. Para cada um que tenta jogar sério, tem um trollzinho que fica mandando o comando oposto. No pico de popularidade, ele chegou a ter mais de cem mil pessoas jogando ao mesmo tempo. Imagine metade disso sacaneando ao invés de jogando, e dá pra ter uma idéia da loucura que é o Pokémon Twitch (ou então leia o FAQ para entender tudo que já aconteceu até agora).

Quando eu fiquei sabendo do Pokémon Twitch, eu fiz o que qualquer pessoa normal faria: criei um modelo matemático do jogo. No meu modelo simplificado, o personagem só anda na vertical, um passo por iteração. O número de jogadores e o número de trolls é o mesmo, então a cada passo ele tem 50% de chance de andar uma unidade para cima ou para baixo. Assuma que ele começa da origem. Nesse modelo, eu faço duas perguntas:

1. Depois de n iterações, onde você espera que o personagem esteja?
2. Depois de n iterações, qual você espera que seja a distância do jogador à origem?

Eu sei, eu sei: "Ricbit, deixa de ser bobo! As duas perguntas são iguais! Suponha que a resposta da primeira pergunta é y. Então a resposta da segunda é y menos 0, ou seja, o mesmo valor y".

Paradoxalmente, isso está errado! A resposta das duas perguntas é diferente! Probabilidade é um tópico que escapa facilmente da nossa intuição, então vale a pena fazer as contas com cuidado para entender a solução.

O valor esperado da posição

Para a resposta da primeira pergunta, vamos calcular qual é o valor esperado da posição. Na iteração 0 ele ainda não executou nenhum movimento, então sabemos a posição deterministicamente: ele está na origem.

E na iteração n+1? Depende de onde ele estava na iteração n. Tem 50% de chance de ter subido, e 50% de ter descido:

Com isso podemos calcular o valor esperado de y(n+1) direto da lei da expectativa total:

Ou seja, o valor esperado em qualquer iteração é sempre zero, na média nós esperamos que ele fique andando em círculos e nunca fique muito longe da origem.

Entretanto, note que só porque o valor esperado é zero não quer dizer que ele vai sempre terminar na origem: a média é zero, mas a variância não é. Um exercício bastante curioso é calcular qual a probabilidade exata do personagem estar na origem após n iterações. Isso é fácil e dá pra resolver com combinatória de colégio. Imagine que você tem n caixinhas:


Cada uma dessas caixinhas você pode preencher com +1 ou -1; se a soma de todas elas for zero, então o personagem termina na origem. De quantas maneiras podemos fazer isso? Note que nós só precisamos escolher as posições dos +1:

+1

+1

+1

+1

Sabendo onde estão os +1, a posição dos -1 restantes fica unicamente determinada. Precisamos então descobrir de quantas maneiras podemos encaixar n/2 números +1 em n caixinhas. Mas isso é a definição do binomial:

Agora é só dividir pelo total. Como cada caixinha tem duas opções possíveis, +1 e -1, então o número total de caixinhas é 2^n. Portanto, a probabilidade dele terminar na origem é:

Essa fórmula tem um problema: embora ela seja exata, é muito ruim de calcular quando o n é grande (tente calcular para n=1000, o número de dígitos da sua calculadora vai acabar rapidinho). Podemos achar um assintótico usando a aproximação do coeficiente binomial central:

Mais fácil né? Agora podemos calcular facilmente que, para n=1000, a chance de terminar na origem é aproximadamente 2.5%.

Se você tem o olho bom, deve ter percebido uma pegadinha na derivação da fórmula acima. Ela só funciona quando n é par! De fato, quando n é ímpar, não tem como o número de +1 e -1 serem iguais, e portanto a chance dele terminar na origem é zero. Olha que curioso: quando o n é ímpar, ele nunca termina em zero, apesar do valor esperado ser zero!

O valor esperado da distância

Vamos agora à segunda pergunta, que é o valor esperado da distância. Por que afinal esse valor é diferente do anterior?

Para exemplificar, vamos supor que n=2. Em uma das amostras, digamos que os jogadores ganharam e ele andou 2 unidades para cima, então a distância para a origem é 2. Em outra amostra, suponha que os trolls ganharam, então ele andou 2 unidades para baixo. Qual a distância da origem? É 2 também! Distâncias são sempre positivas!

Em outras palavras, o valor esperado da distância é o valor esperado do módulo da posição! E quanto é esse valor? As contas aqui são um pouco mais complicadas, então vou precisar da caixa azul:

Dessa vez nós vamos ter que calcular o valor esperado pela definição:

Quando k=0 o somando vale zero também. Vamos então supor k positivo e usar a mesma idéia das caixinhas. A soma era zero quando a quantidade de +1 e de -1 era igual. Dessa vez, nós queremos que a soma seja k, então precisamos ter k números +1 a mais que o número de -1. Portanto, a chance dele terminar em um valor k positivo é:

Quando k é negativo, a fórmula é exatamente a mesma! Afinal, é só trocar todos os +1 por -1 e vice versa. Portanto, o valor esperado da distância é:

Essa fórmula é difícil de manipular, para deixar mais fácil vamos supor que n é par. Se n é par, podemos escrevê-lo como n=2m. Mas olha só, se n for par, então você nunca vai terminar numa distância k que seja ímpar. Se você tirar k ímpar de n par, o que sobra é ímpar, então não tem como ter soma zero. Portanto podemos supor k par, e fazer k=2q. Substituindo:

A somatória é uma soma hipergeométrica, então você pode resolvê-la com o algoritmo de Gosper, chegando assim na forma fechada. Mas se você não souber usar o algoritmo de Gosper, não tem problema, é só provar a identidade abaixo por indução finita:

Substituindo a fórmula temos:

Agora é só usar a aproximação do binomial central:

Chegamos então no resultado final, o valor esperado da distância, que é raiz de 2n sobre pi.

A matemática experimental

Probabilidade tem uma vantagem grande sobre outros ramos da matemática: é muito fácil conferir as contas. Basta escrever uma simulação computacional com uma quantidade razoável de amostras!

Eu escrevi um script em python (o código está no github), e os resultados para n=1000 e cem mil amostras foram os seguintes:

• Posição: calculado 0.00, medido 0.05
• Distância: calculado 25.23, medido 25.27
• Término na origem: calculado 2.52%, medido 2.53%

Ou seja, a matemática experimental concordou com a matemática teórica.

A lição que fica dessa brincadeira é tomar muito cuidado com suas intuições sobre probabilidade, porque ela tem uma tendência a nos enganar bem facilmente :)