Aplicativos

Você conhece os valores de uma função \(f\) em \(n+1\) pontos distintos \(x_0\lt x_1\lt \dots\lt x_n\), os nós, digamos \(y_j=f(x_j)\). Os valores podem vir de uma medição, ou de uma fórmula cara de avaliar. A interpolação polinomial procura um polinômio \(p\) de grau no máximo \(n\) que assuma os mesmos valores nos nós:

Por que polinômios? São as funções mais simples que existem: calcular um exige apenas somas e produtos, e pode-se derivá-lo e integrá-lo exatamente. O teorema de Weierstrass diz que toda função contínua num intervalo fechado é limite uniforme de polinômios. E a interpolação é o motor de muitos outros métodos: a regra de Simpson e todas as fórmulas de quadratura, diferenças finitas, os métodos de Adams para equações diferenciais, o método da secante e a interpolação inversa para raízes, a extrapolação de Richardson e a integração de Romberg.

As perguntas deste curso são as naturais. Esse \(p\) existe, e é único? Como calculá-lo, e qual o melhor modo? Quão grande é \(f-p\) entre os nós? Dá para escolher bem os nós? O que acontece quando \(n\) cresce?

Teorema. Se os nós \(x_0,\dots,x_n\) são distintos, então para quaisquer valores \(y_0,\dots,y_n\) existe exatamente um polinômio \(p\) de grau no máximo \(n\) com \(p(x_j)=y_j\) para todo \(j\).

Unicidade. Se \(p\) e \(q\) interpolam ambos, então \(p-q\) tem grau no máximo \(n\) e se anula nos \(n+1\) nós distintos. Um polinômio não nulo de grau no máximo \(n\) tem no máximo \(n\) raízes, logo \(p-q=0\).

Existência. Escreva \(p(x)=a_0+a_1x+\dots+a_nx^n\). As condições formam um sistema linear \(Va=y\) com a matriz de Vandermonde

A unicidade diz que \(Va=0\) só para \(a=0\), logo \(V\) é invertível e o sistema tem solução. (O próximo grupo constrói a solução explicitamente, o que é uma segunda demonstração.)

Consequências. Se os dados vêm de um polinômio de grau no máximo \(n\), o interpolador é esse mesmo polinômio: escolha “polinômio x⁵ − 2x³ + x/2” no aplicativo e aumente \(n\) de 4 para 5; o erro cai de 0,11 para cerca de \(10^{-15}\). O grau de \(p\) pode ser menor que \(n\) (todos os \(y_j\) iguais dão uma constante). Os nós têm de ser distintos; com nós repetidos é preciso prescrever derivadas (interpolação de Hermite, grupo 6).

Fixe os nós e chame \(\mathcal I_nf\) o polinômio interpolador de \(f\). Três fatos organizam tudo o que vem depois.

• \(\mathcal I_n\) é linear: \(\mathcal I_n(\alpha f+\beta g)=\alpha\,\mathcal I_nf+\beta\,\mathcal I_ng\), porque o sistema \(Va=y\) é linear em \(y\).
• \(\mathcal I_n\) é uma projeção: \(\mathcal I_np=p\) para todo polinômio \(p\) de grau no máximo \(n\), logo \(\mathcal I_n^2=\mathcal I_n\).
• \(\mathcal I_nf\) depende de \(f\) apenas pelos \(n+1\) números \(f(x_j)\), então nada sabe sobre \(f\) entre os nós. O tamanho de \(f-\mathcal I_nf\) ali tem de vir da suavidade (grupo 3), e a sensibilidade a erros nos \(y_j\) é a norma de \(\mathcal I_n\), a constante de Lebesgue (grupo 4).

O polinômio em si não depende de como o escrevemos. Os cinco algoritmos do próximo grupo (Lagrange, Newton, Neville, monomial, baricêntrico) produzem a mesma função, com custos e comportamentos de arredondamento muito diferentes.

Fixe os nós \(x_0,\dots,x_n\) (distintos). Para \(j=0,\dots,n\) chame de \(e_j\) a tabela que vale \(1\) em \(x_j\) e \(0\) nos outros \(n\) nós: \(e_j(x_k)=\delta_{jk}\). Cada uma dessas tabelas tem um polinômio interpolador, e ele é único; chame-o de \(\ell_j\). Toda tabela \(y\) satisfaz \(y=\sum_jy_je_j\) (compare os dois lados nó por nó), e a aplicação tabela \(\mapsto\) polinômio interpolador é linear, logo

O ponto que nunca é demais sublinhar é o primeiro: os polinômios de Lagrange não são só uma fórmula que por acaso funciona, são os polinômios interpoladores das tabelas mais simples possíveis, construídas sobre os mesmos nós. O problema geral foi reduzido a estes \(n+1\) problemas especiais, e os problemas especiais têm resposta explícita.

Procure um polinômio \(q\) de grau no máximo \(n\) com \(q(x_k)=0\) para \(k\ne j\) e \(q(x_j)=1\). Tendo as \(n\) raízes \(x_k\) (\(k\ne j\)) e grau no máximo \(n\), ele tem de ser \(q(x)=c\prod_{k\ne j}(x-x_k)\) para alguma constante \(c\): dividindo por esses \(n\) fatores lineares sobra um polinômio de grau no máximo \(0\). A condição \(q(x_j)=1\) dá \(c=1/\prod_{k\ne j}(x_j-x_k)\), isto é

Então a fórmula não é adivinhada, é deduzida da tabela especial. O mesmo argumento de contagem de raízes prova a unicidade em geral (dois interpoladores diferem por um polinômio de grau \(\le n\) com \(n+1\) raízes), então é o único polinômio que serve. Experimente no cartão acima: escolha a tabela especial \(e_j\) e compare as hastes (os valores \(0,\dots,1,\dots,0\)) com a curva; a nota abaixo do gráfico diz o maior valor de \(|\ell_j|\), que pode ser bem maior que \(1\) entre os nós.

Com \(n\le6\) o cartão mostra a tabela das peças \(y_je_j\): a coluna \(j\) é nula exceto na linha \(j\), onde vale \(y_j\), e a última coluna, a soma de cada linha, é a tabela original \(y\). À direita do cartão as curvas \(y_j\ell_j\) são os polinômios interpoladores dessas colunas, e o controle soma uma a uma: as somas parciais são os polinômios interpoladores das tabelas parciais \(y_0e_0+\dots+y_me_m\), e depois da última peça a curva passa por todos os nós. No modo “seus dados”, arraste um ponto para cima e veja exatamente uma peça mudar de escala, e a soma acompanhar.

Três identidades são a mesma ideia lida com outras tabelas. A tabela de todos uns é \(e_0+\dots+e_n\) e é interpolada pela constante \(1\), logo \(\sum_j\ell_j=1\). A tabela \((x_0^k,\dots,x_n^k)\) é interpolada por \(x^k\) quando \(k\le n\), logo \(\sum_jx_j^k\ell_j(x)=x^k\). E como \(\ell_j\) tem grau \(n\) e as \(n+1\) funções \(\ell_0,\dots,\ell_n\) são independentes (olhe seus valores nos nós), elas formam uma base dos polinômios de grau no máximo \(n\), em que o vetor de coordenadas de \(p\) é simplesmente a tabela \((y_0,\dots,y_n)\).

Escolha uma base \(b_0,\dots,b_n\) dos polinômios de grau no máximo \(n\) e escreva \(p=\sum_ka_kb_k\). As condições de interpolação viram o sistema linear \(Ba=y\) com \(B_{ik}=b_k(x_i)\): a matriz da base avaliada nos nós. As três formas clássicas são três escolhas:

Quanto mais estrutura a base tem em relação aos nós, menos trabalho sobra. A base de Lagrange é o caso extremo. Ela também diz o que é \(V^{-1}\): como \(\ell_j(x)=\sum_ka_{kj}x^k\) é o polinômio da tabela \(e_j\), a \(j\)-ésima coluna de \(V^{-1}\) guarda os coeficientes de \(\ell_j\) nos monômios, e o operador de interpolação, como aplicação de tabelas em vetores de coeficientes, é a multiplicação por \(V^{-1}\). Seu tamanho na norma do sup é a constante de Lebesgue do grupo 4: \(|p(x)|\le\Lambda_n\max_j|y_j|\), pois \(|p(x)|\le\sum_j|y_j||\ell_j(x)|\le(\max_j|y_j|)\sum_j|\ell_j(x)|\).

O nome sugere que Lagrange inventou a fórmula, mas ela é mais antiga. Edward Waring a publicou em 1779, e Euler deduziu a mesma fórmula poucos anos depois (por volta de 1783). A apresentação do próprio Lagrange veio nas suas aulas na École Normale por volta de 1795. O que o nome homenageia é uma maneira de ver: com Lagrange a fórmula aparece dentro de uma teoria geral dos \(\ell_j\) como tijolos de construção, o ponto de vista deste cartão, e foi nessa forma que ela se difundiu. Uma geração depois Gauss, numa conferência de 1812 publicada só depois, mostrou que as fórmulas de Newton, Euler e Waring–Lagrange são um único polinômio escrito de maneiras diferentes, que é o teorema de unicidade do grupo 1 disfarçado. A forma do próprio Newton, mais antiga que todas essas, é o tema do próximo aplicativo desta seção.

Fontes: J. Meijering, A chronology of interpolation: from ancient astronomy to modern signal and image processing, Proc. IEEE 90 (2002) 319–342, para as datas; o ano exato das aulas de Lagrange é o dado ali, até onde pude conferir.

Em vez de resolver um sistema, construa \(p\) diretamente. Para cada \(j\) seja

Cada \(\ell_j\) é um produto de \(n\) fatores lineares, logo tem grau \(n\). Num nó \(x_i\ne x_j\) o numerador contém o fator \(x_i-x_i=0\); em \(x_j\) cada fração vale 1. Portanto \(\ell_j(x_i)=\delta_{ij}\) (delta de Kronecker), e \(p(x_i)=\sum_jy_j\delta_{ij}=y_i\). Esta é a demonstração explícita de existência.

A leitura conceitual, de que \(\ell_j\) é o interpolador da tabela especial \(e_j\) e de que a linearidade monta todo o resto, é o tema do grupo 1b.

Exemplo. Nós \(-1,0,1\): \(\ell_0=\tfrac12x(x-1)\), \(\ell_1=1-x^2\), \(\ell_2=\tfrac12x(x+1)\). Para os valores \(y=(1,0,1)\) obtemos \(p=\ell_0+\ell_2=x^2\).

Duas identidades. A constante 1 interpola os valores \((1,\dots,1)\), então por unicidade \(\sum_j\ell_j(x)=1\) para todo \(x\) (a linha tracejada no cartão da base). Mais geralmente \(\sum_jx_j^k\ell_j(x)=x^k\) para \(k\le n\).

Avaliação. A forma de Lagrange é maravilhosa para a teoria, pois a base \(\ell_j\) está ligada aos nós e os coeficientes são os próprios dados. Como algoritmo é fraca: calcular \(p(x)\) custa \(O(n^2)\), e acrescentar um nó muda todas as \(\ell_j\).

Newton escreve \(p\) na base \(1,\ (x-x_0),\ (x-x_0)(x-x_1),\dots\):

Os coeficientes são as diferenças divididas, definidas por

Elas preenchem uma tabela triangular (veja o cartão Formas; os coeficientes são a linha de cima, destacada). Exemplo. Nós \(0,1,2,3\), valores \(1,2,0,5\): as colunas são \(1,2,0,5\); depois \(1,-2,5\); depois \(-\tfrac32,\tfrac72\); depois \(\tfrac53\). Logo

e de fato \(p(0)=1,\ p(1)=2,\ p(2)=0,\ p(3)=5\) e \(p(1{,}5)=0{,}75\).

Por que é útil. O coeficiente \(c_k\) usa só os primeiros \(k+1\) nós, então acrescentar um nó acrescenta um termo novo e deixa os antigos em paz (Lagrange não consegue isso). A avaliação é aninhada como a regra de Horner, \(p(x)=c_0+(x-x_0)\bigl(c_1+(x-x_1)(c_2+\cdots)\bigr)\), com \(n\) multiplicações. A tabela custa cerca de \(n^2/2\) divisões. As diferenças divididas são simétricas nos argumentos, e para \(f\) suave vale \(f[x_0,\dots,x_k]=f^{(k)}(\xi)/k!\) para algum \(\xi\) (consequência da fórmula do erro do grupo 3); quando todos os nós se fundem a forma de Newton vira o polinômio de Taylor.

Se você precisa de \(p\) em um único ponto \(t\), pode dispensar os coeficientes. Seja \(P_{i,\dots,i+k}(t)\) o valor em \(t\) do polinômio que interpola somente os nós \(x_i,\dots,x_{i+k}\). Partindo de \(P_i=y_i\),

Demonstração. Os dois lados são polinômios em \(t\) de grau no máximo \(k\). Em \(t=x_i\) o lado direito dá \(P_{i,\dots,i+k-1}(x_i)=y_i\); em \(t=x_{i+k}\) dá \(P_{i+1,\dots,i+k}(x_{i+k})=y_{i+k}\); num nó interno as duas parcelas valem \(y\) e os pesos \((t-x_i)-(t-x_{i+k})=x_{i+k}-x_i\) somam o denominador. Por unicidade os dois lados coincidem.

O custo é \(O(n^2)\) por ponto e só um vetor de memória. As colunas sucessivas são os interpoladores de grau crescente, exatamente o que a extrapolação de Richardson e a integração de Romberg usam: extrapolar para \(h=0\) é calcular um polinômio interpolador em 0.

O algoritmo mais óbvio resolve \(Va=y\) por eliminação gaussiana. Custa \(O(n^3)\) e, pior, \(V\) é mal condicionada. O número de condição \(\kappa_\infty(V)\) (nós em \([-1,1]\)) cresce geometricamente:

Com 16 dígitos de aritmética, um número de condição \(\kappa\) custa cerca de \(\log_{10}\kappa\) dígitos nos coeficientes. A visão “Baricêntrica × Vandermonde” do cartão Erro mostra o efeito: dois cálculos do mesmo polinômio, cuja diferença deveria ser zero, diferem por \(3{,}6\cdot10^{-14}\) (n = 10), \(2{,}9\cdot10^{-11}\) (n = 20) e \(3{,}1\cdot10^{-8}\) (n = 30) para a função de Runge em nós de Chebyshev. Para \(e^x\) em nós equiespaçados a diferença é \(8\cdot10^{-13}\) em n = 20, \(5\cdot10^{-10}\) em n = 30, \(3\cdot10^{-7}\) em n = 40 e \(10^{-4}\) em n = 50.

(O aplicativo já ajuda o caminho de Vandermonde: reescala o intervalo para \([-1,1]\) e faz pivoteamento. Num intervalo cru como \([0,100]\) é catastroficamente pior.) A moral: a base monomial é perfeita para demonstrações e para \(n\) minúsculo, e uma base ruim para calcular.

Ponha \(\omega(x)=\prod_k(x-x_k)\) e \(w_j=1/\prod_{k\ne j}(x_j-x_k)\) (que é igual a \(1/\omega'(x_j)\)). Então \(\ell_j(x)=\omega(x)\,\dfrac{w_j}{x-x_j}\), e a forma de Lagrange vira a primeira forma baricêntrica

Interpolar a constante 1 dá \(1=\omega(x)\sum_jw_j/(x-x_j)\). Dividindo a primeira forma por essa identidade eliminamos \(\omega\) e obtemos a segunda forma baricêntrica (a verdadeira)

Em \(x=x_j\) a fórmula é \(0/0\) e simplesmente devolvemos \(y_j\). Como os pesos aparecem no numerador e no denominador, qualquer fator comum é irrelevante, e por isso existem fórmulas fechadas simples:

• nós equiespaçados: \(w_j=(-1)^j\binom nj\);
• pontos de Chebyshev (zeros de \(T_{n+1}\)): \(w_j=(-1)^j\sin\frac{(2j+1)\pi}{2n+2}\);
• pontos de Chebyshev–Lobatto (extremos de \(T_n\), com as pontas): \(w_j=(-1)^j\), pela metade em \(j=0\) e \(j=n\).

Para outros nós os pesos custam \(O(n^2)\) uma única vez. O aplicativo usa as fórmulas fechadas e as confere com a fórmula geral (concordância até \(5\cdot10^{-15}\)).

Propriedades. Avaliar custa \(O(n)\) por ponto; acrescentar um nó custa \(O(n)\); as fórmulas servem para quaisquer nós. Higham (2004) provou que a segunda forma é estável (progressivamente) quando a constante de Lebesgue dos nós é pequena, como nos nós de Chebyshev, o que a torna o método de escolha (é o que o Chebfun usa). Com nós equiespaçados e \(n\) grande o próprio problema de interpolação é mal condicionado (grupo 4), e nenhum algoritmo conserta isso.

O cartão Formas calcula o mesmo \(p(x)\) com cinco algoritmos. Para \(n\) pequeno os cinco valores concordam em 15 ou 16 dígitos; aumente \(n\) com nós equiespaçados e o valor de Vandermonde é o primeiro a errar.

Teorema. Seja \(f\in C^{n+1}[a,b]\), sejam \(x_0,\dots,x_n\in[a,b]\) distintos e seja \(p\) o polinômio interpolador. Para todo \(x\in[a,b]\) existe um ponto \(\xi\) no menor intervalo que contém \(x,x_0,\dots,x_n\) tal que

Demonstração. Se \(x\) é um nó os dois lados são zero. Senão fixe \(x\), e defina \(c=\dfrac{f(x)-p(x)}{\omega(x)}\) e \(g(t)=f(t)-p(t)-c\,\omega(t)\). Então \(g\) se anula nos \(n+1\) nós (ali \(f-p=0\) e \(\omega=0\)) e em \(t=x\) pela escolha de \(c\): são \(n+2\) zeros distintos. Pelo teorema de Rolle \(g'\) tem pelo menos \(n+1\) zeros, \(g''\) pelo menos \(n\), e assim por diante até que \(g^{(n+1)}\) tem pelo menos um zero \(\xi\). Mas \(p^{(n+1)}=0\) e \(\omega^{(n+1)}=(n+1)!\), então

Corolários. \(\|f-p\|_\infty\le\dfrac{\|f^{(n+1)}\|_\infty}{(n+1)!}\,\|\omega\|_\infty\). Se \(f\) é um polinômio de grau no máximo \(n\), o erro é zero. Na interpolação linear (\(n=1\), espaçamento \(h\)) o máximo de \(|\omega|\) é \(h^2/4\), logo o erro é no máximo \(h^2\|f''\|/8\). Há também uma versão que não precisa de derivadas: \(f(x)-p(x)=\omega(x)\,f[x_0,\dots,x_n,x]\), que é a forma de Newton com mais um nó.

O erro é um produto de duas coisas com donos diferentes. O fator \(f^{(n+1)}(\xi)/(n+1)!\) pertence à função; não podemos mudá-lo. O fator \(\omega(x)\) pertence aos nós; nós o escolhemos. Use a visão “Polinômio dos nós” do cartão Erro.

Em \([-1,1]\) o maior \(|\omega|\) é

Em \(n=20\) o \(\omega\) de Chebyshev é cerca de 245 vezes menor. Para nós equiespaçados \(\omega\) é minúsculo no meio e grande perto das duas pontas, então o erro se acumula ali; para nós de Chebyshev os \(n+2\) lóbulos de \(|\omega|\) têm todos a mesma altura.

O fator da derivada também pode ser perigoso. Para \(1/(1+25x^2)\) os coeficientes de Taylor em 0 são \(\pm25^k\), então \(f^{(n+1)}(0)/(n+1)!\) tem tamanho em torno de \(5^{n+1}\). O fator \(\omega\) dos nós equiespaçados encolhe só cerca de 0,7 por passo (8,5·10⁻³ em \(n=10\), 2,3·10⁻⁴ em \(n=20\)), então a cota cresce com \(n\), e o erro real também (grupo 3, Runge).

Tome \(f=e^x\) em \([-1,1]\), onde \(\|f^{(n+1)}\|_\infty=e\). Com nós de Chebyshev \(\|\omega\|=2^{-n}\), então

A cota vale e é justa dentro de um fator de cerca de 2,4 para os nós de Chebyshev. Use-a para escolher \(n\): para garantir erro abaixo de \(10^{-6}\) precisamos de \(e/((n+1)!2^n)\lt 10^{-6}\), o que falha para \(n=6\) (\(8{,}4\cdot10^{-6}\)) e vale para \(n=7\) (\(5{,}3\cdot10^{-7}\)). Como \(e^x\) é inteira, o erro cai mais depressa que qualquer taxa geométrica até chegar ao nível do arredondamento, cerca de \(10^{-15}\) a partir de \(n=15\).

Escolha “Runge” com nós equiespaçados. Em \(n=5\) o polinômio é uma aproximação pobre mas sensata; em \(n=10\) ele balança bem acima da função perto das pontas; em \(n=20\) sai do gráfico. Os erros máximos em \([-1,1]\) para \(1/(1+25x^2)\):

Em \(x=0{,}95\) e \(n=10\) o polinômio vale 1,92 enquanto \(f(0{,}95)=0{,}042\). Mais pontos o deixaram pior. Runge descobriu isso em 1901: os interpoladores convergem para \(f\) só para \(|x|\lt 0{,}726\) e divergem para \(|x|>0{,}726\) (aproximadamente).

A razão está no plano complexo: \(1/(1+25z^2)\) tem polos em \(\pm i/5\), perto do intervalo. A interpolação equiespaçada só converge para funções analíticas numa região em forma de lente em volta de \([-1,1]\) que corta o eixo imaginário em cerca de \(\pm0{,}53\,i\); os polos em \(\pm0{,}2\,i\) estão dentro dela. A mesma história vale para “1/(1 + x²) em [−5, 5]”: é a mesma função após reescalar, e os números são idênticos, o que mostra que só importa a posição relativa da singularidade. Nós de Chebyshev curam isso (próxima lição e grupo 5).

Como \(\omega\) é a parte que controlamos, escolha os nós para tornar \(\|\omega\|_\infty\) o menor possível. Teorema (Chebyshev). Entre todos os polinômios mônicos de grau \(n+1\), a menor norma do máximo em \([-1,1]\) é a de \(T_{n+1}(x)/2^n\), onde \(T_{n+1}(x)=\cos((n+1)\arccos x)\); sua norma é \(2^{-n}\). Logo os melhores nós são os zeros de \(T_{n+1}\),

Demonstração. \(T_{n+1}\) tem coeficiente líder \(2^n\), então \(\omega^\ast=T_{n+1}/2^n\) é mônico, e alterna entre \(+2^{-n}\) e \(-2^{-n}\) nos \(n+2\) pontos \(\cos(k\pi/(n+1))\). Se um mônico \(q\) tivesse \(\|q\|_\infty\lt 2^{-n}\), então \(r=\omega^\ast-q\) teria grau no máximo \(n\) e, nesses \(n+2\) pontos, os sinais de \(\omega^\ast\): \(n+1\) mudanças de sinal, logo \(n+1\) zeros, logo \(r=0\), contradizendo \(\|q\|\lt 2^{-n}\). \(\blacksquare\)

Em \([a,b]\) os nós são levados por uma transformação afim e a cota da lição anterior vira

Os pontos de Chebyshev–Lobatto (extremos de \(T_n\), com as pontas) dão um \(\omega\) de tamanho quase igual e são muitas vezes mais convenientes. Medidos na função de Runge, os erros máximos para nós de Chebyshev em \(n=10,20,40\) são \(0{,}109,\ 0{,}0153,\ 2{,}9\cdot10^{-4}\); para Lobatto \(0{,}132,\ 0{,}0177,\ 3{,}4\cdot10^{-4}\). Uma curiosidade: o erro de Chebyshev em \(n=15\) (0,083) é maior que em \(n=12\) (0,069). Para \(n\) ímpar há um número par de nós e nenhum está em \(x=0\), justamente onde fica o pico da função.

Agora perturbe os dados. Se cada \(y_j\) muda no máximo \(\delta\), o interpolador muda \(\sum_j\delta_j\ell_j(x)\), cujo tamanho é no máximo \(\delta\lambda_n(x)\), onde

A cota é atingida tomando \(\delta_j=\delta\,\mathrm{sinal}\,\ell_j(x)\). Então \(\Lambda_n\), a constante de Lebesgue, é exatamente a amplificação de erros dos dados no pior caso: o número de condição da interpolação. Depende só dos nós. Veja a função de Lebesgue no cartão da base. Medidas em \([-1,1]\):

No aplicativo, tome \(\sin(2\pi x)\), \(n=20\), e some ruído \(10^{-3}\) aos valores: com nós equiespaçados o erro máximo vira cerca de 1,8, com nós de Chebyshev cerca de \(1{,}2\cdot10^{-3}\), o tamanho do próprio ruído.

Seja \(E_n(f)=\min_{q\in\mathcal P_n}\|f-q\|_\infty\) o erro do melhor polinômio de grau no máximo \(n\) (existe e é único, pelo teorema da equioscilação de Chebyshev). A interpolação nunca é muito pior:

Demonstração. Seja \(q\) um polinômio qualquer de grau no máximo \(n\). Como \(\mathcal I_nq=q\) e \(\mathcal I_n\) é linear, \(f-\mathcal I_nf=(f-q)-\mathcal I_n(f-q)\), logo \(\|f-\mathcal I_nf\|\le(1+\|\mathcal I_n\|)\|f-q\|\). Tome \(q\) o melhor polinômio. \(\blacksquare\)

Assim duas coisas controlam a interpolação: quão bem polinômios aproximam \(f\) (\(E_n(f)\), questão de suavidade, próximo grupo) e quão bons são os nós (\(\Lambda_n\)). Com nós de Chebyshev \(1+\Lambda_n\) é 3,9 em \(n=40\) e cerca de 5,4 mesmo em \(n=400\): você perde no máximo um fator pequeno em relação ao melhor polinômio, e o obtém com uma fórmula \(O(n)\) e sem iteração. Com nós equiespaçados o fator \(1+\Lambda_n\) é astronômico e a garantia não vale nada.

As assintóticas, todas visíveis no gráfico “Constante de Lebesgue contra n”:

• nós equiespaçados: \(\Lambda_n\approx\dfrac{2^{n+1}}{e\,n\ln n}\) (cerca de \(1{,}1\cdot10^{4}\) medido em \(n=20\), \(4{,}7\cdot10^{9}\) em \(n=40\)): crescimento exponencial, praticamente dobrando a cada passo;
• nós de Chebyshev–Lobatto: \(\Lambda_n\approx\dfrac2\pi\ln n+0{,}9625\), conferido numericamente com quatro dígitos em \(n=100\) (3,894) e \(n=400\) (4,777); os zeros de Chebyshev são o mesmo a menos de uma constante: crescimento logarítmico.

Existe um conjunto perfeito de nós com \(\Lambda_n\) limitado? Não. Erdős (1961) provou que para toda escolha de nós \(\Lambda_n\ge\frac2\pi\ln(n+1)-C\) para uma constante \(C\). Logo os nós de Chebyshev são ótimos a menos de uma constante aditiva limitada, o melhor que se pode esperar.

Olhe a função de Lebesgue \(\lambda_n(x)\) para nós equiespaçados em \(n=20\): é de ordem 1 a 10 no meio do intervalo e dispara para cerca de \(10^4\) no primeiro e no último subintervalo. Nós afastados demais perto das pontas é o que faz um polinômio por eles balançar loucamente ali; os nós de Chebyshev se aglomeram perto das pontas com densidade \(1/(\pi\sqrt{1-x^2})\) para impedir exatamente isso.

Mesmo que \(f\) seja calculada perfeitamente, cada \(y_j\) é armazenado com erro relativo de até cerca de \(\varepsilon\approx1{,}1\cdot10^{-16}\). Isso é uma perturbação dos dados, e a interpolação pode multiplicá-la por \(\Lambda_n\). Duas medições:

• \(f=e^x\), equiespaçados, \(n=40\): \(\Lambda_n=4{,}7\cdot10^9\), dano esperado \(\sim5\cdot10^{-7}\), erro medido \(3{,}3\cdot10^{-7}\), embora o erro teórico de interpolação seja muito abaixo de \(10^{-30}\) (a função é inteira).
• O polinômio \(x^5-2x^3+x/2\) é reproduzido exatamente para todo \(n\ge5\) em aritmética exata. Em ponto flutuante, nós equiespaçados dão erros \(1{,}1\cdot10^{-15}\) (n = 10), \(5{,}5\cdot10^{-13}\) (n = 20), \(2{,}4\cdot10^{-10}\) (n = 30) e \(1{,}6\cdot10^{-7}\) (n = 40), enquanto nós de Chebyshev ficam abaixo de \(10^{-15}\) o tempo todo.

Isso é uma propriedade do problema, não da fórmula: atinge baricêntrica, Newton e Lagrange por igual. Para \(n\) grande, a interpolação equiespaçada está portanto duplamente errada: imprecisa para funções como a de Runge e ruidosa até para polinômios.

O palpite tentador é que \(\mathcal I_nf\to f\) uniformemente para toda \(f\) contínua se os nós são decentes. É falso. Teorema de Faber (1914): para toda matriz triangular de nós existe uma função contínua \(f\) para a qual \(\|f-\mathcal I_nf\|_\infty\) não tende a 0. A razão é o princípio da limitação uniforme: \(\Lambda_n=\|\mathcal I_n\|\to\infty\) para toda matriz (a cota de Erdős do grupo 4), então alguma \(f\) é atingida pelo crescimento.

Para nós equiespaçados é pior que “alguma função”: Bernstein (1918) mostrou que para \(f(x)=|x|\) os interpoladores divergem em todo ponto de \([-1,1]\) exceto \(-1,0,1\). Tente “|x|” com nós equiespaçados: os erros máximos são 0,66, 95 e \(1{,}5\cdot10^7\) para \(n=10,20,40\). Com nós de Chebyshev são 0,055, 0,029 e 0,015: convergência, devagar.

O lado positivo: com nós de Chebyshev \(\Lambda_n=O(\ln n)\), então \(\|f-\mathcal I_nf\|\le(1+\Lambda_n)E_n(f)\) tende a 0 sempre que \(E_n(f)\ln n\to0\), em particular para toda função lipschitziana, e toda mais suave. (Marcinkiewicz mostrou que para cada \(f\) contínua alguma matriz de nós dá convergência, mas a matriz depende de \(f\), o que não serve na prática.)

Pela propriedade de quase-melhor, a velocidade é a da melhor aproximação \(E_n(f)\), a menos de um fator \(1+\Lambda_n\). E \(E_n(f)\) é governada pela suavidade:

• (Jackson) se \(f\in C^k[-1,1]\) então \(E_n(f)\le C_k\,\|f^{(k)}\|_\infty/n^k\): lei de potência, tanto mais rápida quanto mais suave;
• se \(f\) é analítica numa vizinhança de \([-1,1]\), então \(E_n(f)\) decai geometricamente, \(O(\rho^{-n})\) (próxima lição);
• uma quina custa uma potência: para \(f=|x|\), \(nE_n(f)\to0{,}2802\) (constante de Bernstein). A interpolação de Chebyshev dá \(0{,}0546,\ 0{,}0285,\ 0{,}0146\) em \(n=10,20,40\), cerca do dobro de \(0{,}28/n\), como o fator \(1+\Lambda_n\approx3\) permite.

A conclusão é prática. A precisão que se alcança com um polinômio de grau \(n\) é uma propriedade da função, e a interpolação de Chebyshev chega a um fator pequeno dela. Nada mais esperto bate \(E_n\); só pode evitar o fator \(1+\Lambda_n\).

Seja \(E_\rho\) a elipse com focos \(\pm1\) e semieixos \((\rho+\rho^{-1})/2\) e \((\rho-\rho^{-1})/2\). Teorema (Bernstein; veja Trefethen). Se \(f\) é analítica dentro de \(E_\rho\) com \(|f|\le M\) ali, então o interpolador em \(n+1\) pontos de Chebyshev satisfaz

Convergência geométrica, com taxa fixada pela singularidade mais próxima. Teste com a função de Runge. Seus polos em \(\pm i/5\) estão em \(E_\rho\) quando \((\rho-\rho^{-1})/2=1/5\), isto é \(\rho=\frac{1}{5}+\sqrt{1+\frac1{25}}=1{,}2198\). A previsão: o erro encolhe por \(\rho^{-10}=0{,}137\) a cada 10 passos de \(n\) e por \(\rho^{-20}=0{,}0188\) a cada 20. Medido com nós de Chebyshev: de \(n=10\) a 20 o fator é 0,140, de \(n=20\) a 40 é 0,0189. A teoria prevê a inclinação.

A taxa aparece como uma reta no gráfico logarítmico de convergência. Para \(\sin(2\pi x)\) e \(e^x\), funções inteiras sem singularidade alguma, a curva se curva para baixo mais depressa que qualquer reta (convergência superexponencial) até encontrar o piso do arredondamento perto de \(10^{-15}\). A mesma teoria explica a falha dos equiespaçados: esses nós se comportam como se as singularidades perto do eixo imaginário (os polos da função de Runge) estivessem muito mais perto do que estão.

Erro máximo com nós de Chebyshev, medido no aplicativo (escolha a função e ajuste n):

Lendo as linhas: para o salto, nada ajuda; para a quina e a singularidade na ponta o erro cai pela metade quando \(n\) dobra (\(O(1/n)\)); para o bump cai mais depressa que qualquer potência de \(1/n\) mas mais devagar que qualquer taxa geométrica; para Runge cai geometricamente à taxa \(1{,}22^{-n}\); para funções inteiras cai superexponencialmente até o arredondamento. O gráfico “Erro máximo contra n” mostra exatamente essas formas.

Ponha \(x=\cos\theta\). Os pontos de Chebyshev são \(x_j=\cos\theta_j\) com ângulos igualmente espaçados \(\theta_j=(2j+1)\pi/(2n+2)\): são as sombras, no diâmetro, de pontos igualmente espaçados no círculo unitário. Um polinômio de grau \(n\) em \(x\) é o mesmo que um polinômio trigonométrico par \(\sum_{k\le n}a_k\cos k\theta\), pois \(T_k(\cos\theta)=\cos k\theta\).

Então interpolar em pontos de Chebyshev é interpolar um polinômio trigonométrico em ângulos equiespaçados, e isso é a transformada discreta de Fourier: ortogonal, perfeitamente condicionada, e a fonte do crescimento \(\frac2\pi\ln n\) de \(\Lambda_n\). Pontos equiespaçados em \(x\) não são equiespaçados em \(\theta\), e essa é toda a diferença entre Chebyshev e Runge.

Bônus prático: os coeficientes de Chebyshev \(a_k\) do interpolador podem ser calculados em \(O(n\log n)\) pela transformada rápida de Fourier, que é a base dos métodos espectrais e do sistema Chebfun.

Troque \(f\) por \(p\) e integre: \(\displaystyle\int_a^bf\approx\int_a^bp=\sum_jW_jy_j\) com \(W_j=\int_a^b\ell_j\). É assim que nasce toda regra de quadratura interpolatória. Para nós equiespaçados em \([0,n]\) (espaçamento 1):

Para \(n=8\) os pesos são \(0{,}279,\ 1{,}662,\ -0{,}262,\ 2{,}962,\ -1{,}281,\ 2{,}962,\ -0{,}262,\ 1{,}662,\ 0{,}279\): alguns são negativos, então \(\sum|W_j|>b-a\) e erros de arredondamento nos \(y_j\) são amplificados, a mesma doença da constante de Lebesgue. Por isso as regras de Newton–Cotes de ordem alta são abandonadas. Remédios: regras compostas (grau baixo em muitos painéis pequenos), Clenshaw–Curtis (integrar o interpolador de Chebyshev) e a quadratura de Gauss, em que os nós são escolhidos para tornar a regra exata até o grau \(2n+1\).

Derivar a fórmula do erro dá regras para \(f'\): de \(n=2\) sai a diferença central \((f(x+h)-f(x-h))/(2h)\) com erro \(-h^2f'''(\xi)/6\). Derivar é intrinsecamente mal posto: a norma da matriz de derivação de Chebyshev cresce como \(n^2\), então o ruído é amplificado muito mais que na interpolação. Os métodos de passo múltiplo de Adams–Bashforth para EDOs integram o interpolador de \(f\) pelos últimos passos.

Prescreva valores e primeiras derivadas em \(n+1\) nós: \(2n+2\) condições, e exatamente um polinômio \(H\) de grau no máximo \(2n+1\) as satisfaz. A demonstração da fórmula do erro vale palavra por palavra (cada nó é zero duplo de \(f-H\)) e dá

A interpolação de Hermite se calcula com a tabela de Newton repetindo nós e pondo \(f[x_j,x_j]=f'(x_j)\) nas posições onde haveria divisão por zero. Para \(n=1\) é o interpolador cúbico de Hermite, o tijolo dos splines em computação gráfica e da interpolação monótona “PCHIP”.

O caso extremo é um único nó repetido \(n+1\) vezes: a forma de Newton vira o polinômio de Taylor \(\sum_kf^{(k)}(x_0)(x-x_0)^k/k!\), e a fórmula do erro vira o resto de Lagrange \(f^{(n+1)}(\xi)(x-x_0)^{n+1}/(n+1)!\). Então o teorema de Taylor é o caso especial desta lição em que todos os nós colidem, e a interpolação é sua generalização natural a nós separados, que é por que se comporta muito melhor num intervalo inteiro.

Um único polinômio de grau alto é um objeto global: um solavanco numa ponta pode sacudir a curva toda. A alternativa é grau baixo em muitos pedaços pequenos.

• A interpolação linear por partes é contínua com erro no máximo \(h^2\|f''\|_\infty/8\), pela fórmula do erro com \(n=1\). Estável (\(\Lambda=1\)), mas de suavidade só de primeira ordem.
• Os splines cúbicos são cúbicas por partes com duas derivadas contínuas. Para \(f\) suave o erro é \(O(h^4)\), e a constante de Lebesgue fica limitada quando a malha é refinada, ao contrário da interpolação polinomial global.

Uma regra prática: se \(f\) é analítica (ou muito suave) no intervalo todo e você pode escolher os nós, a interpolação global em pontos de Chebyshev converge geometricamente e vence. Se os dados vêm em pontos arbitrários, ou \(f\) tem quinas e saltos, ou você só precisa de poucos dígitos, pedaços de grau baixo são a escolha segura. No modo de dados livres, ponha doze pontos num zigue-zague suave: o polinômio global oscila entre eles e além deles, o que um spline não faria.

Se os \(y_j\) são medições com erros de tamanho \(\sigma\), a interpolação reproduz fielmente os erros, e pode amplificá-los até \(\Lambda_n\sigma\) entre os pontos. Com dados equiespaçados e \(n\) em torno de 20 isso é uma catástrofe (use o controle de ruído com “equiespaçados”, depois troque a família de nós para ver o contraste).

O que se faz na prática é ajuste por mínimos quadrados: escolher um polinômio de grau baixo \(m\ll n\) minimizando \(\sum_j(p(x_j)-y_j)^2\), ou um spline suavizante. O ajuste não passa pelos pontos, e isso é uma qualidade: ele faz a média do ruído em vez de persegui-lo. A interpolação é a ferramenta certa quando os dados são exatos (valores de uma função conhecida, ou nós que você escolhe) e a pergunta certa é “o que há entre os nós?”.

1. Escolha “Runge: 1/(1 + 25x²)”, nós equiespaçados, e aperte Aumentar n. O erro máximo primeiro fica perto de 0,4 (0,44 em \(n=4\), 0,43 em \(n=5\)) e depois explode: 0,43 em \(n=5\), 1,92 em \(n=10\), 59,8 em \(n=20\), \(2{,}4\cdot10^3\) em \(n=30\), \(1{,}05\cdot10^5\) em \(n=40\).
2. Marque “Ajustar a janela a p” para ver o polinômio inteiro, e abra a visão da função de Lebesgue no cartão da base.
3. Troque para Chebyshev: \(0{,}109\), \(0{,}0153\), \(2{,}9\cdot10^{-4}\) em \(n=10,20,40\). No cartão de convergência a curva de Chebyshev é uma reta de inclinação \(\log_{10}1{,}22^{-1}\) por unidade de \(n\).
4. Escolha “1/(1 + x²) em [−5, 5]” e compare: mesma função após reescalar, números idênticos.

Pergunta: por que o erro de Chebyshev em \(n=15\) supera o de \(n=12\)? (Lição sobre os nós de Chebyshev.)

Com nós de Chebyshev, percorra “|x|”, “√(x + 1)”, “bump suave”, “Runge” e “sen(2πx)” em \(n=10,20,40\) e compare com a tabela da lição do catálogo. No cartão de convergência observe as formas: uma curva quase achatada para a quina (\(O(1/n)\) só é reta num gráfico log–log, aqui parece uma curva suave), uma reta para Runge, uma curva despencando para o seno. Depois troque para nós equiespaçados com \(|x|\): a curva sobe depois de \(n\approx5\), que é a divergência de Bernstein.

Pergunta: a tabela sugere que dobrar \(n\) divide o erro por dois para \(|x|\). Que parte do teorema de Jackson é essa?

Escolha “polinômio x⁵ − 2x³ + x/2”. Para \(n\ge5\) o erro é só arredondamento (cerca de \(10^{-15}\)): unicidade. Para \(n=4\) não é zero, e pode ser previsto exatamente. A função é mônica de grau 5, então \(f^{(5)}/5!=1\) e a fórmula do erro diz \(f-p=\omega(x)\) sem \(\xi\) sobrando. Portanto o erro máximo é \(\|\omega\|_\infty\): para nós equiespaçados 0,113, e para nós de Chebyshev exatamente \(2^{-4}=0{,}0625\), que é o que o aplicativo mostra. Confira ambos na visão “Polinômio dos nós” do cartão Erro.

Escolha \(\sin(2\pi x)\), \(n=20\), ruído \(10^{-3}\). Nós equiespaçados: erro máximo cerca de 1,8, mil vezes o ruído, perto da folga \(\Lambda_{20}\approx1{,}1\cdot10^4\). Nós de Chebyshev: cerca de \(1{,}2\cdot10^{-3}\), o tamanho do ruído. Olhe o polinômio com “Ajustar a janela a p”: com nós equiespaçados o dano está no primeiro e no último subintervalo, exatamente onde a função de Lebesgue é grande. Tente ruído \(10^{-6}\) e \(10^{-2}\) e verifique que o dano escala linearmente com o ruído.

Escolha “Seus dados”. Arraste um ponto e veja o polinômio inteiro se mover: é um objeto global. Ponha pontos à direita do último: o polinômio lá longe à esquerda também muda. Ponha dois pontos bem juntos com alturas diferentes: o polinômio precisa de uma inclinação enorme entre eles e \(\Lambda_n\) dispara. Aperte Pôr um ponto aleatório várias vezes e veja o número de condição de Vandermonde nas métricas. No cartão Formas (para \(n\le8\)) acompanhe a tabela de diferenças divididas enquanto arrasta: a linha de cima, os coeficientes de Newton, só muda do ponto arrastado em diante para a direita.

1. Escolha “Seus dados” e ponha quatro pontos. No cartão A ideia de Lagrange escolha \(e_0\) e olhe o gráfico da esquerda: o polinômio vale \(1\) no primeiro ponto e \(0\) nos outros três. Faça o mesmo com \(e_1,e_2,e_3\).
2. Arraste o controle da direita de \(0\) até a última peça. Qual peça é a maior? Arraste o ponto correspondente para cima e veja aquela peça, e só ela, crescer.
3. Escolha “Sua função f(x)”, digite sin(3x) em \([-1,1]\) com \(n=3\) nós equiespaçados, e compare a soma das quatro peças com \(f\).
4. No script Python abaixo, troque NODES por "cheb" e compare as figuras lagrange_pieces.png das duas famílias: as peças \(y_j\ell_j\) são bem menores com nós de Chebyshev, que é a constante de Lebesgue em ação.

1. Ache o polinômio por \((-1,1),(0,0),(1,1)\) pela fórmula de Lagrange e por unicidade. Resposta: \(\ell_0+\ell_2=x^2\); e \(x^2\) tem grau 2 e serve, logo por unicidade é ele.
2. Mostre que \(\sum_j\ell_j(x)\equiv1\) e \(\sum_jx_j^k\ell_j(x)=x^k\) para \(k\le n\). Resposta: \(1\) e \(x^k\) são polinômios de grau no máximo \(n\), então cada um é o próprio interpolador, cuja forma de Lagrange é a soma dada.
3. Mostre que \(w_j=1/\omega'(x_j)\). Resposta: escreva \(\omega=(x-x_j)\prod_{k\ne j}(x-x_k)\); então \(\omega'(x_j)=\prod_{k\ne j}(x_j-x_k)\).
4. Limite o erro da interpolação linear de \(\ln x\) nos nós 1 e 2, em \(x=1{,}5\). Resposta: \(|f''|=1/x^2\le1\), \(|\omega(1{,}5)|=0{,}25\), então o erro é no máximo \(\tfrac12\cdot0{,}25=0{,}125\). O erro real é \(\ln1{,}5-\tfrac12\ln2=0{,}0589\).
5. Quantos nós de Chebyshev garantem erro abaixo de \(10^{-6}\) para \(e^x\) em \([-1,1]\)? Resposta: \(e/((n+1)!2^n)\lt 10^{-6}\) vale pela primeira vez em \(n=7\), isto é 8 nós.
6. Deduza a regra de Simpson de \(\int\ell_j\). Resposta: nos nós \(0,1,2\): \(\ell_0=\tfrac12(x-1)(x-2)\), \(\ell_1=x(2-x)\), \(\ell_2=\tfrac12x(x-1)\), com integrais em \([0,2]\) iguais a \(\tfrac13,\tfrac43,\tfrac13\). Logo \(\int_0^2f\approx\tfrac13(f_0+4f_1+f_2)\).
7. Prove que o coeficiente de Newton \(c_k\) depende só de \(x_0,\dots,x_k\) e \(y_0,\dots,y_k\). Resposta: \(c_k\) é o coeficiente líder do polinômio de grau no máximo \(k\) que interpola os primeiros \(k+1\) dados, pois os termos de Newton \(c_j\prod_{i\lt j}(x-x_i)\) com \(j\lt k\) têm grau abaixo de \(k\).

Para integrar \(f\) em \([a,b]\), escolha nós \(x_0,\dots,x_n\), seja \(p_n\) o polinômio de grau \(\le n\) pelos pontos \((x_j,f(x_j))\), e tome \(\int_a^b p_n\) como resposta. Na forma de Lagrange \(p_n=\sum f(x_j)\ell_j\), logo

Os pesos dependem só dos nós, nunca de \(f\). Uma regra de Newton–Cotes fechada usa nós igualmente espaçados \(x_j=a+jh\), \(h=(b-a)/n\), com os extremos incluídos. Substituindo \(x=a+th\) obtém-se \(w_j=h\,c_j\) com \(c_j=\int_0^n\prod_{k\ne j}\frac{t-k}{j-k}\,dt\), um número que se calcula uma vez. O painel “Pesos de Newton–Cotes fechados” os lista: \(n=1\) dá \(\tfrac12,\tfrac12\) (trapézio), \(n=2\) dá \(\tfrac13,\tfrac43,\tfrac13\) (Simpson 1/3), \(n=3\) dá \(\tfrac38,\tfrac98,\tfrac98,\tfrac38\) (Simpson 3/8).

Ligação com o aplicativo anterior: isto é exatamente “integrar o interpolador de Lagrange”. O erro da regra é a integral do erro de interpolação, e é assim que se obtém cada termo de erro abaixo.

Uma regra tem grau de precisão \(d\) se é exata para todo polinômio de grau \(\le d\) e não é para algum polinômio de grau \(d+1\). Basta testar nos monômios \(1,x,x^2,\dots\), por linearidade. O painel “Grau de precisão” faz isso em \([0,1]\).

Uma regra interpolatória com \(n+1\) nós tem \(d\ge n\), pois \(p_n=f\) quando \(\deg f\le n\). Para \(n\) par há mais um grau: seja \(m\) o ponto médio de \([a,b]\); os nós são simétricos em torno de \(m\) e os pesos também, logo \(\int (x-m)^{n+1}\) e a regra aplicada a ele se anulam (função ímpar, soma simétrica). Assim

trapézio \(n=1\): \(d=1\); Simpson 1/3 \(n=2\): \(d=3\), não 2; Simpson 3/8 \(n=3\): \(d=3\).

No aplicativo, escolha \(x^3\) em \([0,1]\), Simpson 1/3, simples: o erro é zero (abaixo de \(10^{-16}\)), embora a parábola por três pontos de uma cúbica não seja a cúbica. Para \(x^4\) os erros simples são \(-1/120\approx-8{,}3\cdot10^{-3}\) (Simpson 1/3) e \(-1/270\approx-3{,}7\cdot10^{-3}\) (Simpson 3/8); eles fixarão as constantes dos termos de erro.

Com \(n=1\) o interpolador é a corda, e sua integral é a área de um trapézio:

Escolha “Trapézio”, “Simples”, função \(\sin x\) em \([0,\pi]\): a corda é o eixo \(x\), a regra dá \(0\), a integral é \(2\), o erro é \(2\). Um único painel é uma ferramenta grosseira; o objetivo das próximas lições é quantificar o quanto, e curar isso cortando \([a,b]\) em muitos painéis.

Seja \(f\in C^2[a,b]\). O erro de interpolação por \(a,b\) é \(f(x)-p_1(x)=\tfrac12 f''(\xi_x)(x-a)(x-b)\). Integrando:

O fator \((x-a)(x-b)\) é \(\le0\) em todo \([a,b]\): não muda de sinal. Então o teorema do valor médio para integrais tira \(f''\) para fora, em algum \(\xi\in(a,b)\):

Lendo: \(f\) convexa (\(f''>0\)) dá \(E<0\), a corda fica acima do gráfico e a regra superestima. Exemplo no aplicativo: \(e^x\) em \([0,1]\), trapézio simples. A regra dá \(1{,}8591\), \(E=-0{,}14086\); a cota é \(e/12=0{,}2265\); resolvendo \(-\tfrac1{12}e^{\xi}=E\) obtém-se \(\xi=0{,}525\), dentro de \((0,1)\), como o teorema promete. O cartão “Termo de erro” refaz tudo isso para o caso escolhido.

Corte \([a,b]\) em \(n\) painéis de largura \(h=(b-a)/n\) e aplique a regra simples em cada um. Os nós interiores são contados duas vezes e os extremos uma:

Cada painel contribui \(-\tfrac{h^3}{12}f''(\xi_i)\). Somando, \(-\tfrac{h^3}{12}\sum_{i=1}^n f''(\xi_i)\). A média \(\frac1n\sum f''(\xi_i)\) fica entre o mínimo e o máximo da \(f''\) contínua, então pelo teorema do valor intermediário é igual a \(f''(\xi)\) para algum \(\xi\), e com \(nh=b-a\):

O erro é de ordem 2: dividir \(h\) por 2 divide o erro por 4. Confira no aplicativo com \(\sin x\) em \([0,\pi]\): para \(n=6,12,24,48\) os erros são \(4{,}59\cdot10^{-2}\), \(1{,}144\cdot10^{-2}\), \(2{,}857\cdot10^{-3}\), \(7{,}14\cdot10^{-4}\), e a tabela de convergência mostra ordem \(2{,}00\). No gráfico log–log é uma reta de inclinação \(-2\).

A forma do valor médio esconde o comportamento verdadeiro. Para \(f\) suave o erro tem uma expansão em potências pares de \(h\) (Euler–Maclaurin):

Só aparecem os extremos. Duas consequências. Primeiro, o termo principal é \(-\tfrac{h^2}{12}(f'(b)-f'(a))\), que é a integral de \(f''\) disfarçada, então o erro verdadeiro costuma ser um pouco menor que a cota com \(\max|f''|\). Para \(\sin x\) em \([0,\pi]\) o termo principal é \(h^2/6\), que é o erro observado com três dígitos.

Segundo, se \(f\) é periódica (ou só tem derivadas ímpares iguais em \(a\) e \(b\)) todos os termos se anulam e o trapézio converge mais depressa que qualquer potência de \(h\). Para \(e^{\cos x}\) em \([0,2\pi]\) o erro é \(3{,}4\cdot10^{-2}\) com \(n=4\), \(2{,}8\cdot10^{-4}\) com \(n=6\), \(1{,}3\cdot10^{-6}\) com \(n=8\), \(6{,}5\cdot10^{-12}\) com \(n=12\): nada de reta de inclinação \(-2\) aqui. Por isso o trapézio é o método de escolha para integrandos periódicos, e isso prepara o terreno para a quadratura gaussiana, o próximo tópico.

Como \(I-T_n=c_2h^2+c_4h^4+\cdots\) com as mesmas constantes para \(h\) e \(h/2\), a combinação que elimina o termo em \(h^2\) é

Escrevendo os pesos nos \(2n+1\) pontos, obtém-se \(\tfrac13(1,4,2,4,\dots,4,1)\) vezes o passo mais fino: esta é a regra de Simpson 1/3 composta. Numericamente, para \(\sin x\) em \([0,\pi]\): \(T_6\) e \(T_{12}\) dão \((4T_{12}-T_6)/3=2{,}0000526243\), idêntico a Simpson 1/3 com \(n=12\). Repetir a ideia no termo em \(h^4\) dá a integração de Romberg. A lição: o trapézio e Simpson são uma só família, e a ordem extra de Simpson vem de cancelar o termo principal da série de Euler–Maclaurin.

Tome três nós \(m-h,\,m,\,m+h\) (com \(m=\tfrac{a+b}2\), \(h=\tfrac{b-a}2\)) e a parábola por eles. Em vez de integrar polinômios de Lagrange, use o grau de precisão: a regra \(w_0f_0+w_1f_1+w_2f_2\) deve ser exata em \(1\) e em \((x-m)^2\), e por simetria \(w_0=w_2\). Então

Logo

Em \(\sin x\) sobre \([0,\pi]\) a regra simples dá \(2{,}0944\) (erro \(-0{,}0944\)) contra \(0\) do trapézio simples: três pontos e um ganho enorme. Geometricamente, integramos uma parábola pelo ponto médio em vez de uma corda.

O erro de Lagrange daria \(f'''\), mas a regra tem grau 3, então esperamos \(f^{(4)}\). O truque: seja \(q\) a cúbica com \(q=f\) em \(m-h,m,m+h\) e \(q'(m)=f'(m)\) (Hermite). A regra só enxerga os três valores, então \(Q(f)=Q(q)=\int q\), pois Simpson é exata em cúbicas. Logo \(E=\int(f-q)\), e o erro de Hermite é

Com \(t=x-m\), o peso \((t^2-h^2)t^2\) é \(\le0\) em \([-h,h]\): de novo sinal constante, então

Confira com o teste de \(x^4\): \(f^{(4)}=24\), \(h=\tfrac12\): \(-\tfrac{1}{32\cdot90}\cdot24=-\tfrac1{120}\), o valor da tabela de exatidão. Para \(e^x\) em \([0,1]\) a cota é \(\tfrac{e}{90\cdot32}=9{,}4\cdot10^{-4}\) e o erro real é \(-5{,}8\cdot10^{-4}\).

Tome \(n\) par e aplique a regra simples em cada um dos \(n/2\) pares de subintervalos (cada um de largura \(2h\)):

Cada par tem erro \(-\tfrac{h^5}{90}f^{(4)}(\xi_i)\); somando \(n/2\) deles e usando o valor intermediário como antes:

Ordem 4: dividir \(h\) por 2 divide o erro por 16. No aplicativo, \(e^{-x^2}\) em \([0,2]\) com \(n=6,12\): \(4{,}98\cdot10^{-5}\) e depois \(3{,}14\cdot10^{-6}\), razão \(15{,}9\). E a comparação de custo impressiona: para \(\sin x\) em \([0,\pi]\) com \(n=12\), trapézio e Simpson usam os mesmos 13 pontos, mas os erros são \(1{,}1\cdot10^{-2}\) e \(5{,}3\cdot10^{-5}\), um fator de 200.

Quatro nós, espaçamento \(h=\tfrac{b-a}3\), interpolador cúbico. A simetria dá \(w_0=w_3\), \(w_1=w_2\). Com \(m\) o ponto médio os nós são \(m\pm\tfrac{3h}2,\ m\pm\tfrac h2\). Exatidão em \(1\) e \((x-m)^2\):

Seu grau de precisão é 3 (não 4: o bônus da simetria vale só para \(n\) par), então \(E=K\,h^5f^{(4)}(\xi)\) com constante de sinal \(-\) (o núcleo de Peano da regra tem sinal constante, fato que citamos). A constante é fixada pelo teste \(f=x^4\): em \([0,1]\) a tabela dá \(E=-\tfrac1{270}\), com \(h=\tfrac13\), \(f^{(4)}=24\), então \(K\,\tfrac1{243}\cdot24=-\tfrac1{270}\), o que dá

A regra composta exige \(n\) múltiplo de 3, tem coeficientes \(1,3,3,2,3,3,2,\dots,3,3,1\), e erro \(-\tfrac{(b-a)h^4}{80}f^{(4)}(\xi)\) (há \(n/3\) grupos, cada um contribuindo \(-\tfrac{3h^5}{80}f^{(4)}\)).

As duas regras são de ordem 4; as constantes dos erros compostos são \(\tfrac1{180}\) e \(\tfrac1{80}\). Com o mesmo \(n\) (logo o mesmo \(h\) e o mesmo número de avaliações) o erro da 3/8 é \(\tfrac{180}{80}=2{,}25\) vezes maior. Observado com \(\sin x\) em \([0,\pi]\), \(n=12\): \(5{,}26\cdot10^{-5}\) para a 1/3 contra \(1{,}19\cdot10^{-4}\) para a 3/8, razão \(2{,}27\).

Logo, quando \(n\) é par, use a 1/3. A 3/8 se justifica quando \(n\) é ímpar: tome \(n-3\) (par) subintervalos com Simpson 1/3 e os três últimos com Simpson 3/8; a regra combinada tem a mesma ordem. Exemplo: \(n=9\) é uma 3/8 mais três painéis 1/3. As duas regras são também os primeiros membros de um padrão: a regra de Boole (\(n=4\), grau 5) é a próxima com \(n\) par.

Para garantir \(|E|\le\varepsilon\) resolva a cota para \(n\). Tome \(\int_0^\pi\sin x\,dx\) com \(\varepsilon=10^{-8}\) e \(\max|f''|=\max|f^{(4)}|=1\):

Experimentos diretos mostram o quanto a cota é pessimista: o trapézio atinge \(10^{-8}\) com cerca de \(12\,830\) subintervalos (o termo de Euler–Maclaurin \(h^2/6\) prevê \(\pi/\sqrt{6\cdot10^{-8}}\approx12\,826\)) e Simpson com \(102\). A razão \(16\,075/114\) já diz por que Simpson é a escolha padrão para integrandos suaves. Na prática faz-se também o inverso: dobra-se \(n\) até que dois valores sucessivos concordem na tolerância.

Os termos de erro supõem \(f\in C^2\) ou \(C^4\). Teste o aplicativo em duas funções que violam isso.

\(\sqrt x\) em \([0,1]\): \(f''\) explode em 0, então a cota é infinita. Ordens medidas: \(1{,}46\) a \(1{,}49\) para o trapézio, \(1{,}50\) para as duas regras de Simpson. A singularidade \(x^{1/2}\) no extremo limita toda regra à ordem \(3/2\), e a precisão extra de Simpson é desperdiçada. Remédio: substituir \(x=t^2\) para eliminar a singularidade, ou dividir o intervalo.

\(|x-1/\pi|\) em \([0,1]\): uma quina dentro do intervalo. O erro é irregular porque depende de onde a quina cai em relação aos nós (uma quina exatamente sobre um nó daria uma resposta enganosamente exata, por isso a quina está no irracional \(1/\pi\)). A tabela de convergência mostra ordens que saltam e podem até ficar negativas. A lição: uma regra de ordem alta precisa de suavidade em todo o intervalo; ponha um nó em cada quina, isto é, integre por partes.

As ordens 2 e 4 valem para uma \(f\) suave genérica; funções especiais podem render mais. Tome \(1/(1+x^2)\) em \([0,1]\), cuja integral é \(\pi/4\). O trapézio tem ordem \(2{,}00\), como esperado (\(f'(1)=-\tfrac12\ne0\)). Mas Simpson 1/3 e 3/8 mostram ordem \(6{,}0\): erros \(2{,}2\cdot10^{-7}\), \(3{,}3\cdot10^{-9}\), \(5{,}2\cdot10^{-11}\) para \(n=6,12,24\).

A razão está na série de Euler–Maclaurin de Simpson, cujo coeficiente de \(h^4\) é proporcional a \(f'''(b)-f'''(a)\). Aqui \(f'''(x)=-24x(1-x^2)/(1+x^2)^4\) se anula em \(0\) e em \(1\), então o termo em \(h^4\) some e \(h^6\) lidera. A cota de erro com \(f^{(4)}\) continua verdadeira, só que longe de justa. Consequência prática: estime a ordem a partir dos dados (a tabela de convergência), não confie cegamente no número do livro. De bônus, \(4\,S_6\) com este integrando é \(3{,}14159178\), seis dígitos corretos de \(\pi\) com sete valores da função.

Por que parar em painéis de \(n=3\)? Porque as regras fechadas maiores herdam a fraqueza da interpolação em nós igualmente espaçados (o fenômeno de Runge do aplicativo anterior). A tabela de pesos mostra isso. Para \(n=6\) (vizinha de Boole, grau 7) todos os pesos são positivos. A partir de \(n=8\) os pesos mudam de sinal:

e para \(n=10\) o maior peso é \(7{,}14\) e a quantidade \(\sum|c_j|/n\) cresceu de \(1\) para \(3{,}06\) (já é \(1{,}45\) em \(n=8\)). Como \(\sum c_j=n\) sempre, pesos grandes positivos e negativos se cancelam, e qualquer erro de arredondamento em \(f\) é amplificado por esse fator: a regra fica instável. A alternativa estável é o que fizemos: manter uma regra de grau baixo (trapézio, Simpson) e refinar \(h\). A única maneira de obter grau alto sem perder estabilidade é abandonar os nós igualmente espaçados, assunto da próxima aula, a quadratura gaussiana.

Aula 1: da interpolação ao trapézio.

• Lições 1.1–1.2: quadratura interpolatória, pesos \(w_j=\int\ell_j\), grau de precisão, bônus da simetria (abra os painéis de pesos e de exatidão).
• Lições 2.1–2.3: trapézio simples, demonstração do erro com o teorema do valor médio para integrais, regra composta e ordem 2. No aplicativo: seno, trapézio, simples e depois composta, leia a ordem na tabela de convergência.
• Se sobrar tempo, 2.4 (Euler–Maclaurin e funções periódicas) como observação.

Aula 2: Simpson e prática.

• Lições 3.1–3.3: Simpson 1/3, a demonstração de Hermite de \(-h^5f^{(4)}/90\), composta de ordem 4. A lição 2.5 (Richardson) liga com a primeira aula.
• Lições 3.4–3.5: Simpson 3/8, comparação a custo igual, a combinação para \(n\) ímpar.
• Lições 4.1–4.4 como discussão: escolha de \(n\), funções ásperas (\(\sqrt x\), a quina), a afortunada \(1/(1+x^2)\), pesos negativos como ponte para a quadratura gaussiana.

1. Deduza os pesos \(\tfrac38,\tfrac98,\tfrac98,\tfrac38\) resolvendo as equações dos momentos para \(1,x,x^2,x^3\) em \([0,3]\). Resposta: o sistema é o do painel de pesos; a solução é \(c=(\tfrac38,\tfrac98,\tfrac98,\tfrac38)\) em unidades de \(h\), somando \(3\).
2. Para \(\int_0^1e^x\,dx\) com o trapézio composto e \(n=10\), dê a cota e compare com o aplicativo. Resposta: cota \(e/1200=2{,}27\cdot10^{-3}\); real \(-1{,}43\cdot10^{-3}\). Com Simpson 1/3, \(n=10\): cota \(e/(180\cdot10^4)=1{,}51\cdot10^{-6}\); real \(-9{,}5\cdot10^{-7}\).
3. Mostre que Simpson 1/3 integra \(x^3\) exatamente em qualquer \([a,b]\). Resposta: escreva \(x^3=(x-m)^3+3m(x-m)^2+3m^2(x-m)+m^3\); as potências ímpares de \(x-m\) se cancelam por simetria e o resto tem grau \(\le2\).
4. Estime \(\pi=4\int_0^1\frac{dx}{1+x^2}\) com \(n=6\). Resposta: trapézio \(3{,}13696\) (erro \(4{,}6\cdot10^{-3}\)); Simpson 1/3 \(3{,}14159178\) (erro \(8{,}7\cdot10^{-7}\)).
5. Integre \(e^x\) em \([0,1]\) com \(n=9\) usando Simpson 1/3 nos seis primeiros subintervalos e Simpson 3/8 nos três últimos, e depois o contrário. Resposta: erros \(-2{,}27\cdot10^{-6}\) e \(-1{,}87\cdot10^{-6}\); a 3/8 composta pura com \(n=9\) dá \(-3{,}26\cdot10^{-6}\). Pôr o bloco 3/8 (de erro maior) onde \(f^{(4)}\) é menor, perto de \(x=0\), é melhor.
6. Explique por que o trapézio com \(n=8\) em \(e^{\cos x}\), \([0,2\pi]\), tem precisão de cerca de \(10^{-6}\) enquanto em \(e^x\), \([0,1]\), precisa de \(n\approx 380\) (trapézio) para o mesmo. Resposta: lição 2.4, periodicidade.

Uma regra \(Q(f)=\sum_{j=1}^n w_jf(x_j)\) para \(\int f\,w\) tem \(2n\) números livres: \(n\) pesos e \(n\) nós. Newton–Cotes fixa os nós e usa só os pesos, e chega ao grau \(n-1\) (ou \(n\)). Se os nós também são livres, podemos exigir exatidão em \(1,x,\dots,x^{2n-1}\), que são \(2n\) equações:

n = 1 (Legendre, \([-1,1]\)): \(w_1=2\) e \(w_1x_1=0\), logo \(x_1=0\): a regra do ponto médio, uma avaliação, grau 1. n = 2: por simetria \(w_1=w_2=1\), e \(2x^2=\mu_2=\tfrac23\) dá \(x=\pm\tfrac1{\sqrt3}=\pm0{,}5774\). Duas avaliações e grau 3, o mesmo que Simpson usa com três. O sistema é não linear nos nós, e por isso as próximas lições precisam de estrutura.

Nenhuma regra faz melhor que \(2n-1\). Tome \(q(x)=\prod_j(x-x_j)^2\), de grau \(2n\): \(q\ge0\), então \(\int q\,w>0\), enquanto \(Q(q)=0\).

Teorema. Sejam \(x_1,\dots,x_n\) distintos, \(\omega(x)=\prod(x-x_j)\), e sejam os pesos os da quadratura interpolatória (\(w_j=\int\ell_jw\)). A regra tem grau de precisão \(\ge2n-1\) se e somente se \(\int\omega\,q\,w=0\) para todo polinômio \(q\) de grau \(\le n-1\).

Demonstração. (⇐) Seja \(\deg f\le2n-1\). Divida: \(f=q\,\omega+r\) com \(\deg q,\deg r\le n-1\). Então \(\int fw=\int q\omega w+\int rw=\int rw\) por ortogonalidade. Do outro lado, \(f(x_j)=r(x_j)\) pois \(\omega(x_j)=0\), logo \(Q(f)=Q(r)=\int rw\), o último passo porque a regra é interpolatória e \(\deg r\le n-1\). (⇒) Tome \(f=q\,\omega\) com \(\deg q\le n-1\), então \(\deg f\le2n-1\): a exatidão dá \(\int q\omega w=Q(q\omega)=0\), pois \(\omega\) se anula nos nós. ∎

Assim o problema de escolher os nós vira: achar o polinômio mônico \(\omega\) de grau \(n\) ortogonal a todos os graus menores, e tomar suas raízes. Esse polinômio é \(p_n\), o \(n\)-ésimo polinômio ortogonal para o peso \(w\).

Para um peso \(w>0\) em um intervalo (com momentos finitos) ponha \(\langle f,g\rangle=\int fg\,w\). Aplicando Gram–Schmidt a \(1,x,x^2,\dots\) obtém-se uma única sequência de polinômios mônicos \(p_k\), \(\deg p_k=k\), com \(\langle p_j,p_k\rangle=0\) para \(j\ne k\). Equivalentemente, \(p_k\perp\) todos os polinômios de grau \(\lt k\) (eles geram o mesmo espaço que \(p_0,\dots,p_{k-1}\)).

A tabela de Gram do aplicativo confere isso numericamente: para Legendre a diagonal é \(2,\ 0{,}6667,\ 0{,}4,\ 0{,}2857,\ 0{,}2222\), isto é, \(2/(2k+1)\), e toda entrada fora da diagonal é zero a menos de arredondamento.

Os polinômios ortogonais satisfazem, para \(p_k\) mônicos,

Demonstração. \(xp_k-p_{k+1}\) tem grau \(\le k\), então \(xp_k=p_{k+1}+\sum_{j\le k}c_jp_j\) com \(c_j=\langle xp_k,p_j\rangle/h_j\). Para \(j\lt k-1\): \(\langle xp_k,p_j\rangle=\langle p_k,xp_j\rangle=0\), porque \(\deg(xp_j)=j+1\lt k\). Só sobrevivem \(j=k\) e \(j=k-1\). ∎

Os coeficientes são explícitos: \(\alpha_k=\langle xp_k,p_k\rangle/h_k\) e \(\beta_k=h_k/h_{k-1}\). Nas normalizações usuais as recorrências são, para Legendre \((k+1)P_{k+1}=(2k+1)xP_k-kP_{k-1}\); Chebyshev \(T_{k+1}=2xT_k-T_{k-1}\); Laguerre \((k+1)L_{k+1}=(2k+1-x)L_k-kL_{k-1}\); Hermite \(H_{k+1}=2xH_k-2kH_{k-1}\). O código desta página avalia as quatro com um só laço; a recorrência é também a melhor forma de calcular \(p_n(x)\) (é barata e estável).

Afirmação. \(p_n\) tem \(n\) raízes reais distintas, todas no intervalo aberto \((a,b)\).

Demonstração. Sejam \(t_1,\dots,t_m\) os pontos de \((a,b)\) onde \(p_n\) muda de sinal (raízes de multiplicidade ímpar); \(m\le n\). Ponha \(q=\prod_{i\le m}(x-t_i)\). Então \(p_nq\) não muda de sinal em \((a,b)\), logo \(\int p_nq\,w\ne0\). Se \(m\lt n\), a ortogonalidade de \(p_n\) a polinômios de grau \(\lt n\) daria \(0\): contradição. Portanto \(m=n\): as \(n\) raízes são reais, simples e estão em \((a,b)\). ∎

É isto que torna a regra utilizável: os nós são reais, dentro do domínio de integração, nunca repetidos. Veja o gráfico: para Legendre com \(n=4\) as raízes estão em \(\pm0{,}3400\) e \(\pm0{,}8611\), e os nós se aglomeram perto dos extremos de \([-1,1]\) quando \(n\) cresce (experimente \(n=30\)); para Laguerre e Hermite as raízes se espalham em \([0,\approx4n]\) e \([-2\sqrt n,2\sqrt n]\).

Entre duas raízes consecutivas de \(p_n\) há exatamente uma raiz de \(p_{n-1}\), e \(p_{n-1}\) não tem raiz fora de \([x_1,x_n]\): as raízes se entrelaçam. O gráfico mostra: os pontos pequenos (raízes de \(p_{n-1}\)) ficam nos intervalos entre os pontos grandes (raízes de \(p_n\)). Isso segue da identidade de Christoffel–Darboux, ou da recorrência por indução: \(p_{n+1}(t_i)=-\beta_np_{n-1}(t_i)\) em uma raiz \(t_i\) de \(p_n\), então \(p_{n+1}\) e \(p_{n-1}\) têm sinais opostos nas raízes de \(p_n\), o que força uma mudança de sinal de \(p_{n+1}\) em cada intervalo.

Uso prático: para achar todas as raízes de \(p_n\) precisa-se de um intervalo para cada uma. Os intervalos entre as raízes de \(p_{n-1}\), mais os dois intervalos das pontas, são \(n\) intervalos disjuntos, cada um com exatamente uma raiz. O código calcula as linhas de raízes de \(p_1,p_2,\dots,p_n\) sucessivamente e em cada intervalo roda o método de Newton mantido dentro do intervalo (bisseção quando um passo de Newton escapa). Nunca se precisa de chute inicial e nenhuma raiz pode ser perdida ou achada duas vezes, que é a falha usual dos localizadores gerais de raízes polinomiais em grau 20 ou 30.

Com os nós fixados, os pesos são os da quadratura interpolatória, \(w_j=\int\ell_j\,w\), onde \(\ell_j\) é o polinômio de base de Lagrange. Dois fatos:

Positividade. \(\ell_j^2\) tem grau \(2n-2\le2n-1\), então a regra o integra exatamente: \(0\lt \int\ell_j^2w=\sum_kw_k\ell_j(x_k)^2=w_j\), pois \(\ell_j(x_k)=\delta_{jk}\). Todo peso é positivo, e \(\sum w_j=\int w=\mu_0\).

Fórmula de Christoffel. \(w_j=\dfrac{k_n\,h_{n-1}}{k_{n-1}\,p_n'(x_j)\,p_{n-1}(x_j)}\), com \(k_m\) o coeficiente líder de \(p_m\). É o que o código usa, com \(p_n'\) obtido derivando a recorrência. Confira no aplicativo (cartão Nós e pesos): Legendre \(n=3\) dá pesos \(\tfrac59,\tfrac89,\tfrac59\) em \(0,\pm0{,}7746\); Chebyshev dá \(\pi/n\) em todo nó.

Pesos positivos são a diferença profunda em relação a Newton–Cotes de ordem alta (cujos pesos, a partir de \(n=8\), mudam de sinal): sem cancelamento, a regra é numericamente estável para todo \(n\).

Seja \(f\in C^{2n}\) e seja \(q\) o interpolador de Hermite: grau \(\le2n-1\), \(q(x_j)=f(x_j)\), \(q'(x_j)=f'(x_j)\). A regra só enxerga os valores, então \(Q(f)=Q(q)=\int qw\) (exata, grau \(2n-1\)). Assim \(E=\int(f-q)\,w\), e a fórmula do erro de Hermite é

Agora \(\omega^2\ge0\): sinal constante, como na demonstração de Simpson. Pelo teorema do valor médio para integrais

Formas fechadas: Legendre \(\frac{(b-a)^{2n+1}(n!)^4}{(2n+1)((2n)!)^3}\); Chebyshev \(\frac{2\pi}{2^{2n}(2n)!}\); Laguerre \(\frac{(n!)^2}{(2n)!}\); Hermite \(\frac{n!\sqrt\pi}{2^n(2n)!}\). Para \(n=2\), Legendre: \(E=f^{(4)}(\xi)/135\) em \([-1,1]\). Confira com \(f=x^4\) (\(f^{(4)}=24\)): a regra dá \(2\cdot\tfrac19=\tfrac29\), a integral é \(\tfrac25\), a diferença \(\tfrac8{45}=\tfrac{24}{135}\). ✓. Com \(e^x\): a constante para \(n=2\) é \(1/135\), a cota \(e/135=2{,}0\cdot10^{-2}\), o erro real \(7{,}7\cdot10^{-3}\); para \(n=3\) a constante é \(6{,}35\cdot10^{-5}\) e o erro \(6{,}5\cdot10^{-5}\); o aplicativo imprime ambos, e o valor implícito \(f^{(6)}(\xi)=1{,}03\in[1,e]\).

Como os pesos são positivos e \(\sum w_j=\mu_0\), para \(f\) contínua em um intervalo finito e qualquer polinômio \(r\) de grau \(\le2n-1\): \(|I-Q|\le|\int(f-r)w|+|Q(r-f)|\le2\mu_0\max|f-r|\). Escolhendo o melhor \(r\) (Weierstrass), a quadratura de Gauss converge para toda \(f\) contínua, e tão depressa quanto a melhor aproximação polinomial de grau \(2n-1\): geometricamente para funções analíticas. Newton–Cotes não pode prometer isso (pesos negativos).

Números (Legendre, \(e^x\) em \([-1,1]\), mesmo número de avaliações): 3 pontos: Gauss \(6{,}5\cdot10^{-5}\), Simpson \(1{,}2\cdot10^{-2}\); 5 pontos: Gauss \(8{,}2\cdot10^{-10}\), Simpson \(7{,}9\cdot10^{-4}\); com \(n=6\) Gauss já está em \(1{,}6\cdot10^{-12}\), enquanto Simpson com 11 pontos está em \(2{,}1\cdot10^{-5}\). Outro exemplo, do aplicativo de Newton–Cotes: \(\int_0^\pi\sin x\,dx\) com três avaliações: Simpson \(2{,}0944\) (erro \(-0{,}094\)), Gauss–Legendre \(2{,}00139\) (erro \(-1{,}4\cdot10^{-3}\)). Abra o cartão de convergência: a curva de Gauss se curva para baixo mais depressa que qualquer reta, enquanto a de Simpson decai como \(n^{-4}\).

Três maneiras de obter a regra para um dado \(n\). (1) Tabelas, boas para \(n\le5\) à mão. (2) O método desta página: raízes de \(p_n\) por Newton protegido dentro dos intervalos de entrelaçamento, pesos pela fórmula de Christoffel. (3) Golub–Welsch (1969): a recorrência diz \(xp_k=p_{k+1}+\alpha_kp_k+\beta_kp_{k-1}\), isto é, os nós são os autovalores da matriz de Jacobi tridiagonal simétrica com diagonal \(\alpha_k\) e fora da diagonal \(\sqrt{\beta_k}\), e \(w_j=\mu_0\,v_{j,1}^2\) onde \(v_{j,1}\) é a primeira componente do autovetor normalizado. É o padrão nas bibliotecas.

Limites da precisão dupla. Os nós e pesos são precisos perto do arredondamento para todo \(n\le30\) aqui (Legendre \(n=20\): \(\sum w_j=2\) a \(4\cdot10^{-14}\)). Os pesos de Laguerre e Hermite decaem super-exponencialmente: para \(n=10\) o menor peso de Laguerre é \(9{,}9\cdot10^{-13}\) no nó \(29{,}9\), e o menor peso de Hermite é \(7{,}6\cdot10^{-6}\). É por isso que a tabela de “grau de precisão” só é mostrada para \(n\le10\): os monômios \(x^k\), \(k\) perto de \(2n\), são uma base ruim e o próprio teste, não a regra, perde dígitos.

Peso \(1\) em \([-1,1]\); raízes do polinômio de Legendre \(P_n\). Casos pequenos explícitos:

Para \([a,b]\): \(x=\frac{b-a}2t+\frac{a+b}2\), \(\int_a^bf=\frac{b-a}2\sum w_jf(x(t_j))\), e a constante do erro vira \(\frac{(b-a)^{2n+1}(n!)^4}{(2n+1)((2n)!)^3}\) (para \(n=2\): \((b-a)^5/4320\)). É a regra para todo uso com integrandos suaves em intervalo finito. No aplicativo, Legendre, \(e^x\), \(n=1,2,3,4,5\): erros \(3{,}5\cdot10^{-1}\), \(7{,}7\cdot10^{-3}\), \(6{,}5\cdot10^{-5}\), \(2{,}9\cdot10^{-7}\), \(8{,}2\cdot10^{-10}\): cada nó multiplica a precisão por um fator entre \(45\) e \(360\). A função de teste \(x^6\) é exata para \(n\ge4\) (grau \(7\ge6\)) e tem erro \(4{,}6\cdot10^{-2}\) para \(n=3\), que é \(720\cdot6{,}35\cdot10^{-5}\), exatamente o termo de erro, já que \(f^{(6)}=720\).

Peso \((1-x^2)^{-1/2}\) em \([-1,1]\). Aqui as raízes de \(T_n\) são explícitas, porque \(T_n(\cos\theta)=\cos n\theta\):

Por que pesos iguais: com \(x=\cos\theta\), \(\int_{-1}^1\frac{f(x)}{\sqrt{1-x^2}}dx=\int_0^\pi f(\cos\theta)\,d\theta\), e a regra é a do ponto médio em \(\theta\) para a função par e periódica \(f(\cos\theta)\), o caso em que as regras do trapézio e do ponto médio são espectralmente precisas (lição de Euler–Maclaurin do aplicativo de Newton–Cotes). Use-a em integrais com o fator \(1/\sqrt{1-x^2}\), que torna inofensivas as singularidades nos extremos. Exemplo no aplicativo: \(\int_{-1}^1\frac{e^x}{\sqrt{1-x^2}}dx=\pi I_0(1)=3{,}97746326\); \(n=4\) dá \(3{,}97746264\), erro \(6{,}3\cdot10^{-7}\), abaixo da cota \(6{,}09\cdot10^{-7}\cdot e=1{,}7\cdot10^{-6}\); para \(n=2\): nós \(\pm0{,}7071\), valor \(\pi\cosh0{,}7071=3{,}9603\), erro \(1{,}7\cdot10^{-2}\).

Peso \(e^{-x}\) em \([0,\infty)\): \(\int_0^\infty e^{-x}f(x)\,dx\approx\sum w_jf(x_j)\), \(E=\frac{(n!)^2}{(2n)!}f^{(2n)}(\xi)\). Para \(n=2\): nós \(2\mp\sqrt2=0{,}5858,\ 3{,}4142\), pesos \(\frac{2\pm\sqrt2}4=0{,}8536,\ 0{,}1464\). A regra tem grau 3: \(\int_0^\infty x^3e^{-x}=6\) é reproduzida (\(5{,}99999\ldots\)) enquanto para \(x^4\) dá \(20\) em vez de \(24\). A função de teste \(x^5\) é exata para \(n\ge3\).

Use-a quando o integrando decai como \(e^{-x}\) vezes algo suave. Para integrar \(\int_0^\infty g\) tome \(f=ge^{x}\), o que só faz sentido se \(g\) já decai pelo menos como \(e^{-x}\). No aplicativo: \(\int_0^\infty e^{-x}\sin x\,dx=\tfrac12\): o erro oscila de sinal e cai a \(-1{,}9\cdot10^{-8}\) para \(n=12\) e \(1{,}5\cdot10^{-11}\) para \(n=16\). Com uma função que tem singularidade no extremo, \(\sqrt x\) (exato \(\sqrt\pi/2\)) ou \(\ln x\) (exato \(-\gamma\)), a convergência é só algébrica: para \(\sqrt x\) o erro é \(1{,}3\cdot10^{-2},\ 4{,}6\cdot10^{-3},\ 1{,}6\cdot10^{-3}\) para \(n=4,8,16\), um fator \(2^{3/2}\) a cada dobra.

Peso \(e^{-x^2}\) em \((-\infty,\infty)\), \(E=\frac{n!\sqrt\pi}{2^n(2n)!}f^{(2n)}(\xi)\). Casos pequenos:

Os nós são simétricos, então funções ímpares são integradas exatamente para qualquer \(n\). Com \(n=2\): \(\int x^2e^{-x^2}=\tfrac{\sqrt\pi}2\) é exata (grau 3), \(x^4\) dá \(0{,}4431\) contra \(\tfrac{3\sqrt\pi}4=1{,}3293\). Teste \(x^6\) (exato \(\tfrac{15\sqrt\pi}8\)): erro \(1{,}33\) para \(n=3\), zero para \(n\ge4\). Exemplo com função transcendente: \(\int e^{-x^2}\cos x\,dx=\sqrt\pi\,e^{-1/4}=1{,}38038845\); os erros para \(n=4,6,8\) são \(5{,}9\cdot10^{-5}\), \(3{,}7\cdot10^{-8}\), \(1{,}2\cdot10^{-11}\), e para \(n=12\) nível de arredondamento. Use Gauss–Hermite para esperanças sob uma distribuição normal (com a reescala \(x\mapsto\sqrt2\,x\)), integrais em mecânica quântica e qualquer integrando com fator gaussiano. Para \(1/(1+x^2)\) (exato \(1{,}3433\)) a convergência é lenta e algébrica: os polos em \(\pm i\) estão perto demais do eixo real.

• A suavidade decide. Para \(|x|\) em \([-1,1]\) (uma quina no meio) os erros para \(n=2,4,8,16,32\) são \(1{,}5\cdot10^{-1}\), \(4{,}3\cdot10^{-2}\), \(1{,}2\cdot10^{-2}\), \(3{,}0\cdot10^{-3}\), \(7{,}8\cdot10^{-4}\): só ordem \(n^{-2}\), a mesma do trapézio composto. Divida o intervalo na quina e use Gauss em cada pedaço. Para \(1/(1+25x^2)\) os erros são \(3{,}4\cdot10^{-1}\), \(1{,}8\cdot10^{-1}\), \(4{,}1\cdot10^{-2}\), \(1{,}8\cdot10^{-3}\), \(3{,}1\cdot10^{-6}\): Gauss converge aqui (a interpolação em nós equiespaçados não), mas tarde, pois os polos em \(\pm i/5\) estão perto do intervalo.
• Oscilação exige nós. Para \(\cos10x\) em \([-1,1]\): \(n=2,4\) são lixo (\(-1{,}85\), \(1{,}63\)), depois \(5{,}3\cdot10^{-3}\) em \(n=8\), \(3{,}8\cdot10^{-5}\) em \(n=10\), \(1{,}6\cdot10^{-13}\) em \(n=16\). Regra prática: cerca de \(\pi\) nós por comprimento de onda antes de começar o regime geométrico.
• Gauss composta. Para intervalo longo ou função áspera, corte em painéis e use um \(n\) fixo (3 a 7) em cada um, como Simpson composta mas com ordem \(2n\). Os nós de Gauss não são aninhados, então refinar \(n\to2n\) não reaproveita avaliações; as regras de Gauss–Kronrod e de Clenshaw–Curtis resolvem isso e são o que as bibliotecas adaptativas usam.
• Escolha a família pelo peso, não pelo intervalo. Se o integrando é \(e^{-x^2}g(x)\) com \(g\) suave, Gauss–Hermite com \(f=g\) é exponencialmente melhor que truncar a integral e usar Legendre.

Aula 1: por que os nós são raízes de polinômios ortogonais.

• 1.1–1.2: contagem de parâmetros, grau \(2n-1\), o teorema da ortogonalidade (a demonstração é curta e é o coração do assunto).
• 2.1–2.4: produto interno, Gram–Schmidt, recorrência de três termos, raízes reais e simples, entrelaçamento. No aplicativo: mude \(n\) e observe as raízes e os pontos pequenos; leia a tabela de Gram.
• 3.1 (pesos, positividade) e 4.1 (Gauss–Legendre, \(n=2,3\) à mão, depois o aplicativo).

Aula 2: erro, as outras três famílias, prática.

• 3.2 (demonstração do termo de erro por interpolação de Hermite; reaproveite a demonstração de Simpson da aula anterior), 3.3 (convergência e Gauss contra Simpson no cartão de convergência).
• 4.2–4.4: Chebyshev (nós explícitos), Laguerre, Hermite: para cada uma, o peso, um cálculo explícito com \(n=2\) e um experimento no aplicativo.
• 3.4 e 5.1 como discussão: Golub–Welsch, Gauss composta, funções ásperas.

1. Deduza a regra de Gauss–Legendre com \(n=2\) a partir das quatro equações de momentos. Resposta: \(w_1=w_2=1\), \(x=\pm1/\sqrt3\) (lição 1.1).
2. Mostre que ela não é exata para \(x^4\) e calcule o erro de duas maneiras. Resposta: \(\int x^4=\tfrac25\), \(Q=\tfrac29\), \(E=\tfrac8{45}\); a fórmula dá \(24/135=\tfrac8{45}\).
3. Estime \(\int_0^\pi\sin x\,dx\) com Gauss–Legendre \(n=2\) e \(n=3\) (mude as variáveis para \([-1,1]\)). Resposta: \(1{,}93582\) (erro \(6{,}4\cdot10^{-2}\)) e \(2{,}00139\) (erro \(-1{,}4\cdot10^{-3}\)).
4. Com a regra de Chebyshev, \(n=2\), estime \(\int_{-1}^1\frac{e^x}{\sqrt{1-x^2}}dx\). Resposta: \(\pi\cosh(1/\sqrt2)=3{,}9603\); o valor exato é \(3{,}97746\), erro \(1{,}7\cdot10^{-2}\).
5. Verifique que a regra de Laguerre com \(n=2\) é exata para \(x^3\) mas não para \(x^4\). Resposta: \(0{,}8536\cdot0{,}5858^3+0{,}1464\cdot3{,}4142^3=6\); para \(x^4\) dá \(20\neq24\).
6. Mostre que Gauss–Hermite com \(n=2\) integra \(x^2e^{-x^2}\) exatamente. Resposta: \(2\cdot\tfrac{\sqrt\pi}2\cdot\tfrac12=\tfrac{\sqrt\pi}2\).
7. Prove que nenhuma regra de quadratura com \(n\) nós é exata para todos os polinômios de grau \(2n\). Resposta: \(q=\prod(x-x_j)^2\) (lição 1.1).

A ponte Tacoma Narrows, no estado de Washington, tinha vão principal de 2.800 pés (\(853\,\mathrm{m}\)) e um tabuleiro com apenas cerca de 39 pés (\(12\,\mathrm{m}\)) de largura, esbelto para o comprimento. Foi aberta em 1.º de julho de 1940 e, desde os primeiros dias, balançava com vento moderado. Em 7 de novembro de 1940 as oscilações verticais deram lugar a um grande movimento de torção, com as duas bordas do tabuleiro subindo e descendo em oposição, e em cerca de uma hora o vão principal caiu. Ninguém morreu; um cachorro dentro de um carro foi a única vítima. O vento era de cerca de 64 km/h (40 mph). As fontes divergem em alguns números (relatos antigos dão o dia 10 de novembro), então trate os valores como aproximados.

Não no sentido dos livros-texto. A história popular diz que o vento empurrava na frequência natural da ponte, como uma criança num balanço. Billah e Scanlan (1991) mostraram que livros de física repetiram isso e que não é o que a evidência sustenta. Os vórtices desprendidos pelo tabuleiro tinham frequência bem diferente da frequência de torção que apareceu. A explicação aceita pela maioria dos engenheiros é o flutter aeroelástico: quando o tabuleiro torce, o escoamento do ar ao redor muda de modo a empurrá-lo ainda mais, uma realimentação autoexcitada que se comporta como amortecimento negativo acima de uma velocidade crítica do vento. A ressonância exige um ritmo externo; o flutter cria o seu próprio.

Outros mecanismos foram propostos, e levantamentos publicados afirmam que nenhuma explicação de todos os detalhes é universalmente aceita. Este aplicativo trata de um deles, não do veredito.

Porque é um exemplo limpo de mecanismo não linear que acopla dois movimentos que a teoria linear mantém separados. Lazer e McKenna (1990) observaram que um cabo estica mas não comprime, então um pendural é uma mola que se desliga quando frouxa. McKenna e Tuama (2001) mostraram num modelo simples que a força vertical pode então criar uma resposta de torção. Críticos, como Green e Unruh, argumentam que os pendurais de Tacoma de fato não afrouxaram e que os parâmetros necessários são irreais. Portanto trate o cartão como um laboratório para o mecanismo e para a numérica, não como reconstituição de 7 de novembro de 1940.

A seção do tabuleiro é uma barra rígida de meia largura \(\ell\), massa \(m\) por metro. Seu centro desce \(y\) e ela gira \(\theta\), de modo que as bordas esquerda e direita descem \(y\mp\ell\sin\theta\). Cada pendural é alongado por \(f=y_0+y\mp\ell\sin\theta\), em que \(y_0=mg/(2K)\) é o alongamento estático, e puxa com \(K f^{+}\), com \(f^+=\max(f,0)\). Newton para o centro dá \(m\ddot y=-K(f_L^++f_R^+-2y_0)\) mais amortecimento e força externa. Para a rotação, o torque de uma força vertical \(F\) à distância \(\ell\cos\theta\) do centro é \(F\ell\cos\theta\), e o momento de inércia de uma barra em torno do centro é \(\int_{-\ell}^{\ell}(m/2\ell)\,x^2\,dx=m\ell^2/3\). Isso dá a segunda equação do cartão-resumo.

Então \(f^+=f\), e para \(\theta\) pequeno as equações viram \(\ddot y+\delta\dot y+\tfrac{2K}{m}y=\varepsilon\sin\mu t\) e \(\ddot\theta+\delta\dot\theta+\tfrac{6K}{m}\theta=0\): dois osciladores independentes. A força vertical excita só \(y\); a torção \(\theta\) não tem força externa e qualquer torção inicial simplesmente se extingue. Com \(K=8000\,\mathrm{N/m^2} e\) \(m=2500\,\mathrm{kg/m}\), \(\omega_y=\sqrt{2K/m}=2{,}53\) rad/s (período \(2{,}48\,\mathrm{s}\)) e \(\omega_\theta=\sqrt3\,\omega_y=4{,}38\) rad/s (período \(1{,}43\,\mathrm{s}\)). Rode “Mesma rajada, pendurais sempre esticados”: a torção fica em cerca de 0,06°, a semente que lhe demos. A razão \(\sqrt3\) não é coincidência dos números: vem só da geometria, inércia \(m\ell^2/3\) contra duas molas à distância \(\ell\).

Em repouso os dois pendurais sustentam o peso: \(2Ky_0=mg\), logo \(y_0=mg/(2K)=2500\cdot9{,}8/16000=1{,}53\) m. Esse número importa para o afrouxamento: um pendural afrouxa quando seu alongamento \(y_0+y\mp\ell\sin\theta\) fica negativo, isto é, quando o centro do tabuleiro sobe mais que \(y_0\) acima da posição de repouso (com \(\theta=0\)), ou antes, quando o tabuleiro está inclinado. Um pendural mais rígido (\(K\) maior) dá \(y_0\) menor, então o afrouxamento ocorre com movimento menor, mas \(\omega_y\) também sobe. Tentativas anteriores com \(K=1000\) deram um afundamento irreal de \(12\,\mathrm{m}\) e torções de dezenas de radianos, por isso o padrão é \(K=8000\).

No modelo esticado o torque é \(-2K\ell^2\sin\theta\cos\theta\) e não depende de \(y\). Com afrouxamento ele passa a ser \(K\ell\cos\theta\,(f_L^+-f_R^+)\). Se o tabuleiro sobe o bastante, os dois pendurais afrouxam e o torque some, liberando a torção da sua mola restauradora. Quando um lado afrouxa antes do outro, o que acontece assim que \(\theta\neq0\), o torque não é simétrico e depende do deslocamento vertical. Assim o movimento vertical muda a rigidez do modo de torção a cada ciclo, e uma rigidez que varia no ritmo certo pode bombear energia numa oscilação (excitação paramétrica, como uma criança que bombeia o balanço ficando de pé e agachando). Essa é uma leitura plausível, não uma prova: os experimentos mostram o que acontece; a matemática de quando é difícil.

Rode “Rajada na frequência vertical, pendurais frouxos” e observe a curva da elevação da borda: uma torção minúscula de 0,06° é amplificada até o limiar de 5,7° e além.

Com RK4 e \(h=0{,}02\) s a varredura dá, por exemplo, \(\varepsilon_{crit}\approx6{,}6\,\mathrm{m/s^2}\) em \(\mu/\omega_y=0{,}6\), 3,2 em 0,8, 0,78 em 1,0, 1,4 em 1,04, 7,4 em 1,08 e nenhuma torção até \(8\,\mathrm{m/s^2}\) de 1,16 em diante. Três lições. Primeira, há um vale perto da frequência vertical, como a ressonância sugere. Segunda, ele é assimétrico: basta uma força menor do lado das baixas frequências que do lado das altas, o que combina com a imagem de uma mola que amolece quando os pendurais afrouxam (uma mola que amolece entra em ressonância em frequência menor quando excitada com força); é uma leitura plausível, não algo que o aplicativo prove. Terceira, a curva é irregular: com acoplamento não linear a resposta não é monótona em \(\varepsilon\), então “o menor \(\varepsilon\) que torce” depende de como se procura. Varremos para cima numa grade geométrica, delimitamos o primeiro cruzamento e fazemos bisseção dentro dele; um limiar escondido entre dois pontos da varredura pode passar despercebido.

O método de Euler transforma \(y''=-\omega^2y\) numa aplicação cuja amplitude cresce a cada passo por \(\sqrt{1+z^2}\), \(z=\omega h\): uma oscilação sem amortecimento espirala para fora qualquer que seja o \(h\). Para o modo de torção com \(h=0{,}02\) s, \(z=0{,}088\) e a amplitude cresce 0,4% por passo, cerca de 21% por segundo simulado, e a energia cerca de 47% por segundo no início. Experimente “Queda livre, sem vento: Euler”: a curva de energia sobe e a linha de status avisa que a energia cresceu sem vento. A mesma trajetória com RK4 e o mesmo \(h\) fica \(a 3\cdot 10^{-4}\) da energia inicial em \(30\,\mathrm{s}\), e Heun \(a 5\cdot 10^{-3}\) (deriva relativa máxima em \(30\,\mathrm{s}\): Euler 4,9, Heun \(4{,}7\cdot 10^{-3}\), RK4 \(2{,}9\cdot 10^{-4}\) com \(h=0{,}02\); com \(h=0{,}005\): 1,1, \(1{,}3\cdot 10^{-4}\), \(1{,}8\cdot 10^{-5}\).

A consequência para o vento crítico: com Euler a varredura devolve limiares próximos de zero em todo lugar, porque o próprio método excita a torção.

Euler tem erro local \(O(h^2)\) e global \(O(h)\); Heun (a regra do trapézio explícita) global \(O(h^2)\); o método clássico de Runge–Kutta global \(O(h^4)\). No cartão do passo, numa trajetória suave (pendurais sempre esticados), as inclinações medidas são 0,95, 2,02 e 4,02. Quando os pendurais ligam e desligam, o lado direito tem quinas, então a expansão de Taylor por trás das provas de ordem falha perto de cada troca: Euler e Heun mantêm aproximadamente sua ordem (0,99 e 2,04), mas o RK4 cai para cerca de 2,4 e sua curva de erro é irregular, dependendo de onde as quinas caem entre os passos. Moral: um método de ordem alta não dá alta precisão num problema não suave; para recuperá-la seria preciso localizar cada instante de troca e reiniciar o integrador ali (tratamento de eventos).

Para \(y''=-\omega^2y\) um passo de RK4 multiplica o estado por \(R(z)=1-\tfrac{z^2}{2}+\tfrac{z^4}{24}+i\bigl(z-\tfrac{z^3}{6}\bigr)\), \(z=\omega h\), e \(|R|\le1\) apenas para \(z\le2\sqrt2\approx2{,}83\). Com \(h=0{,}7\) s o modo de torção tem \(z=4{,}38\cdot0{,}7=3{,}07\), logo além do limite, e \(|R|=1{,}75\): a oscilação cresce 75% por passo. Rode “RK4 com um passo enorme”: a energia sobe por um fator de várias centenas, embora haja amortecimento e não haja vento. No modelo não linear a torção não cresce sem limite (o torque \(\sin\theta\cos\theta\) satura), então você vê um movimento errático de grande amplitude em vez de um estouro: a curva de energia é o alarme confiável. Abaixo do limite o RK4 é estável e amortece levemente as altas frequências (em \(z=2\) a amplitude perde cerca de 25% por passo), então um passo grande demais também pode fazer uma oscilação real morrer só por amortecimento numérico.

Nunca se sabe com certeza; coleta-se evidência. (1) Verifique um caso com resposta conhecida: o modelo esticado tem períodos \(2\pi/\omega_y=2{,}484\) s e \(2\pi/\omega_\theta=1{,}434\) s, e a simulação os reproduz com quatro dígitos. (2) Verifique um invariante: sem vento e sem amortecimento a energia deve ser constante. (3) Divida o passo ao meio e veja se a resposta muda por cerca de \(2^{-p}\); para RK4 em trechos suaves, um fator 16. (4) Compare dois métodos: Heun e RK4 dão os mesmos limiares críticos com três dígitos (0,776 contra 0,78 em \(\mu/\omega_y=1\)), enquanto Euler não, o que diz em qual desconfiar. Concordância não é prova, mas discordância é um aviso certo.

Para ensino: com passo fixo, as fórmulas do fator de amplificação e da ordem se aplicam diretamente, e dá para ver como o erro depende de \(h\). Um código de produção usaria um par embutido adaptativo, como Dormand–Prince, para escolher \(h\) a partir de uma estimativa de erro, e um localizador de eventos para pousar exatamente em cada troca frouxo/esticado. Passos adaptativos encolhem sozinhos perto das trocas, que é onde o RK4 de passo fixo perde sua ordem.

Comece com “Rajada na frequência vertical”. Mantenha a intensidade em 1,2 e mova o controle de frequência de 0,90 a 1,10 em passos de 0,02. Em nossas execuções o tabuleiro só torce para \(0{,}94\lesssim\mu/\omega_y\lesssim1{,}04\), chegando a 12° a 22° em dois minutos, e leva cerca de \(25\,\mathrm{s}\) para cruzar o limiar de 5,7° em 1,00. Fora dessa faixa a torção fica na semente. Depois baixe a intensidade para 0,5: nenhuma frequência torce o tabuleiro, e o menor limiar de toda a varredura é \(0{,}78\,\mathrm{m/s^2}\). Compare com a curva do vento crítico (rode a varredura antes): a faixa em que ela fica abaixo de 1,2 é a faixa que você achou.

Mude \(K\) de 8000 para 2000 e depois para 30000, com a frequência da rajada fixa em 1,0 vez \(\omega_y\). Preveja primeiro: como \(\omega_y\), \(\omega_\theta\) e \(y_0\) variam com \(K\)? (Respostas: \(\omega\propto\sqrt K\), caindo pela metade, para 1,26 e \(2{,}53\,\mathrm{rad/s}\), quando \(K\) é dividido por 4, e \(y_0\propto1/K\), crescendo para \(6{,}1\,\mathrm{m}\).) Depois confira a leitura ao lado do controle. Qual ajuste torce com mais facilidade? Note que o menu de experimentos vai para “Personalizado” e a simulação reinicia, porque mudar \(K\) muda o que significa “repouso”.

Escreva \(y''=-\omega^2y\) como sistema de primeira ordem \(u'=Au\), \(u=(y,y')\), e aplique um passo de cada método. Como \(A\) tem autovalores \(\pm i\omega\), basta aplicar o método à equação escalar de teste \(u'=\lambda u\), \(\lambda=i\omega\). Euler: \(R=1+\lambda h=1+iz\), logo \(|R|=\sqrt{1+z^2}\). Heun: \(R=1+\lambda h+\tfrac12(\lambda h)^2=1-\tfrac{z^2}2+iz\), e \(|R|^2=(1-\tfrac{z^2}{2})^2+z^2=1+\tfrac{z^4}{4}\). RK4: \(R\) é o polinômio de Taylor de \(e^{\lambda h}\) até a ordem 4: \(1+\lambda h+\tfrac{(\lambda h)^2}{2}+\tfrac{(\lambda h)^3}{6}+\tfrac{(\lambda h)^4}{24}\). Mostre que \(|R|=1\) é atingido em \(z=2\sqrt2\) (ali \(R=-\tfrac13-\tfrac{2\sqrt2}{3}i\)) e que perto de \(z=0\), \(|R|^2=1-z^6/72+\ldots\), de modo que o RK4 amortece levemente. Confira com o gráfico.

Liste o que o modelo deixa de fora e diga o que mais importaria. Uma lista parcial: a força do vento depende do ângulo de ataque e, portanto, de \(\theta\) e \(\dot\theta\) (isto é o flutter); os cabos e as torres têm movimento próprio; a ponte é uma estrutura longa com muitos modos, não uma seção; os pendurais têm comprimentos e tensões diferentes; o amortecimento não é uma única constante viscosa. Depois decida qual termo novo tornaria o amortecimento \(\delta\) negativo acima de uma velocidade crítica do vento, que é a assinatura do flutter. (Dica: uma força proporcional a \(+\dot\theta\) na equação do torque.)

Andar empurra o tabuleiro de lado com uma força que alterna esquerda e direita uma vez a cada dois passos, a cerca de \(1\,\mathrm{Hz}\), perto das frequências laterais de muitas passarelas. Se os empurrões de pessoas diferentes são aleatórios, eles em grande parte se cancelam e a ponte sente uma força da ordem de \(\sqrt N\). Se as pessoas estão em passo entre si e com o balanço, as forças somam-se até a ordem \(N\) e empurram a ponte em fase com a sua própria velocidade, que é exatamente o que significa amortecimento negativo.

É uma malha de realimentação. Mais balanço faz mais pedestres ajustarem o passo, mais pedestres em passo fazem mais balanço. Se a malha amplifica ou se extingue depende de o empurrão da multidão superar o amortecimento da ponte. Abaixo do limiar qualquer balanço pequeno é amortecido e os pedestres ficam nos seus passos; acima dele a malha dispara até saturar, porque pedestres de passo natural distante do da ponte não conseguem travar. Teste os experimentos “abaixo do crítico” (um empurrão se extingue) e “acima do crítico” (uma semente minúscula cresce em dois minutos).

A multidão crítica é proporcional ao coeficiente de amortecimento: \(N_c=4\zeta K\gamma/(GC\sin\alpha)\propto B\). Dobrar o amortecimento dobra o limiar. A rigidez entra só por \(K\); triplicar o limiar enrijecendo exigiria três vezes a rigidez, movendo a frequência por um fator \(\sqrt3\), e uma ponte mais rígida continua a ser excitada por uma multidão na própria frequência. Acrescentar amortecimento eleva o limiar diretamente e a baixo custo. No aplicativo, o botão “Instalar amortecedores” leva \(\zeta\) a 3 %, de modo que \(N_c\) vai de 143 a 428 e uma multidão de 250 que balançava tem o balanço levado abaixo de 1 cm em cerca de \(35\,\mathrm{s}\).

Um modo lateral da ponte é um oscilador \(M\ddot X+B\dot X+KX=G\sum\sin\Theta_i\), em que \(\Theta_i\) é a fase do ciclo de balanço do \(i\)-ésimo pedestre e \(G\) a amplitude da força lateral de um pedestre. Cada pedestre obedece a \(\dot\Theta_i=\omega_i+CA\sin(\Psi-\Theta_i+\alpha)\): seu ritmo próprio \(\omega_i\) mais uma atração para a fase \(\Psi\) da ponte, proporcional à amplitude \(A\) da ponte. Com \(X=A\sin\Psi\), \(\dot X=\Omega_0A\cos\Psi\), a amplitude e a fase são calculadas de \(X\) e \(\dot X\) a cada passo. Com \(\alpha=\pi/2\) a força dos pedestres fica um quarto de volta à frente do deslocamento e assim empurra em fase com a velocidade. As \(N+2\) equações são integradas juntas. É um modelo de Kuramoto em que o acoplamento não é fixo, mas carregado pela própria ponte.

Procure um estado estacionário com amplitude de balanço \(A\). No referencial que gira com a ponte, um pedestre com dessintonia \(\nu=\omega_i-\Omega_0\) obedece a \(\dot\varphi=\nu-CA\sin\varphi\). Os de \(|\nu|\lt CA\) travam em \(\sin\varphi=\nu/CA\); os demais derivam e, em média, não contribuem para a força. Os travados dão o parâmetro de ordem \(r=\int_{-CA}^{CA}P(\nu)\sqrt{1-\nu^2/(CA)^2}\,d\nu\), e a ponte, excitada na ressonância, tem seu amortecimento equilibrado pela parte da força da multidão em fase com a sua velocidade: \(B\Omega_0A=GNr\sin\alpha\). Para \(A\) pequeno, \(r\approx\tfrac\pi2P(0)\,CA\); substituindo, só há \(A\) não nulo se \(N>N_c=\dfrac{2B\Omega_0}{\pi\,G\,C\,P(0)\sin\alpha}\). Com uma lorentziana, \(P(0)=1/(\pi\gamma)\) e \(N_c=2\gamma B\Omega_0/(GC\sin\alpha)\). Equivalentemente \(N_c=4\zeta K\gamma/(GC\sin\alpha)\), proporcional à razão de amortecimento \(\zeta\).

Para uma lorentziana a integral sai em forma fechada: \(r=(\sqrt{a^2+\gamma^2}-\gamma)/a\) com \(a=CA\). O balanço \(B\Omega_0A=GNr\sin\alpha\) vira então \(\sqrt{a^2+\gamma^2}-\gamma=q\,a^2\), cuja solução é \(a^2=(1-2\gamma q)/q^2\); daí \(A_\ast=\dfrac{G\sin\alpha}{B\Omega_0}\sqrt{N(N-N_c)}\) e \(r_\ast=\sqrt{1-N_c/N}\). Uma lorentziana verdadeira tem caudas que chegam a frequências de caminhada negativas; por isso o aplicativo a corta em \(\pm T=\pm\min(5\gamma,0{,}8\Omega_0)\) e renormaliza; a fração mantida é \(Z=\tfrac2\pi\arctan(T/\gamma)=0{,}874\) nos ajustes iniciais. Tudo acima vale então com \(N_c\) multiplicado por \(Z\), e \(A_\ast\) e \(r_\ast\) divididos por \(Z\), enquanto \(CA_\ast\le T\).

Com os ajustes iniciais a multidão crítica prevista é \(N_c=142{,}6\). Simulando \(1500\,\mathrm{s}\) e tomando a média dos últimos \(150\,\mathrm{s}\), para \(N=150,175,200,250,300,400\) os balanços são 2,15, 5,15, 7,40, 11,36, 15,09 e 22,31 cm contra os previstos 2,33, 5,28, 7,51, 11,48, 15,23 e 22,49 cm, com diferença de 1 a 8 %; a discrepância é maior perto do limiar, onde o crescimento é lento. O parâmetro de ordem concorda a menos de 0,005 para \(N\ge200\). Ajustando \(A^2/N\) contra \(N\) aos balanços simulados de \(N=175\) a \(300\) e lendo o intercepto, obtém-se \(N_c=143{,}7\), a 1 % da fórmula. Os pedestres são amostrados por quantis, então não há aleatoriedade: as simulações são exatamente reproduzíveis.

O estado tem \(N+2\) números: deslocamento e velocidade da ponte e \(N\) fases. A parte rápida é a oscilação da ponte em \(\Omega_0\) e das fases dos pedestres nos seus ritmos, tudo dentro de alguns rad/s, e a parte lenta é o crescimento ou decaimento da amplitude ao longo de minutos. Um método de quarta ordem permite que um passo cubra boa fração de um período, de modo que 25 minutos de tempo simulado custam menos de 20 000 passos. O custo por passo é proporcional a \(N\) (quatro avaliações, cada uma um seno por pedestre), então uma simulação com 250 pedestres leva menos de um segundo.

Dois limites. Para a ponte sozinha, o fator de amplificação por passo de cada método é função de \(z=\Omega_0h\): \(|R|=\sqrt{1+z^2}\) para Euler, \(\sqrt{1+z^4/4}\) para Heun, ambos maiores que 1 para todo \(z>0\), de modo que injetam energia à taxa \(\Omega_0^2h/2\) (Euler) e \(\Omega_0^4h^3/8\) (Heun) por segundo. O amortecimento verdadeiro é \(\zeta\Omega_0\), então um método fabrica um balanço quando sua injeção o supera: \(h>2\zeta/\Omega_0=0{,}004\) s para Euler e \(h>(8\zeta/\Omega_0^3)^{1/3}=0{,}086\) s para Heun. Isso concorda com a tabela acima (o Euler já cresce em \(0{,}004\,\mathrm{s}\), o Heun em \(0{,}09\,\mathrm{s}\)). Para o sistema acoplado com RK4, o erro no balanço final em relação a \(h=0{,}01\) é \(-76\%\) em \(h=0{,}2\), \(-3{,}4\%\) em 0,1, \(-0{,}12\%\) em 0,05 e abaixo de 0,005 % em 0,02.

Frequências aleatórias tornariam cada simulação diferente, e um limiar medido em algumas centenas de pedestres carregaria ruído de vários por cento. Colocar o pedestre \(i\) no quantil \((i+\tfrac12)/N\) da distribuição faz da soma sobre os pedestres uma regra do ponto médio para a integral que define \(r\), de modo que a simulação é determinística e o erro de \(N\) finito é suave. As fases iniciais são uma sequência da razão áurea \(2\pi\{i\varphi\}\), que se espalha uniformemente (o parâmetro de ordem inicial é pequeno), e as duas sequências são embaralhadas para que os primeiros \(n\) pedestres, os que já entraram quando a multidão chega aos poucos, ainda representem toda a distribuição.

Escolha o experimento “abaixo do crítico”, aperte “Empurrar a ponte” algumas vezes e veja o balanço morrer. Aumente \(N\) de 10 em 10 e repita. Em que \(N\) um empurrão de 1 cm deixa de morrer e passa a crescer? Como o crescimento perto do limiar é lento, espere alguns minutos de tempo simulado (velocidade \(32\times ) e\) compare com o \(N_c\) previsto na linha de métricas. Depois dobre \(\zeta\) e veja como o limiar se move.

Rode “Pedestres com passos quase iguais”. A multidão crítica cai de 143 para 54, porque \(N_c\) é proporcional à dispersão \(\gamma\): quando todos caminham naturalmente quase no mesmo ritmo, mais gente está ao alcance da frequência da ponte e uma multidão menor basta. Com os ajustes desse experimento e \(N=100\), o balanço previsto é 4,8 cm e o simulado 4,7 cm. Agora baixe \(C\) para 8 e confira que \(N_c\) dobra: é inversamente proporcional à sensibilidade.

Rode “Sem pedestres, Euler explícito”. A ponte está vazia e amortecida, mas a amplitude cresce. Baixe \(h\) para \(0{,}004\,\mathrm{s}\) e veja o crescimento parar; baixe mais e veja o decaimento que deveria estar ali. Troque para Heun com \(h=0{,}1\): também crescimento. Depois troque para RK4. O que um modelo da ponte de Tacoma, uma ponte real cujo comportamento dependia de um amortecimento pequeno, teria dito se o autor tivesse usado Euler explícito?

Na relatividade geral a gravidade é a curvatura do espaço-tempo, e a luz, que não tem massa para ser empurrada, simplesmente segue os caminhos mais retos possíveis nesse espaço-tempo curvo, as geodésicas nulas. Perto de uma massa esses caminhos deixam de ser retas no espaço que desenhamos. Um argumento newtoniano que trata a luz como partículas de velocidade \(c\) também prevê um desvio, mas pela metade: \(2M/b\) contra o \(4M/b\) de Einstein. A outra metade vem da curvatura do próprio espaço, e não só da parte temporal da métrica.

O horizonte \(r=2M\) é o ponto sem volta: nada que esteja dentro consegue sair, nem a luz. Em \(r=3M\), a esfera de fótons, a luz pode viajar em círculo em torno do buraco, mas a órbita é instável como uma bola equilibrada no topo de uma colina: o menor empurrão manda a luz para dentro ou para fora. Um raio que chega de longe é capturado se o seu parâmetro de impacto for menor que \(b_c=3\sqrt3\,M\approx5{,}196\,M\), que não é \(3M\) por causa da curvatura do raio na chegada. O disco escuro que um observador distante enxerga, a sombra, tem raio \(b_c\).

Imagine o raio vindo de muito longe ao longo de uma reta. O parâmetro de impacto \(b\) é a distância entre essa reta e a paralela que passa pelo buraco: a altura do raio na vista de cima. É a razão \(L/E\) de duas grandezas conservadas, o momento angular e a energia, então é o mesmo em todos os pontos do raio e sozinho decide a órbita inteira. Por isso o app só precisa de um número para lançar um raio: \(b=8\) dá sempre a mesma curva.

Escreva as grandezas conservadas \(E=(1-2M/r)\dot t\) e \(L=r^2\dot\varphi\) (pontos são derivadas no parâmetro afim) e a condição nula \(0=-(1-2M/r)\dot t^2+\dot r^2/(1-2M/r)+r^2\dot\varphi^2\). Eliminando \(t\) obtém-se \(\dot r^2=E^2-(L^2/r^2)(1-2M/r)\). Com \(u=1/r\) e \(du/d\varphi=-\dot r/L\) isto vira \((du/d\varphi)^2=1/b^2-u^2+2Mu^3\), a integral primeira. Derivando uma vez em \(\varphi\) e dividindo por \(2u'\) vem \(u''+u=3Mu^2\).

Só as duas últimas linhas são necessárias neste curso; a dedução a partir da métrica é um bom exercício para uma aula de relatividade.

Ponha \(y_1=u\) e \(y_2=u'\). Então \(y_1'=y_2\) e \(y_2'=3y_1^2-y_1\) (com \(M=1\)): um sistema de primeira ordem \(\mathbf y'=\mathbf f(\mathbf y)\) em que o ângulo \(\psi\) não aparece à direita (um sistema autônomo). Euler: \(\mathbf y_{n+1}=\mathbf y_n+h\,\mathbf f(\mathbf y_n)\). Heun: média da inclinação no início e na previsão de Euler. RK4: quatro inclinações com pesos \(1,2,2,1\) sobre 6. O arquivo bh-core.js tem os três em umas vinte linhas, e o script Python no fim da página também.

O horizonte está em \(u=1/2\). Depois de cada passo o programa verifica se \(u\ge1/2\); se sim, o raio é capturado, o sino toca e o último passo é encurtado para que o caminho termine exatamente no horizonte. O raio escapa quando \(u'<0\) e \(u\) volta ao valor inicial \(u_0=1/r_0\). Se nenhuma das duas coisas acontece depois de \(14\pi\) radianos, o raio é declarado indeciso, ainda circulando. Localizar o último passo é um problema de busca de raiz: algumas iterações de Newton sobre a fração do passo, feitas com o mesmo método, de modo que o método mantém sua ordem.

Sem gravidade o raio varreria o ângulo polar de \(\psi_0=\arcsin(b/r_0)\) a \(\pi-\psi_0\). O raio curvado varre um ângulo maior \(\psi_{\text{saída}}-\psi_0\), e a diferença é o desvio \(\alpha=(\psi_{\text{saída}}-\psi_0)-\pi+2\psi_0\). O programa começa em \(r_0=1000M\), longe o bastante para que o resultado fique próximo do valor de um raio vindo do infinito. Em primeira ordem \(\alpha=4M/b\); a expansão continua como \(4M/b+15\pi M^2/(4b^2)+128M^3/(3b^3)+\dots\), e falha perto de \(b_c\), onde \(\alpha\) diverge.

Linearize a equação da órbita em torno da órbita circular de fótons \(u=1/3\): com \(u=1/3+\delta\) obtém-se \(\delta''=\delta+3\delta^2\), então pequenos erros crescem como \(e^{\psi}\) (uma sela, não um centro). Um raio com \(b-b_c\) igual a \(10^{-6}\) gasta cerca de \(\ln10^{6}\approx14\) radianos, mais de duas voltas, perto da esfera de fótons antes de se decidir. Um erro numérico desse tamanho é amplificado do mesmo modo, então localizar \(b_c\) com muitos dígitos exige um método cujo erro esteja bem abaixo dos dígitos desejados. O RK4 com \(h=0{,}02\) tem isso; o Euler não.

A equação é um oscilador não linear, \(u''=-u+3u^2\). O método de Euler aplicado a um oscilador multiplica a amplitude por \(\sqrt{1+h^2}\) a cada passo: a solução numérica ganha energia. Aqui a energia é a integral primeira \(1/b^2\), e ganhá-la significa se comportar como um raio de \(b\) menor: o Euler puxa para dentro raios que deveriam passar. O fator de amplificação do Heun é \(\sqrt{1+h^4/4}\), menor por duas potências de \(h\), e o do RK4 é cerca de \(1-h^6/144\) (perde um pouco), então esses dois são muito melhores. O mesmo estudo de estabilidade é feito no app Tacoma.

No tempo coordenado de Schwarzschild um raio que se aproxima do horizonte desacelera e nunca chega: \(dr/dt\to0\). Um método que marcha em \(t\) se arrastaria e nunca veria a captura. No ângulo \(\psi\) a equação \(u''=3u^2-u\) é suave em todo lugar e o horizonte é simplesmente o valor \(u=1/2\). Para o filme o app usa o parâmetro afim, \(d\psi/d\lambda=bu^2\), que longe é quase o tempo comum e mantém o movimento vivo perto do buraco.

Quando \(b\to0\) o raio entra quase em linha reta e o ângulo \(\psi\) quase não muda, então um passo fixo em \(\psi\) cobre quase toda a queda: o controle para em \(|b|=1\) por isso. Para raios radiais usa-se outra formulação, por exemplo \(r(\lambda)\) com \((dr/d\lambda)^2=1-b^2(1-2M/r)/r^2\), ou um passo adaptativo. Escolher bem a variável independente faz parte da modelagem.

Escolha “Euler com passo grande mente” (\(b=5{,}3\), \(h=0{,}1\)). A resposta exata é que o raio escapa, desviado de 204°, mas o Euler o faz cair, porque o erro da integral primeira empurra o valor crítico efetivo para \(b_c\approx5{,}356\). Aperte “Rodar a bisseção” para ver esse número. Depois diminua \(h\) no controle: abaixo de cerca de \(h=0{,}07\) o raio volta a escapar, e a estimativa de \(b_c\) tende a \(5{,}196\) por cima, aproximadamente em proporção a \(h\). Faça o mesmo com Heun e com RK4 em \(h=0{,}1\): os dois acertam o destino.

Lance \(b=6\), \(5{,}3\) e \(5{,}2\) em sequência com o menu de experimentos e leia o desvio na linha de status: 99°, 204° e 390°. Em \(b=5{,}2\), apenas \(0{,}004\) acima de \(b_c\), o raio dá mais de uma volta em torno do buraco antes de sair. Pela linearização da seção anterior cada volta completa a mais exige que o raio esteja \(e^{2\pi}\approx535\) vezes mais perto de \(b_c\), então o passo de \(0{,}05\) do controle é grosso demais para uma segunda volta; use a bisseção do script Python. A lei logarítmica \(\alpha\approx-\ln(b-b_c)+\text{const}\) se ajusta à curva verde perto de \(b_c\)?

1. Mostre pela integral primeira que a captura ocorre exatamente para \(b<3\sqrt3\,M\): o raio não pode voltar quando \(F(u)=1/b^2-u^2+2u^3\) permanece positiva em \([0,1/2]\); o mínimo de \(F\) ali está em \(u=1/3\), onde \(F=1/b^2-1/27\).
2. Escreva Euler, Heun e RK4 para \(u''+u=3u^2\) na sua linguagem e reproduza a tabela da fronteira de captura.
3. Use \(\alpha\approx4M/b\) para estimar o desvio da luz na borda do Sol (\(M=1{,}48\,\mathrm{km}\), \(b=6{,}96\times10^{5}\,\mathrm{km}\)) em segundos de arco e compare com \(1{,}75''\).
4. Tome \(r\) como variável independente e o parâmetro afim como tempo, e tente integrar um raio com \(b\) logo acima de \(b_c\). O que acontece no ponto de retorno, onde \(dr/d\lambda=0\)?
5. Com a linearização em torno de \(u=1/3\), estime quantas voltas em torno do buraco dá um raio com \(b-b_c=10^{-8}M\).

A equação \(y^2=k\) é fácil: \(y=\pm\sqrt k\). O problema de \(x^2+px=q\) é que \(x\) aparece duas vezes, uma ao quadrado e outra sozinha. Completar o quadrado reescreve o lado esquerdo como um só quadrado, \((x+p/2)^2\), de modo que \(x\) aparece uma vez e a raiz quadrada pode descascá-lo. A figura geométrica mostra por que funciona: um quadrado mais uma faixa de mesma altura é quase um quadrado maior, e o canto que falta tem área \((p/2)^2\).

O primeiro passo de verdade divide por \(a\). Se \(a=0\) a equação \(bx+c=0\) é linear, com a única solução \(x=-c/b\), e a fórmula dividiria por zero. Há um limite agradável: quando \(a\to0\) com \(b>0\), a raiz \(\frac{-b+\sqrt{b^2-4ac}}{2a}\) tende a \(-c/b\) enquanto a outra foge para o infinito. A forma estável \(x=c/q\) mostra esse limite de imediato, o que a forma ingênua esconde atrás de \(0/0\). Teste o exemplo com \(a\) minúsculo.

Os dois sinais dão uma solução: se \(y^2=k\) então \((-y)^2=k\) também. Assim a fórmula é uma maneira curta de escrever dois números, \(x_1=\frac{-b+\sqrt\Delta}{2a}\) e \(x_2=\frac{-b-\sqrt\Delta}{2a}\), que são imagens espelhadas em torno do vértice \(x=-b/2a\). Se \(\Delta=0\) eles coincidem (raiz dupla, a parábola só toca o eixo); se \(\Delta<0\) são complexos conjugados e a parábola não alcança o eixo. Na animação o sinal aparece no passo da raiz quadrada e sobrevive até o fim.

Dividir por \(a\ne0\), somar o mesmo número aos dois lados e reescrever uma expressão como outra igual podem ser desfeitos, então toda equação da cadeia tem exatamente as mesmas soluções que a anterior. Só o passo da raiz quadrada exige cuidado: de \(y^2=k\) para \(y=\pm\sqrt k\) ele mantém as duas soluções. Uma verificação concreta: para qualquer \(x\), a diferença entre os dois lados de cada um dos sete primeiros quadros vale \(f(x)/a\) com \(f(x)=ax^2+bx+c\) (ou o próprio \(f(x)\) no primeiro quadro), então os sete se anulam exatamente nos mesmos \(x\). O script em Python abaixo confere isso com frações exatas.

\(\Delta=b^2-4ac\) é o número sob a raiz, e a raiz é real exatamente quando \(\Delta\ge0\). Geometricamente, a parábola \(y=ax^2+bx+c\) tem vértice em \(x=-b/2a\) com altura \(-\Delta/4a\), então cruza o eixo duas vezes quando \(a\) e \(-\Delta\) têm o mesmo sinal (\(\Delta>0\)), toca uma vez quando \(\Delta=0\) e não o alcança quando \(\Delta<0\). A distância entre as duas raízes é \(\sqrt\Delta/|a|\). Para coeficientes inteiros \(\Delta\) é inteiro, e as raízes são racionais exatamente quando \(\Delta\) é um quadrado perfeito.

Some e multiplique as duas raízes, com \(s=\sqrt\Delta\): \(x_1+x_2=\frac{-b+s}{2a}+\frac{-b-s}{2a}=-\frac ba\) e \(x_1x_2=\frac{(-b)^2-s^2}{4a^2}=\frac{b^2-(b^2-4ac)}{4a^2}=\frac ca\). São as relações de Viète, e elas dizem que \(a(x-x_1)(x-x_2)=ax^2+bx+c\), o que é o mesmo que dizer que \(x_1\) e \(x_2\) são as raízes. A relação do produto também é a chave da forma estável da próxima seção.

Para \(a=1\), \(p>0\), \(q>0\) é: um quadrado de lado \(x\) e uma faixa \(p\times x\) têm área \(x^2+px=q\); cortar a faixa em duas metades \(\frac p2\times x\) e colá-las em dois lados do quadrado dá um L de área \(q\); somar o canto \((p/2)^2\) faz um quadrado de lado \(x+p/2\) e área \(q+(p/2)^2\), de modo que \(x+\frac p2=\sqrt{q+p^2/4}\). É exatamente a fórmula da raiz positiva, \(x=-\frac p2+\sqrt{\frac{p^2}{4}+q}\). Para outros sinais as peças teriam comprimentos negativos e a figura deixa de ser figura, enquanto a álgebra não se importa.

Quando \(b^2\gg|4ac|\), \(\sqrt{b^2-4ac}\) é quase \(|b|\), e para um dos dois sinais a fórmula calcula \(-b+\sqrt\Delta\) (para \(b>0\)) como diferença de dois números quase iguais. Cada um carrega um erro de arredondamento de cerca de \(\varepsilon|b|\), com \(\varepsilon=2^{-53}\approx1{,}1\times10^{-16}\), e a subtração cancela os dígitos principais e deixa esse erro na frente: o erro absoluto do resultado continua cerca de \(\varepsilon|b|/(2|a|)\), enquanto a raiz mesma vale só \(|c/b|\). O erro relativo é da ordem de \(\varepsilon\,b^2/|ac|\). Para \(x^2+10^8x+1\): a raiz pequena é \(-10^{-8}\), a fórmula em precisão dupla devolve \(-7{,}45\times10^{-9}\) (erro de 25 %), e para \(b\ge2^{27}\approx1{,}34\times10^8\) devolve exatamente 0. A outra raiz, \(-b-\sqrt\Delta\), não tem cancelamento e está certa.

Calcule primeiro a raiz que não tem cancelamento, escolhendo o sinal que soma módulos: \(q=-\tfrac12\bigl(b+\operatorname{sgn}(b)\sqrt\Delta\bigr)\), de modo que \(q\) tem o tamanho de \(b\) e nenhum dígito se perde, e depois \(x_1=q/a\). A outra raiz vem do produto de Viète \(x_1x_2=c/a\): \(x_2=\frac{c/a}{q/a}=c/q\). Ambas são divisões e somas de números de mesmo sinal, então cada uma tem erro relativo de poucas unidades de \(\varepsilon\). Na tabela acima a coluna estável está correta até o último dígito impresso em todas as linhas. É o que as bibliotecas fazem: as rotinas-padrão para quadráticas usam essa forma, não a do livro-texto.

Primeiro, estouro: \(b^2\) passa do maior double (cerca de \(1{,}8\times10^{308}\)) para \(|b|>1{,}3\times10^{154}\) mesmo que as raízes sejam modestas; multiplicar todos os coeficientes por uma potência de dois resolve sem mudar as raízes. Segundo, um discriminante quase nulo: \(b^2\) e \(4ac\) quase iguais, e \(\Delta\) perde dígitos; perto de uma raiz dupla as raízes são mal condicionadas, uma perturbação \(\delta\) dos coeficientes as move cerca de \(\sqrt\delta\), então \(10^{-16}\) vira \(10^{-8}\) não importa como se escreva a fórmula. Um cálculo compensado de \(b^2-4ac\) (com multiplicação-soma fundida) recupera os dígitos perdidos quando os coeficientes são exatos. Terceiro, \(a\) minúsculo: com \(a=10^{-9}\) no exemplo, a fórmula para a raiz que tende a \(-c/b\) perde dígitos, e a forma estável não.

Ponha \(p=4\) e \(q=5\) na figura de áreas. A raiz positiva é \(x=-2+\sqrt{4+5}=1\). Confira com a figura: o quadrado tem lado 1, as duas faixas são \(2\times1\), o canto é \(2^2=4\), o quadrado grande tem lado \(x+2=3\) e área \(9=q+4\). Agora tome \(p=3\), \(q=10\) e calcule \(x\) e o quadrado grande à mão antes de olhar a leitura. Por que a leitura nunca pede \(q<0\)?

No laboratório, digite \(a=1\), \(c=1\) e \(b=10^6\), depois \(10^7\), \(10^8\) e \(1{,}5\times10^8\). Compare a fórmula como escrita com a raiz exata e note onde o erro relativo passa de 1 % e onde passa a ser exatamente 1. Confira a regra \(\varepsilon\,b^2/|ac|\): para \(b=10^6\) ela prevê \(10^{-4}\); para \(b=10^8\), cerca de \(1\). Depois tente \(a=1\), \(b=-10^8\), \(c=1\): qual raiz é a ruim agora, e por quê?

Tome \(x^2+2x+5\), cujo discriminante é \(-16\). O laboratório mostra \(x=-1\pm2i\); confira que \((x+1)^2=-4\) na forma do quinto quadro, e que \((-1+2i)+(-1-2i)=-2=-b/a\) e \((-1+2i)(-1-2i)=5=c/a\). Qual é o vértice da parábola, e como a altura dele diz que as raízes não são reais?