Matemática Computacional

Resumo da matéria [Versão 3.1]

[Capítulo 1] [Capítulo 2] [Capítulo 3] [Capítulo 4] [Capítulo 5]

Capítulo IV

Aproximação de Funções

Este capítulo será dividido em três secções principais:

(i) Interpolação
(ii) Mínimos Quadrados
(iii) Integração Numérica

A aproximação de funções pode ser considerada para substituir o cálculo moroso de uma função por funções mais simples. Num outro caso, quando não dispomos da função original, estes métodos permitem considerar uma aproximação da função baseada apenas em alguns pontos.

Interpolação

A interpolação polinomial é uma noção simples que apareceu naturalmente com o desenvolvimento do cálculo. Até há 50 anos atrás, os processos de cálculo eram limitados ao engenho manual, e a interpolação polinomial surgia como uma via simples na aproximação de funções complexas que eram normalmente listadas por tabelas, em alguns valores. Associada a Newton está uma fórmula recursiva para polinómios interpoladores, mas foi com Lagrange e Hermite que estas noções levaram a alguma sistematização. Nomeadamente, no final do séc. XIX, Runge e Chebyshev consideraram questões de instabilidade da aproximação polinomial.
Convém ainda referir que a interpolação para funções com mais que uma variável é ainda um problema actual, que surge em várias aplicações. Num exemplo simples de design automóvel foram desenvolvidas nos anos 60, na Renault, as curvas de Bézier. Grande parte da modelação actual assenta em elementos finitos enquanto base para a interpolação 2D ou 3D.

Iremos concentrar-nos apenas na interpolação dos valores de uma função (interpolação de Lagrange), sem ter em conta eventuais valores das suas derivadas (interpolação de Hermite).
Começamos por considerar um quadro geral, ainda assim apenas linear e para funções de uma variável. Concentramo-nos depois nos casos particulares de interpolação polinomial e trigonométrica.

Interpolação geral

O caso geral aqui considerado assume que dispomos informação num número finito de pontos,

(x₀, y₀), ..., (x_n, y_n)∈ R².

pretendendo encontrar uma função g num certo subespaço, tal que o gráfico de g intersecte estes pontos.

Os valores das abcissas x₀, ..., x_n são designados nós de interpolação e são todos distintos.
Os valores das ordenadas são os valores dados e que podem corresponder aos valores que uma função f, ou seja
y₀ = f(x₀), ..., y_n = f(x_n).

O objectivo será:

Dado um subespaço de funções
S = < φ₀, ⋯ ,φ_m >
(as funções φ₀, ⋯ ,φ_m devem constituir uma base);
Encontrar uma função g∈ S, ou seja,
g(x) = a₀ φ₀(x) + ⋯ +a_m φ_m(x)
tal que g verifique
g(x₀) = y₀, ..., g(x_n) = y_n,
ou seja, g interpola os pontos dados − ou ainda, g interpola f nos nós dados.

O problema é linear e representado pelo sistema

╭
│	a₀ φ₀(x₀) + ⋯ +a_m φ_m(x₀)	= y₀
┤	⋮	⋮
│	a₀ φ₀(x_n) + ⋯ +a_m φ_m(x_n)	= y_n
╰

onde há n+1 equações e os m+1 coeficientes a₀, ⋯, a_m são as incógnitas.
Para que o sistema seja bem determinado consideramos m=n (nº funções base = nº pontos) e podemos representar matricialmente

┌				┐
│	φ₀(x₀)	⋯	φ_n(x₀)	│
│	⋮	⋱	⋮	│
│	φ₀(x_n)	⋯	φ_n(x_n)	│
└				┘

┌		┐
│	a₀	│
│	⋮	│
│	a_n	│
└		┘

┌		┐
│	y₀	│
│	⋮	│
│	y_n	│
└		┘

Se as funções base forem linearmente independentes no conjunto dos nós dados, então a matriz é invertível, e haverá solução única.

Relativamente ao caso geral o problema resume-se à resolução deste sistema (usualmente designado sistema de Vandermonde). No entanto, podemos ver que a solução deste sistema pode ser dada de forma explícita no caso da interpolação polinomial (as funções base são monómios).

Interpolação Polinomial

Considerando como funções base monómios,

φ₀(x) = 1, φ₁(x) = x, φ₂(x) = x², ⋯, φ_m(x) = x^m;

o subespaço S será um espaço dos polinómios de grau ≤ m,

g(x) = a₀ + a₁x + a₂x² + ⋯ + a_mx^m .

Assim, o sistema é dado por

┌					┐
│	1	x₀	⋯	(x₀)ⁿ	│
│	1	x₁	⋯	(x₁)ⁿ	│
│	⋮	⋮	⋱	⋮	│
│	1	x_n	⋯	(x_n)ⁿ	│
└					┘

┌		┐
│	a₀	│
│	a₁	│
│	⋮	│
│	a_n	│
└		┘

┌		┐
│	y₀	│
│	y₁	│
│	⋮	│
│	y_n	│
└		┘

e podemos ver que no caso polinomial a matriz de Vandermonde é sempre invertível, desde que os nós de interpolação sejam distintos.

Teorema:
Dados n+1 nós de interpolação distintos, existe um e um só polinómio interpolador de grau ≤ n.

dem: Como o número de equações é igual ao número de incógnitas, resta ver que há unicidade para concluir que a matriz de Vandermonde é invertível.
Supondo que existiam dois polinómios interpoladores p_n e q_n com grau ≤ n. Então

p_n(x₀) = y₀ = q_n(x₀)	⇒ x₀ é raiz de p_n − q_n
⋮	⋮
p_n(x_n) = y_n = q_n(x_n)	⇒ x_n é raiz de p_n − q_n

Isto significa que o polinómio p_n − q_n (com grau ≤ n) tem n+1 raízes, logo só pode ser nulo e assim p_n ≡ q_n.

Fórmula de Lagrange

Podemos obter uma expressão explícita para o polinómio interpolador sem ser necessário resolver o sistema de Vandermonde, poupando ainda no número de operações.

Dado o conjunto dos nós de interpolação x₀, ..., x_n definimos os polinómios base de Lagrange L_i: polinómios de grau n que verificam

L_i(x_j) = δ_ij =

╭	1	se i=j
┤
╰	0	se i≠ j

Pela definição, o polinómio L_i terá x₀, ..., x_n como raízes − à excepção de x_i. Isto significa que

L_i(x) = C

n
∏
j=0, j≠ i

(x − x_j)

A condição no índice i permite retirar o valor de C, pois

1 = L_i(x_i) = C

n
∏
j=0, j≠ i

(x_i − x_j) ⇒ C = 1 /

n
∏
j=0, j≠ i

(x_i − x_j)

Concluímos assim que

L_i(x) =

n
∏
j=0, j≠ i

x − x_j

x_i − x_j

Proposição: O polinómio interpolador é dado pela Fórmula de Lagrange:

p_n(x) =

n
∑
i=0

y_i L_i(x)

ou seja,

p_n(x) =

n
∑
i=0

y_i

n
∏
j=0, j≠ i

x − x_j

x_i − x_j

dem: Resulta da expressão dos L_i, uma vez que

p_n(x_j) =

n
∑
i=0

y_i L_i(x_j) =

n
∑
i=0

y_i δ_ij = y_j,

e portanto p_n interpola os valores dados.

Vimos dois processos de obter o polinómio interpolador, e iremos ver ainda um terceiro, através de diferenças dividas − a designada Fórmula de Newton.

Diferenças divididas

Uma diferença dividida de uma função f num conjunto de nós {x₀, ⋯, x_n} é denotada

f[x₀,⋯,x_n]

e dizemos que tem ordem n.

Podemos encarar como ordem 0, o valor da função no ponto f[x₀] = f(x₀).
Para ordem 1, corresponde à noção de razão incremental, e para ordens superiores generalizam esta noção.

Dados dois nós x₀, x₁ a diferença dividida de ordem 1 de f é
f[x₀,x₁] = f[x₀] − f[x₁]
x₀ − x₁
dados três nós x₀, x₁, x₂ a diferença dividida de ordem 2 de f é
f[x₀,x₁,x₂] = f[x₀,x₁] − f[x₁,x₂]
x₀ − x₂
... e assim recursivamente:
Dados n+1 nós x₀, ⋯, x_n a diferença dividida de ordem n de f é

f[x₀,⋯,x_n] = f[x₀,⋯,x_n-1] − f[x₁,⋯,x_n]
x₀ − x_n

A diferença divida de ordem n pode ainda ser caracterizada por ser o coeficiente de maior grau num polinómio de grau n.

Proposição: Para um qualquer polinómio p_n=c₀+⋯+c_nxⁿ temos

p_n[x₀,⋯, x_n] = c_n

(quaisquer que sejam os nós distintos x₀,⋯, x_n).

dem: Quando n=0 é trivial. Argumentamos por indução para n ≥ 1.
Consideramos dois polinómios interpoladores de p_n:

q_n-1^(a)(x) = c_n-1^(a) x^n-1+⋯ polinómio de grau ≤ n-1 interpolador de p_n nos nós x₀,⋯, x_n-1;
q_n-1^(b)(x) = c_n-1^(b) x^n-1+⋯ polinómio de grau ≤ n-1 interpolador de p_n nos nós x₁,⋯, x_n;

Primeiro vemos que

p_n(x) =

(x − x_n)q_n-1^(a)(x) − (x − x₀)q_n-1^(b)(x)

x₀ − x_n

pois a igualdade verifica-se para x = x₀, ⋯, x_n e o polinómio interpolador é único. Por hipótese de indução,

c_n-1^(a) = q_n-1^(a)[x₀,⋯, x_n-1] = p_n[x₀,⋯, x_n-1]
c_n-1^(b) = q_n-1^(b)[x₁,⋯, x_n] = p_n[x₁,⋯, x_n]

porque q_n-1^(a,b) interpolam p_n. É assim claro que

c_n =

c_n-1^(a)-c_n-1^(b)

x₀ − x_n

p_n[x₁,⋯, x_n] − p_n[x₀,⋯, x_n-1]

x₀ − x_n

= p_n[x₀,⋯, x_n]

Observações:
(1) Quando p_n é o polinómio interpolador de f nos nós x₀,⋯, x_n temos

c_n = p_n[x₀,⋯, x_n] = f[x₀,⋯, x_n]

pois os valores de f coincidem com os valores p_n nesses nós.
(2) Consequentemente, e como o polinómio interpolador é único, é indiferente a ordem pela qual se consideram os nós nas diferenças divididas, ou seja

f[x₀,⋯, x_n] = f[x_σ(0),⋯, x_σ(n)]

onde σ é uma permutação em { 0, ⋯ , n } (por exemplo, f[x₀,x₁, x₂] = f[x₂,x₀, x₁]).

Para simplificar, assumimos uma ordenação dos nós x₀ < ⋯ < x_n.

Teorema: Seja f ∈ Cⁿ[x₀, x_n], então existe ξ ∈ [x₀, x_n] :

f[x₀, ⋯ , x_n] =

f⁽ⁿ⁾(ξ)

dem: Consideramos p_n o polinómio interpolador de f nos nós x₀,⋯, x_n.
Então E_n = f − p_n tem n+1 zeros: x₀ < ⋯ < x_n.
Pelo T. Rolle, E_n' tem n zeros: ξ₀ < ⋯ < ξ_n em [x₀, x_n]
... e assim sucessivamente, até que (E_n)⁽ⁿ⁾ tem um zero: ξ∈ [x₀, x_n].
Ou seja,

∃ ξ∈ [x₀, x_n]: 0 = (E_n)⁽ⁿ⁾(ξ) = f⁽ⁿ⁾(ξ) − c_n n!

onde c_n é coeficiente de maior grau de p_n, e como p_n é interpolador temos c_n = f[x₀,⋯, x_n].

Observação:
Este resultado torna evidente a relação entre as diferenças divididas de ordem n e as derivadas de ordem n. Para além disso, é conveniente notar que no limite faz sentido considerar o mesmo ponto, pois

f[z,z] = lim_{y→ z} f[y,z] = f'(z)

e assim é possível utilizar diferenças divididas mesmo com nós repetidos.
Para n+1 repetições de x, obtemos

f[z,⋯, z] = f⁽ⁿ⁾(z)/n!.

Fórmula de Newton

Como consequência das noções anteriores, é fácil obter um outro processo para escrever o polinómio interpolador − a denominada Fórmula de Newton.

De facto, se considerarmos p_n-1 o polinómio interpolador para os nós x₀, ⋯, x_n-1 podemos obter:

p_n(x) = p_n-1(x) + f[x₀, ⋯, x_n-1](x − x₀)⋯(x − x_n-1)

porque

p_n(x_i) = p_n-1(x_i) = y_i para i = 0, ⋯, n-1.
o coeficiente de maior grau é f[x₀, ⋯, x_n-1];

Assim, este processo pode ser aplicado recursivamente até p₀(x)=f[x₀].
Abreviadamente, temos assim a Fórmula de Newton:

p_n(x) =

n
∑
k=0

f[x₀,⋯,x_k] (x − x₀)⋯ (x − x_k-1)

Observação:
Reparamos que se trata de uma expansão − associada à expansão de Taylor. De facto, se considerarmos que x_k→ x₀ então

p_n(x) =

n
∑
k=0

f[x₀,⋯,x₀] (x − x₀)^k =

n
∑
k=0

f⁽ⁿ⁾(x₀)/n! (x − x₀)^k

usando a relação entre derivadas e diferenças divididas. Vemos que p_n é o polinómio de Taylor correspondente à expansão de f em x₀.

Erro de interpolação polinomial

Conhecendo f, o erro de interpolação é definido como uma função E=f − p_n , definida num qualquer ponto de avaliação z∈ [x₀, x_n]

E(z) = f(z) − p_n(z).

Convém referir que podemos ainda considerar z∉ [x₀, x_n], ou seja z está fora do intervalo onde estão os nós − este processo é normalmente designado extrapolação.
Escrevemos [x₀, x_n; z] para designar o intervalo que contém os nós de interpolação e o ponto de avaliação.

Teorema: Temos as seguintes expressões para o erro de interpolação:

E(z) = f[x₀,⋯, x_n,z] (x − x₀)⋯(x − x_n)
Se f ∈ Cⁿ⁺¹[x₀, x_n; z], então existe ξ ∈ [x₀, x_n; z] :

E(z) = f⁽ⁿ⁺¹⁾(ξ)
(n+1)! (x − x₀)⋯(x − x_n)
Conclui-se a estimativa de erro:
|E(z)| ≤ max_{x∈[x₀, x_n; z]} |f⁽ⁿ⁺¹⁾(x)|
(n+1)! |z − x₀|⋯|z − x_n|

dem:
Sendo p_n o polinómio interpolador nos nós x₀,⋯, x_n, acrescentado z como novo nó, obtemos pela Fórmula de Newton:

p_n+1(x) = p_n(x) + f[x₀,⋯, x_n,z](x − x₀)⋯(x − x_n)

e como p_n+1(z) = f(z), temos

f(z) = p_n(z) + f[x₀,⋯, x_n,z](z − x₀)⋯(z − x_n)

concluindo-se a primeira igualdade.
Quando f ∈ Cⁿ⁺¹[x₀, x_n; z], então pela relação entre as diferenças divididas e as derivadas

∃ ξ ∈ [x₀, x_n; z] : f[x₀,⋯, x_n, z] = f⁽ⁿ⁺¹⁾(ξ)/(n+1)!

ficando demonstrada a segunda igualdade.

Exemplo:
Considere a função dada por uma série

f(x) =

∞
∑
k=1

1/(k³+x)

que assume os seguintes valores

x	1	2	3	4	5
f(x)	0.686503	0.507161	0.413262	0.354294	0.313232

Esta função não é fácil de calcular, pois envolve uma soma infinita em que a rapidez de convergência não é muito significativa.
Esta é uma situação em que pode ser aconselhável substituir o cálculo de f por uma aproximação dada por um polinómio interpolador.
Neste caso, o polinómio interpolador p₄ é dado por

p₄(x) = 1.03529 − 0.469887 x + 0.142074 x² − 0.0223729 x³ + 0.0013954 x⁴

podendo ser obtido pela fórmula de Lagrange ou pela fórmula de Newton.
A derivação termo a termo, permite obter

f⁽⁵⁾(x) =

∞
∑
k=1

-120(k³+x)^-6 (a série é absolutamente convergente)

e temos

max_{x∈ [1,5]} |f⁽⁵⁾(x)| = |f⁽⁵⁾(1)| ≈ |-1.88| ≤ 2.

Conclui-se assim que

|E(z)| ≤

|z − x₀|⋯|z − x_n|

Por exemplo, para z=3.5 obtemos |E(z)| ≤ 1.40625/60 = 0.02344, mas podemos mesmo constatar que a aproximação é melhor (tem sempre erros relativos inferiores a 0.5% nesse intervalo). Ver Figura 1.

Figura 1: Erro relativo entre a função f e a aproximação p₄.
Notamos ainda que, conforme esperado, o erro é nulo nos nós de interpolação, e por outro lado tem maior amplitude próximo das extremidades do intervalo [1,5] do que no centro.

Observação complementar − Nós de Chebyshev

Quando procuramos aproximar uma função num intervalo [a,b], podemos considerar a interpolação usando n+1 nós igualmente espaçados,

x_k = a + kh, com h=(b-a)/n,

mas esta escolha não se revela a mais eficaz.
Isto é ilustrado pelo exemplo de Runge, para a função

f(x) =

1+(5x)²

no intervalo [-1,1]

Como podemos ver na Figura 2A, a aproximação com n=16 é boa no centro mas detriora-se rapidamente nas extremidades do intervalo [-1,1], pois o polinómio de grau 16 apresenta grandes oscilações.

Figura 2A: Exemplo de Runge. Utilizando nós igualmente espaçados, o polinómio p₁₆ (a vermelho) apresenta uma má aproximação da função f (a preto) nas extremidades do intervalo [-1,1].

O problema da má aproximação nas extremidades tem a ver com o comportamento da expressão que intervem no erro de interpolação

|W(x)| = |x − x₀| ⋯ |x − x_n|

Esta expressão não depende da função f, depende apenas da escolha dos nós, pelo que (conforme Chebyshev) constata-se que o valor |W| pode ser minimizado com uma escolha apropriada de nós no intervalo [-1,1]:

c_k = cos (

(2k+1)π

2n+2

) (k=0,⋯,n)

são os denominados nós de Chebyshev, que são as n+1 raízes de T_n+1, polinómios de Chebyshev, cuja expressão pode ser dada trigonometricamente

T_n(x) = cos( n arccos(x)).

É possível mostrar que

T_n+1(x) = 2ⁿ (x − c₀) ⋯ (x − c_n)

e assim concluir que para os nós de Chebyshev,

|W(x)| = |x − c₀| ⋯ |x − c_n| ≤ 2^-n|T_n+1(x)| ≤ 2^-n

sendo esta escolha de nós a que permite minimizar |W| na expressão do erro (no intervalo [-1,1] − a passagem para outro intervalo pode ser feita trivialmente por mudança de variável).
Na Figura 2B vemos a aproximação obtida considerando os nós de Chebyshev, que suprime o problema de instabilidade nas extremidades (por outro lado, piorando ligeiramente a aproximação central).

Figura 2B: Exemplo de Runge. Utilizando nós de Chebyshev, o polinómio p₁₆ passando pelos nós de Chebyshev (verde) corrige a aproximação da função f (a preto) nas extremidades do intervalo [-1,1].

Observação complementar − Interpolação trigonométrica

A interpolação trigonométrica consiste numa aproximação em que as funções base são senos e co-senos, ao invés de monómios. No entanto, dada a relação trigonométrica nos complexos,

x = exp(It) = cos(t) + I sin(t) ⇒ x^m = exp(Imt) = cos(mt) + I sin(mt)

a sua dedução pode ser facilitada considerando o contexto complexo.
(Para evitar confusões com índices i, escrevemos aqui o número imaginário como I)

O intervalo de referência será agora [0, 2π[ e as funções a aproximar são entendidas como periódicas (ie. f(0)=f(2π)), ainda que isso não seja necessário (trabalhamos apenas com valores discretos).
Consideramos

Os nós t₀< ⋯ < t_n-1, em [0, 2π[
As funções base monomiais complexas
φ_k = x^k = exp(I kt)
procurando uma aproximação da forma
q(t) = a₀φ₀(t) + ⋯ + a_n-1φ_n-1(t) = a₀ + a₁ e^It ⋯ + a_n-1 e^I(n-1)t.

A matriz de Vandermonde fica assim

┌					┐
│	1	x₀	⋯	(x₀)^n-1	│
│	1	x₁	⋯	(x₁)^n-1	│
│	⋮	⋮	⋱	⋮	│
│	1	x_n-1	⋯	(x_n-1)^n-1	│
└					┘

┌					┐
│	1	exp(I t₀)	⋯	exp(I (n-1)t₀)	│
│	1	exp(I t₁)	⋯	exp(I (n-1)t₀)	│
│	⋮	⋮	⋱	⋮	│
│	1	exp(I t_n-1)	⋯	exp(I (n-1)t_n-1)	│
└					┘

Os nós igualmente espaçados são aqui apropriados (... as funções base são outras),

t_k=2kπ/n

pois a matriz de Vandermonde W fica simétrica

W_ij = (x_i)^j = exp(I t_ij) = exp(2π I ij/n)

ou seja,

┌					┐
│	1	1	⋯	1	│
│	1	exp(2π I /n)	⋯	exp(2π I (n-1)/n)	│
│	⋮	⋮	⋱	⋮	│
│	1	exp(2π I (n-1)/n)	⋯	exp(2π I (n-1)²/n)	│
└					┘

e sendo a=(a₀,⋯,a_n-1), y=(y₀,⋯,y_n-1) obtemos o sistema linear

W a = y

A matriz W tem uma propriedade muito útil que relaciona a sua inversa com a adjunta W^∗ = W^T.

Teorema: W W^∗ = n I (em que I é a matriz identidade).

dem:
Basta reparar que

[W W^∗]_ij =

n-1
∑
k=0

W_ik W_kj =

n-1
∑
k=0

exp(2π I ik/n) exp(-2π I kj/n) =

n-1
∑
k=0

exp(2π I k(i-j)/n)

definindo z_ij = exp(2π I (i-j)/n) ≠ 1, se i ≠ j, temos a soma da progressão geométrica

[W W^∗]_ij =

n-1
∑
k=0

(z_ij)^k =

1-z_ijⁿ

1-z_ij

= 0

reparando que z_ijⁿ=exp(2π I (i-j))=1.
Quando i=j temos z_ij = 1 e os termos da diagonal [W W^∗]_ii = n× 1 = n.

Concluímos assim que a matriz W é invertível com

W^-1 = n^-1 W^∗.

Associada a esta noção temos uma outra:

Transformada de Fourier Discreta (TFD ou DFT): A transformada de Fourier discreta F transforma um vector
v=(v₀, ⋯, v_n-1) num vector com as mesmas dimensões dado pelas componentes

[ Fv ]_k =

n-1
∑
j=0

v_j e^{-2π I jk/n}, (k=0,⋯,n-1).

Ou seja,

Fv = W^∗ v.

... e por outro lado concluímos que a inversa da TFD é simplesmente

F^-1v = n^-1 W v;

(alguns autores partilham o termo n^-1, considerando n^-1/2 na definição de directa e inversa).

Conclusão: A solução do problema de interpolação trigonométrica pode ser dada aplicando directamente a transformação de Fourier discreta sobre o vector dos dados y:

a = W^-1y = n^-1 W^∗y = n^-1 Fy.

Observação:
1) Ao considerar o contexto complexo, os coeficientes são complexos. Há que ter algum cuidado extra na passagem para uma aproximação que envolva valores só reais, apenas com senos e co-senos (nomeadamente considerar um número ímpar de pontos...).
2) O número de operações na DFT é O(n²) mas pode ser reduzido para O(n log(n)) através do algoritmo FFT (Fast Fourier Transform) de Cooley-Tukey (apresentado nos anos 60, mas já era do conhecimento de Gauss...).

Método dos Mínimos Quadrados

Regressamos ao contexto da interpolação geral em que dispomos informação num número finito de pontos,

(x₀, y₀), ..., (x_n, y_n)∈ R².

pretendendo encontrar uma função g∈ S,

S = < φ₀, ⋯ ,φ_m >

ou seja,

g(x) = a₀ φ₀(x) + ⋯ + a_m φ_m(x),

mas em que m < n.

A grande diferença é que como m < n agora não podemos exigir que a função g passe pelos pontos dados. Apenas pretendemos minimizar (tanto quanto possível) a distância a esses pontos.

Para esse efeito temos que definir uma distância, e isso será feito através de um produto interno, relembrando que

dist(f,g) = ∥f − g∥ = √ < f − g, f − g >,

onde <⋅,⋅> representa o produto interno.

No caso dos mínimos quadrados discretos, usamos o produto interno habitual (para valores reais)

< u, v > =

n
∑
k=0

u(x_k) v(x_k)

em que x₀, ⋯, x_n são os nós de colocação (distinguimos colocação de interpolação, pois não exigimos que haja interpolação).

Devemos ainda considerar que os valores y₀, ⋯, y_n podem resultar de uma função f, ou seja

y_k = f(x_k) (para k=0,⋯,n)

e além disso, considerar mesmo a possibilidade de minimizar não apenas num conjunto finito de pontos, mas uma função f em todo um intervalo [a,b].
Nesse sentido falaremos do caso contínuo, e o produto interno será em L²[a,b], dado por

< u, v > = ∫_[a,b] u(x) v(x) dx

Observação: No caso complexo, o produto interno deve incluir o conjugado complexo, ou seja:

( i) caso discreto:

< u, v > =

n
∑
k=0

u(x_k) v(x_k),

(ii) caso contínuo:
< u, v > = ∫_[a,b] u(x) v(x) dx.

Sistema Normal

Em qualquer das situações (caso discreto ou contínuo) podemos desenvolver o processo de minimização de forma abstracta, sem ter em conta o produto interno específico, que está em causa.

Sendo

g = a₀ φ₀ + ⋯ + a_m φ_m,

o objectivo será minimizar a distância de f a g, ou seja, minimizar o resíduo

Q = dist(f,g) = ∥f − g∥ = √ < f − g, f − g >,

em termos dos coeficientes a₀, ⋯, a_m.

Para esse efeito, pretendemos encontrar o ponto de mínimo de Q² (será também de Q) que é ponto crítico, verificando
∇ Q² = 0, ou seja,

∂ Q²

∂ a₀

= 0 , ... ,

∂ Q²

∂ a_m

= 0

Isto corresponde a um sistema de equações

∂

∂ a_i

< f − g, f − g > = 0 , (i = 0, ... m) ,

(onde a dependência nos a_i está exclusivamente em g.
Pela regra da derivação do produto temos

< 0 − ∂ g/∂ a_i, f − g > + < f − g, 0 − ∂ g/∂ a_i > = 0 , (i = 0, ... m) ,

e como ∂ g/∂ a_i = φ_i e o produto comuta (nos reais), temos

-2 < f − g , φ_i > = 0 , (i = 0, ... m) ,

ou ainda

< g , φ_i > = < f , φ_i > , (i = 0, ... m) .

Assim, explicitando g obtemos o sistema de equações (denominado sistema normal)

a₀< φ₀, φ_i > + ⋯ + a_m< φ_m, φ_i > = < f , φ_i > , (i = 0, ... m) ,

que pode ser representado na forma matricial

┌				┐
│	<φ₀,φ₀>	⋯	<φ_m,φ₀>	│
│	⋮	⋱	⋮	│
│	<φ₀,φ_m>	⋯	<φ_m,φ_m>	│
└				┘

┌		┐
│	a₀	│
│	⋮	│
│	a_m	│
└		┘

┌		┐
│	< f, φ₀>	│
│	⋮	│
│	< f, φ_m>	│
└		┘

em a matriz S do sistema verifica:

Proposição: Se as funções base forem linearmente independentes, a matriz do sistema normal é simétrica e definida positiva.

dem:
A comutatividade do produto interno real implica a simetria.
A matriz S verifica (para qualquer v≠ 0)

v^T S v =

n
∑
i,j=0

v_iS_ijv_j =

n
∑
i,j=0

v_i <φ_i,φ_j> v_j = < g , g > , com g=

n
∑
k=0

v_kφ_k.

Portanto a matriz é definida positiva:

v^T S v = ∥g∥² > 0,

notando ainda que as funções base são linearmente independentes, e g=0 implicaria v=0.

Caso Discreto

No caso discreto, o sistema normal pode ainda ser obtido através do sistema sobredeterminado associado à interpolação.
Com efeito, o sistema sobredeterminado será

M a = y ⇔

┌				┐
│	φ₀(x₀)	⋯	φ_m(x₀)	│
│	⋮	⋱	⋮	│
│	⋮	⋱	⋮	│
│	φ₀(x_n)	⋯	φ_m(x_n)	│
└				┘

┌		┐
│	a₀	│
│	⋮	│
│	a_m	│
└		┘

┌		┐
│	y₀	│
│	⋮	│
│	⋮	│
│	y_n	│
└		┘

e multiplicando em ambos os lados por M^∗ obtemos um sistema m × m:

M^∗M a = M^∗y

Com efeito a matriz S é M^∗M, pois

[M^∗M]_ij =

n
∑
k=0

φ_i(x_k)φ_j(x_k) = < φ_i, φ_j > = S_ij

e da mesma forma

[M^∗y]_i =

n
∑
k=0

φ_i(x_k) f(x_k) = < f, φ_i > .

Exemplo (Regressão Linear):
Consideramos os pontos dados por

(x_k, f_k) (com k = 0,⋯,N)

e procuramos aproximar estes valores por um polinómio de 1º grau

g(x) = a₀ + a₁ x ⇒ φ₀(x)=1, φ₁(x)=x

Obtemos o sistema

┌			┐
│	< 1,1>	< x,1>	│
│	< 1,x>	< x,x>	│
└			┘

┌		┐
│	a₀	│
│	a₁	│
└		┘

┌		┐
│	< f, 1>	│
│	< f, x>	│
└		┘

onde

< 1, 1 > = N+1, < 1, x > = < x, 1 > =

N
∑
k=0

x_k = (N+1)E(x)

(a notação E(x) designa a média de x),

< x, x > =

N
∑
k=0

x_k² = (N+1)E(x²)

e ainda

< f, 1 > =

N
∑
k=0

f_k = (N+1)E(f), < f, x > =

N
∑
k=0

x_k f_k = (N+1)E(xf) .

Com as novas notações, o sistema fica

┌			┐
│	1	E(x)	│
│	E(x)	E(x²)	│
└			┘

┌		┐
│	a₀	│
│	a₁	│
└		┘

┌		┐
│	E(f)	│
│	E(xf)	│
└		┘

e a solução do sistema é simples, tendo-se

(a₀, a₁) =

( E(x²)E(f) − E(x)E(xf) , E(xf) − E(x)E(f) )

E(x²) − E(x)²

em que o denominador é a variância σ² = Var(x) = E(x²) − E(x)².

Caso Contínuo

Este caso é em tudo semelhante ao anterior, variando apenas o produto interno.

Exemplo 1: Um caso simples é aquele em pretendemos aproximar uma função por um polinómio, em que as funções base são definidas pelos monómios no intervalo [0,1]:

φ_k(x) = x^k, com k=0,⋯,m

Então obtemos

S_ij = < φ_i , φ_j > = ∫_[0,1] xⁱ x^j dx =

i+j+1

que é a denominada matriz de Hilbert:

┌				┐
│	1	⋯	1/m	│
│	⋮	⋱	⋮	│
│	1/m	⋯	1/(2m+1)	│
└				┘

As matrizes de Hilbert são mal condicionadas, tendo-se já para m=3 um valor elevado para o número de condição, Cond₁(S) = 28375, (e para m=4 já atinge 943656).
Os termos do segundo membro são dados pelos integrais

< f, xⁱ> =∫_[0,1] xⁱ f(x) dx.

Dependendo da complexidade de f poderá ser necessário calcular estas componentes através de integração numérica, se a primitivação não for exequível.

Exemplo 2: Suponhamos que pretendemos aproximar a função

f(x) = exp(x), no intervalo [0,1]

usando um polinómio de grau 1, g(x)=a+bx. Então basta resolver um sistema 2 × 2

┌				┐
│	1	⋯	1/2	│
│	1/3	⋯	1/3	│
└				┘

┌		┐
│	a	│
│	b	│
└		┘

┌		┐
│	∫_[0,1] exp(x) dx = e − 1	│
│	∫_[0,1] x exp(x) dx = 1	│
└		┘

A solução deste sistema dá a melhor aproximação

g(x) = 4e − 10 + (18 − 6e)x

que minimiza o resíduo Q, ou seja:

Q ≥ ( ∫_[0,1] (f(x) − g(x))² dx ) ^1/2 = 0.06277…

Observação: (Polinómios Ortogonais)
Como vimos há um problema de condicionamento nas matrizes do sistema normal, pelo que é conveniente poder evitar a resolução desse sistema. Um processo consiste em definir a matriz já na forma diagonal, bastando para isso que as funções base sejam ortogonais.
No entanto, se encontrar uma base de funções ortogonal é sempre possível pelo processo de ortogonalização de Gram-Schmidt, isso não é sempre praticável, sendo preferível ter um cálculo prévio.
Isso é exequível em certas situações particulares, como é o caso dos Polinómios de Legendre que são ortogonais para o produto interno contínuo no intervalo [-1,1}.
Os Polinómios de Legendre P_k podem ser dados de forma recursiva pela expressão

P₀(x) = 1 , P₁(x) = x,
P_n+1(x) = (2n+1) x P_n(x) − n P_n-1(x)
n+1

cujos termos seguintes são:

P₂(x) = (3x² − 1)/2, P₃(x) = (5x³ − 3x)/2, ⋯

É ainda possível definir outros polinómios ortogonais, adaptados a outros produtos internos definidos com pesos de ponderação W.
Por exemplo, considerando uma função peso W positiva, definimos um produto interno

< u , v >_W = ∫_[-1,1] W(x) u(x) v(x) dx.

O caso mais interessante é considerar como peso

W(x) = 1 / √ 1 − x²

pois isso leva aos polinómios ortogonais de Chebyshev, que também podem ser descritos recursivamente

T₀(x) = 1, T₁(x) = x,

T_n+1(x) = 2 x T_n(x) − T_n-1(x)

tendo-se sucessivamente T₂(x) = 2x² − 1, T₃(x) = 4x³ − 3x, ...
As raízes destes polinómios vão dar os nós de Chebyshev...
Referimos ainda que para outras funções de peso W é possível encontrar uma fórmula recursiva geral que gera uma base de polinómios ortogonais.
Sinteticamente, essa formula que gera polinómios ortogonais G_k é dada por

G_k+1(x) = (x − < x G_k, G_k >_W ∥G_k∥_W^-2 ) G_k(x) − ∥G_k∥_W² ∥G_k-1∥_W^-2 G_k-1(x).

Integração Numérica

Um outro problema na aproximação de funções, diz respeito à aproximação de funções definidas por integrais. Primeiro, porque poderá acontecer que o integral não tenha uma primitiva explícita ou conhecida, depois porque esse cálculo pode ser computacionalmente moroso e menos eficaz do que a aplicação de regras numéricas.

Seja f∈ C[a,b]. Para a aproximação do integral

I(f) = ∫_[a,b] f(x) dx,

definimos uma fórmula de integração (ou regra de quadratura)

Q(f) = w₀ f(z₀) + ⋯ + w_m f(z_m)

em que designamos w_k como pesos de integração e z_k como nós de integração.
Estes valores − pesos e nós − podem ser calculados de forma a que o resultado Q(f) seja exacto, igual a I(f), para certas funções base f.
Normalmente usam-se como funções base os monómios, pretendendo que a fórmula seja exacta para polinómios de um determinado grau, levando à seguinte definição.

Definição: Dizemos que uma regra de quadratura Q tem pelo menos grau p se for exacta para polinómios de grau ≤ p.

Ou seja, a regra tem pelo menos grau p se for exacta para monómios x^k,

I(x^k) = Q(x^k) (para k=0,⋯,p).

Terá exactamente grau p se tiver pelo menos grau p e não tiver grau p+1, ou seja, se verificar I(x^p+1)≠ Q(x^p+1).

Esta definição leva ainda a um método para definir as regras de quadratura − o denominado Método dos Coeficientes Indeterminados.

Método dos Coeficientes Indeterminados

A verificação de que uma regra tem grau p leva a p+1 equações:

(P){

I(1) = Q(1) ⇔	b − a	= ∫_[a,b] 1 dx = w₀ + ⋯ + w_m
I(x) = Q(x) ⇔	(b² − a²)/2	= ∫_[a,b] x dx = w₀z₀ + ⋯ + w_mz_m
⋮	⋮
I(x^p) = Q(x^p) ⇔	(b^p+1 − a^p+1)/(p+1)	= ∫_[a,b] x^p dx = w₀z₀^p + ⋯ + w_pz_m^p

onde podemos distinguir duas situações:

os nós de integração são conhecidos − sistema linear nos pesos;
os nós e os pesos de integração são desconhecidos − sistema não linear em ambos;

Caso Linear.
Começamos por analisar o caso do sistema linear, supondo que os nós

z₀, ⋯, z_m são conhecidos/dados.

Então, para que a regra de quadratura tenha pelo menos grau m, obtemos o sistema linear

┌					┐
│	1	1	⋯	1	│
│	z₀	z₁	⋯	z_m	│
│	⋮		⋱	⋮	│
│	(z₀)^m	(z₁)^m	⋯	(z_m)^m	│
└					┘

┌		┐
│	w₀	│
│	w₁	│
│	⋮	│
│	w_m	│
└		┘

┌		┐
│	b − a	│
│	(b² − a²)/2	│
│	⋮	│
│	(b^m+1 − a^m+1)/(m+1)	│
└		┘

em que a matriz (m+1) × (m+1) é a transposta da matriz de Vandermonde (e portanto é invertível).

Exemplo − Regra dos Trapézios
Considerando z₀=a, z₁=b obtemos um sistema linear 2 × 2:

┌ ┐
│ 1 1 │
│ a b │
└ ┘

┌ ┐
│ w₀ │
│ w₁ │
└ ┘
=
┌ ┐
│ b − a │
│ (b² − a²)/2 │
└ ┘
com solução w₀=w₁=(b − a)/2 o que dá

T(f) = b − a
2 ( f(a) + f(b) ) [Regra dos Trapézios].

Exemplo − Regra de Simpson
Considerando z₀=a, z₁=(a+b)/2, z₂=b obtemos um sistema linear 3 × 3:

┌				┐
│	1	1	1	│
│	a	(a+b)/2	b	│
│	a²	(a+b)²/4	b²	│
└				┘

┌		┐
│	w₀	│
│	w₁	│
│	w₂	│
└		┘

┌		┐
│	b − a	│
│	(b² − a²)/2	│
│	(b³ − a³)/3	│
└		┘

com solução w₀=(b − a)/6, w₁=4w₀, w₂=w₀ o que dá

S(f) =

b − a

( f(a) + 4 f((a+b)/2)+ f(b) ) [Regra de Simpson].

Por construção, a Regra dos Trapézios tem pelo menos grau 1, e como I(x²)≠ T(x²), tem exactamente grau 1.
A construção da Regra de Simpson garante pelo menos grau 2, mas pode verificar-se que I(x³)= T(x³), verificando-se que tem grau 3.

Caso não Linear.
No caso de não serem dados os nós de integração, estes são considerados como m+1 novas variáveis, e o sistema (P) é não linear.
Havendo 2m+2 graus de liberdade, podemos considerar 2m+2 equações e como tal a verificação de que a regra de quadratura terá grau 2m+1.
O caso mais simples será considerar apenas um nó z₀

Q(f) = w₀ f(z₀)

e vemos que é possível obter w₀ e z₀ de forma a que a regra tenha grau 1:

I(1) = Q(1) ⇔	b − a	= w₀
I(x) = Q(x) ⇔	(b² − a²)/2	= w₀ z₀

ou seja, obtemos w₀ = b − a e z₀ = (b + a)/2,

Q(f) = (b − a) f (

b + a

) [Regra do Ponto-Médio]

- O caso não linear será aprofundado mais tarde, no contexto das Fórmulas de Gauss.

- Começamos por considerar o caso linear, em que os nós são igualmente espaçados (denominadas fórmulas de Newton-Cotes). Em particular, iremos abordar os casos da Regra dos Trapézios e da Regra de Simpson, começando por analisar as regras simples e depois as regras compostas. Vamos distinguir duas noções

Regras Simples : são as regras deduzidas para um intervalo de referência.
Regras Compostas : são regras resultantes de dividir o intervalo inicial em sub-intervalos, aplicando as regras simples aos sub-intervalos.

Integração da Fórmula de Lagrange

Um processo de obter os pesos sem recorrer à resolução do sistema linear, é considerar a integração do polinómio interpolador. Dados os nós de integração z₀, ⋯, z_m, e sendo p_m o polinómio interpolador da função f, obtemos pela fórmula de Lagrange:

p_m(x) = f(z₀) L₀(x) + ⋯ + f(z_m) L_m(x)

portanto

I(p_m) = f(z₀) I(L₀) + ⋯ + f(z_m) I(L_m).

Para que a fórmula de quadratura seja exacta para um polinómio interpolador nestes nós, teremos então

w₀ = I(L₀) = ∫_[a,b]L₀(x)dx, ⋯, w_m = I(L_m) = ∫_[a,b]L_m(x)dx.

Exemplos:

Regra dos Trapézios. Os nós são z₀=a, z₁=b, e portanto usando o polinómio interpolador p₁ obtemos
T(f) = I(p₁) = f(a) ∫_[a,b] x-b
a-b dx + f(b) ∫_[a,b] x-a
b-a dx = f(a) (b-a)/2 + f(b) (b-a)/2
Regra de Simpson. Os nós são z₀=a, z₁=c, z₂=b, (com c=(a+b)/2), e portanto usando o polinómio interpolador p₂ obtemos
S(f) = I(p₂) = f(a) ∫_[a,b] x-c
a-c x-b
a-b dx + f(c) ∫_[a,b] x-a
c-a x-b
c-b dx + f(b) ∫_[a,b] x-a
b-a x-c
b-c dx

= f(a) (b-a)/6 + 4 f(c) (b-a)/6 + f(b) (b-a)/6

Convém ainda notar que como a Regra de Simpson tem grau 3, é ainda válida a igualdade

I(p₃) = S(p₃) = S(f)

para qualquer polinómio de grau ≤ 3 que interpole a função f nos nós a, c, b. Isto será utilizado na dedução da fórmula do erro na Regra de Simpson.

Iremos de seguida analisar o erro de integração, sendo necessário o seguinte resultado:

Teorema (valor intermédio para integrais)
Consideremos duas funções F ∈ C[a,b] e G integrável não mudando de sinal em [a,b]. Existe ξ∈[a,b]:

∫_[a,b] F(x) G(x) dx = F(ξ) ∫_[a,b] G(x) dx.

dem: Sem perda de generalidade, seja G≥ 0:

(min_[a,b] F) ∫_[a,b] G(x) dx ≤ ∫_[a,b] F(x) G(x) dx ≤ (max_[a,b] F) ∫_[a,b] G(x) dx.

Pelo Teorema do valor intermédio, para a função

H(x)= F(x) I(G), com I(G) = ∫_[a,b] G(x) dx

há pelo menos um ξ∈[a,b]:

H(ξ) = ∫_[a,b] F(x) G(x) dx.

Erro de Integração (Regras Simples)

O erro de integração é dado pela diferença entre o valor exacto e aproximado:

E(f) = I(f) − Q(f).

Tendo estabelecido uma relação com um polinómio interpolador, da forma

Q(f) = I(p_m)

em que m é o grau da fórmula, podemos obter:

E(f) = I(f) − I(p_m) = I(f − p_m)

reduzindo o erro de integração à integração do erro de interpolação.

Erro da Regra dos Trapézios (Simples).
Neste caso, os dois nós são a, b, e o polinómio interpolador é p₁. Relembrando o erro de interpolação:
f(x) − p₁(x) = f[a,b,x](x-a)(x-b)
obtemos
E_T(f) = ∫_[a,b] f(x) − p₁(x) dx = ∫_[a,b] f[a,b,x](x-a)(x-b) dx
e como G(x)=(x-a)(x-b) nunca é positiva em [a,b], obtemos pelo Teorema do valor intermédio para integrais:
∃ ξ ∈ [a,b] : E_T(f) = f[a,b,ξ] ∫_[a,b] (x-a)(x-b) dx
Quando f ∈ C²[a,b] podemos relacionar com a 2ª derivada, e como I(G)=(a-b)³/6 obtemos

E_T(f) = − f ''(ξ)/12 (b − a)³
Erro da Regra de Simpson (Simples).
Neste caso, há três nós a, c, b, e o polinómio interpolador é p₂. Relembrando o erro de interpolação:
f(x) − p₂(x) = f[a,c,b,x](x-a)(x-c)(x-b)
se consideramos G(x)=(x-a)(x-c)(x-b) vemos que há uma mudança de sinal no ponto médio c.
Por isso adoptamos um polinómio p₃ que interpola em a,c,b mas também num outro qualquer ponto d→ c
f(x) − p₃(x) = f[a,c,d,b,x](x-a)(x-c)(x-d)(x-b) → f[a,c,c,b,x](x-a)(x-c)²(x-b)
Assim, G(x)=(x-a)(x-c)²(x-b) não muda de sinal, tendo-se
E_S(f) = ∫_[a,b] f(x) − p₃(x) dx = f[a,c,c,b,ξ]∫_[a,b](x-a)(x-c)²(x-b) dx
pelo Teorema do valor intermédio para integrais (ξ ∈ [a,b]).
Quando f ∈ C⁴[a,b] relacionamos com a 4ª derivada, e de I(G)=(a-b)⁵/120 obtemos

E_S(f) = − f⁽⁴⁾(ξ)/2880 (b − a)⁵

Regras Compostas

Aplicamos agora as fórmulas anteriores a sub-intervalos, dividindo o intervalo original em N partes com um espaçamento igual

h = (b − a)/N, definindo os nós z_k = a + k h.

Regra dos Trapézios Composta
Notando que
∫_[a,b] = N
∑
k=1 ∫_{[z_k-1,z_k]}
aplicamos T_{[z_k-1,z_k]} (a R. Trapézios simples em cada sub-intervalo [z_k-1,z_k]) e definimos
T_N(f) = N
∑
k=1 T_{[z_k-1,z_k]} (f).
Como o espaçamento é igual podemos obter
T_N(f) = N
∑
k=1 (f(z_k-1)+f(z_k)) h/2
ou ainda

T_N(f) = h ( f(a)+f(b)
2 + N-1
∑
k=1 f(z_k) )
Regra de Simpson Composta
Seja N par. Subdividindo em intervalos de comprimento 2h
∫_[a,b] = N/2
∑
k=1 ∫_{[z_2k-2,z_2k]}
Aplicamos S_{[z_2k-2,z_2k]} (a R. Simpson simples em cada sub-intervalo [z_2k-2,z_2k]) e definimos
S_N(f) = N/2
∑
k=1 S_{[z_2k-2,z_2k]} (f)
Como o espaçamento é igual podemos obter
S_N(f) = N/2
∑
k=1 (f(z_2k-2)+4f(z_2k-1)+f(z_2k)) h/3
ou ainda

S_N(f) = h/3 ( f(a)+f(b) + 2 N/2 − 1
∑
k=1 f(z_2k) + 4 N/2
∑
k=1 f(z_2k-1) )

Notamos que a expressão de T₁ é a da regra dos trapézios simples, da mesma forma que a expressão de S₂ é a da regra de Simpson simples.

Erro das Regras Compostas

Como consequência do anterior, as fórmulas de erro para as regras compostas, resultam directamente duma soma das fórmulas de erro das regras simples.

Erro da Regra dos Trapézios Composta. Sendo f∈ C²[a,b] temos
E(T_N)(f) = I(f) − T_N(f) = N
∑
k=1 I_{[z_k-1,z_k]}(f) − N
∑
k=1 T_{[z_k-1,z_k]} (f)

= N
∑
k=1 I_{[z_k-1,z_k]}(f) − T_{[z_k-1,z_k]} (f) = − N
∑
k=1 f '' (ξ_k) h³/12
em que ξ_k∈[z_k-1,z_k] e como por aplicação do teorema do valor intermédio
∃ ξ ∈ [a,b] : f ''(ξ) = (1/N) N
∑
k=1 f '' (ξ_k),
e N h = b − a, obtemos

E(T_N)(f) = − f '' (ξ) h² (b − a)/12
Erro da Regra de Simpson Composta. Sendo f∈ C⁴[a,b] temos
E(T_N)(f) = I(f) − S_N(f) = N/2
∑
k=1 I_{[z_2k-2,z_2k]}(f) − N/2
∑
k=1 S_{[z_2k-2,z_2k]} (f)

= N/2
∑
k=1 I_{[z_2k-2,z_2k]}(f) − S_{[z_2k-2,z_2k]} (f) = − N/2
∑
k=1 f⁽⁴⁾(ξ_k) h⁵/90
em que ξ_k∈[z_2k-2,z_2k] e como por aplicação do teorema do valor intermédio
∃ ξ ∈ [a,b] : f⁽⁴⁾(ξ) = (2/N) N/2
∑
k=1 f⁽⁴⁾(ξ_k),
e ainda de N h = b − a, obtemos

E(S_N)(f) = − f⁽⁴⁾(ξ) h⁴ (b − a)/180

Fórmulas de Gauss-Legendre

Iremos agora verificar que é possível encontrar, de forma muito mais simples, os valores dos nós z₁, ..., z_n e dos pesos w₁, ..., w_n que permitem obter uma fórmula de quadratura de grau 2n-1, aproveitando os 2n graus de liberdade existentes.

Teorema:
Seja P_n o polinómio de Legendre de grau n, e z₁, ..., z_n as suas raízes. A fórmula de quadratura

Q(f) = w₁ f(z₁ )+ ... + w_n f(z_n ),

tem grau 2n-1, desde que os pesos w_k sejam obtidos pelo método dos coeficientes indeterminados ou pela integração dos polinómios de Lagrange

w_k = ∫_[-1,1] L_k(x) dx.

dem: Pretendemos ver que, sendo p_2n-1 um polinómio de grau ≤ 2n-1, temos

I(p_2n-1) = Q(p_2n-1).

Para isso consideramos a divisão de p_2n-1 por P_n, o que nos dá

p_2n-1 (x) = q_n-1(x) P_n(x) + r_n-1(x),

em que q_n-1 é o quociente e r_n-1 é o resto (ambos polinómios de grau ≤ n-1).
Como P_n(z_k) = 0, obtemos imediatamente Q(p_2n-1) = Q(r_n-1).
Por outro lado, como r_n-1 tem grau n-1, o polinómio interpolador em quaisquer n pontos coincide com r_n-1 (em particular, considerando os pontos z₁, ..., z_n). Assim, r_n-1 pode ser escrito de acordo com a fórmula de Lagrange

r_n-1(x) = r_n-1(z₁) L₁(x) + ... + r_n-1(z_n) L_n(x).

Portanto,

I(r_n-1) = r_n-1(z₁) I(L₁) + ... + r_n-1(z_n) I(L_n) = Q(r_n-1),

já que w_k = I(L_k).
Concluímos assim que Q(p_2n-1) = Q(r_n-1) = I(r_n-1).
Resta ver que I(r_n-1) = I(p_2n-1) para concluirmos o resultado. Ora,

I(p_2n-1) = I(q_n-1 P_n)+ I(r_n-1),

e como I(q_n-1 P_n) representa o produto interno entre q_n-1 e P_n, que são ortogonais, concluímos que I(q_n-1 Pn) = 0.
Falta apenas verificar que não é exacta para um polinómio de grau 2n. Para esse efeito, basta observar que

Q(P_n²) = 0 ≠ I(P_n²) = ∥P_n∥² .

Observação: Existem ainda outras Fórmulas de Gauss, aplicadas a integrais com peso W que define outros produtos internos.
Por exemplo, para o cálculo de

I(f) = ∫_[-1,1] W(x) f(x) dx

com W(x) = 1 / √ 1 − x² obtemos as fórmulas de Gauss-Chebyshev considerando os zeros dos polinómios ortogonais de Chebyshev (que são os nós de Chebyshev).
Este tipo de fórmulas revela-se especialmente útil para o cálculo deste tipo de integrais em que há uma parte singular (reparar que |W(x)|→∞ quando x = ± 1).

Observação − ordens de convergência

No capítulo sobre resolução de equações não lineares, encontrámos uma definição de ordem de convergência para sucessões que não deve ser confundida com a noção de ordem de convergência que aparece nas regras de integração compostas.

Assim, podemos dizer que uma regra de integração Q tem ordem de convergência p se verificar

E_h(f) = I(f) − Q_h(f) = O(h^p),

tendo-se verificado que isto ocorria para regras de grau p.
No entanto, à partida isto significa apenas que existe uma constante C>0:

|E_h(f)| ≤ C h^p, ou ainda formalmente |E_h(f)| ≈ C h^p.

A relação com as ordens de convergência das sucessões não é directa, pois o parâmetro h é livre e não está definida nenhuma sucessão.

Nota: A relação entre as duas noções é por vezes feita considerando uma sucessão de valores h,

h_n+1 = h_n/2 ⇒ h_n= 2^-n h₀

e dessa forma podemos considerar que

|E_{h_n}| ≈ C (h_n)^p = C h₀^p (2^-pn).

Perante uma sucessão destas a convergência será sempre linear com um coeficiente assimptótico K₁=2^-p.
É claro que perante uma outra sucessão de h, por exemplo h_n= 2^{-2^n} , obteríamos sempre convergência quadrática (para qualquer ordem p na convergência do método O(h^p))... e por outro lado, para h_n = 1/n^r teríamos sempre convergência logarítmica.

Convém por isso ter em atenção que há uma nomenclatura semelhante com significado diferente.

[Capítulo 1] [Capítulo 2] [Capítulo 3] [Capítulo 4] [Capítulo 5]

C J S Alves (2008, 2009)

Matemática Computacional Resumo da matéria [Versão 3.1]

Capítulo IV Aproximação de Funções

Interpolação

Interpolação geral

Interpolação Polinomial

Fórmula de Lagrange

Diferenças divididas

Fórmula de Newton

Erro de interpolação polinomial

Observação complementar − Nós de Chebyshev

Observação complementar − Interpolação trigonométrica

Método dos Mínimos Quadrados

Sistema Normal

Caso Discreto

Caso Contínuo

Integração Numérica

Método dos Coeficientes Indeterminados

Integração da Fórmula de Lagrange

Erro de Integração (Regras Simples)

Regras Compostas

Erro das Regras Compostas

Fórmulas de Gauss-Legendre

Observação − ordens de convergência

Matemática Computacional

Resumo da matéria [Versão 3.1]

Capítulo IV

Aproximação de Funções