Fenômeno de Runge

                Fala pessoal! Hoje nós vamos discutir um tópico mais avançado, que envolve Cálculo Numérico, basicamente interpolações polinomiais. Mas vamos abordar também em alguns pontos tópicos de computação. Recomendo a leitura deste artigo se você está no ensino superior e já estudou um pouco sobre cálculo numérico.

1) Introdução

                O fenômeno de Runge pode ser definido como um problema que aparece quando desejamos interpolar a função f(x)=1/(1+25x²) num intervalo I=[-1,1] por um polinômio de grau n em n+1 pontos eqüidistantes. A partir dos gráficos de P₄(x), P₈(x) e P₁₀(x) teremos uma ferramenta gráfica para verificar que, à medida que o grau do polinômio aumenta, mais ele se distancia da função perto das extremidades do intervalo I.

O teorema de de Weierstrass afirma que existe um polinômio P_n(x) que satisfaz:

Logo, seria normal pensar que o polinômio interpolador dos pontos eqüidistantes aproximaria-se melhor da função à medida que seu grau subisse. Mas esse pensamento é errado, porque o teorema nos garante que existe um polinômio apenas, ele não garante que o polinômio obtido por uma interpolação polinomial seja a solução. Além do mais, o fenômeno de Runge nos mostra que este tipo de pensamento é incorreto. O que polinômio da interpolação polinomial nos garante é que ele será igual a função nos n+1 pontos interpolados, apenas. Isso se assemelha a funções analíticas apenas em um intervalo, como por exemplo, a soma de uma série geométrica, que só converge se |x| for menor que 1 e diverge fora desse intervalo aberto.

                No nosso caso, temos que o polinômio que satisfaz o teorema de Weierstrass para o intervalo a=-1/5 e b=1/5 é o obtido por aproximação por séries. Com a soma de série geométrica:

mas fazendo a mudança de variável de x por -25u² e fazendo com que |-25u²|<1, ou seja |u|<0,2.  Interpolando a função por polinômios, chegaremos a mesma soma de série, que é polinomial. Mas ela é convergente apenas para |x|<0,2 logo, ao se aproximar das bordas -1 e 1, a aproximação polinomial também torna-se pior, haja vista que a série é divergente para esses x's. Outro modo de explicar o Fenômeno de Runge, seria a partir da prova de que f(x)-P_n(x), o erro,  tende ao infinito quando n tende ao infinito.

2) Método da Matriz

                O método da matriz consiste em elaborar uma matriz n+1 por n+1 do sistema linear que nos permite obter os coeficientes do polinômio P_n(x) que passa por todos os pontos. A matriz resulta das n+1 equações resultantes das n+1 igualdades P_n(xi)=f(xi).

A partir dessa matriz achamos os coeficientes a_n que multiplicam xⁿ e logo, o polinômio interpolados P_n(x). Neste exercício, para resolver este sistema linear vamos utilizar o método LU com pivoteamento parcial.

3) Decomposição LU com pivoteamento parcial

                A decomposição LU consiste em transformar uma matriz A numa multiplicação de matriz na forma LU, em que L é uma matriz nula acimada diagonal e U é nula abaixo da diagonal. Assim, o sistema Ax=b é mais facilmente resolvido, pois fazemos Ux=y e Ly=b e encontramos a solução imediatamente.

                A dedução para existência e unicidade das matrizes LU é feita utilizando as matrizes de escalanamento(M) da matriz A. A matriz U é obtida a partir da matriz A escalonada e a matriz L é obtida com as inversas em ordem contrária das matrizes M, multiplicadas. Um algoritmo que consome pouca memória é fazer o escalonamento com Gauss armazenando nas posições zeradas os multiplicadores. A matriz U será a matriz escalonada e a matriz L será nula acima da diagonal, 1 na diagonal e com os multiplicadores nas devidas posições, abaixo da diagonal.

                Para realizar a decomposição LU com pivoteamento parcial, basta armazenarmos um vetor de permutação P e realizar essa permutação depois no vetor b, bem como realizar as permutações dos multiplicadores na matriz escalonada por Gauss. Esse foi o método utilizado no algoritmo feito em Linguagem Python que será mostrado abaixo.

4) O algoritmo

                O algoritmo para obtenção das matrizes LU foi o descrito acima, sem realizar a inversão das matrizes M, para poupar o processamento.

#Funções Auxiliares

def permuta(matriz, linha1, linha2, b):

    """Recebe uma matriz e dois números referentes as linhas

e realiza a troca entre elas. Já realiza a troca em b, para não armazenar

o vetor permutação."""

    A=matriz

    blinha=b

    A[linha1]=matriz[linha2]

    A[linha2]=matriz[linha1]

    blinha[linha1]=b[linha2]

    blinha[linha2]=b[linha1]

    return None

def pivoteamento(matriz, coluna, b):

    """Recebe uma matriz e uma coluna e realiza o pivoteamento

da coluna, permutando as linhas da matriz e de b."""

    maximo=matriz[coluna][coluna]

    subs=0

    for linha in range(coluna+1,n):

        if abs(matriz[linha][coluna])>maximo:

               subs=linha

               maximo=abs(matriz[linha][coluna])

    permuta(matriz, coluna, subs,b)

    return None

def resolveU(U,b):

    """Recebe uma matriz U nula abaixo da diagonal e um vetor b e

resolve o sistema Ux=b, retornando o vetor x."""

    x=[]

    i=n

    soma=0

    while i>=0:

        for j in range(i+1,n):

            soma=soma+(A[i][j])*x[j]

        xi=(b[i]-soma)/A[i][i]

        x.insert(0,xi)

        i=i-1

    return (x)

def resolveL(L,b):

    """Recebe uma matriz L nula acima da diagonal e um vetor b e

resolve o sistema Lx=b, retornando o vetor x."""

    x=[]

    cont=0

    soma=0

    for i in range(0,n):

        for j in range(0,i-1):

            soma=soma+((A[i][j])*x[j])

        xi=(b[i]-soma)/A[i][i]

        x.append(xi)

    return (x)

#O programa

#1)Cálculo da Matriz

#1.1)Declaração de Variáveis

n=int(input("Digite o grau do polinômio interpolador Pn(x):"))

delta=(2/n)

xo=-1

cont1=cont2=0

A=linha=[]

#1.2)Obtenção da Matriz A

while cont1<=n:

    xo=-1+cont1*delta

    while cont2<=n:

        b=xo**cont2

        linha.append(b)

        cont2=cont2+1

    A.append(linha)

    cont1=cont1+1

#1.3)Obtenção de b

b=[]

cont=0

while cont<=n:

    x=-1+cont*(2/n)

    valor=(1/(1+(25*(x**2))))

    b.append(valor)

    cont=cont+1

#2)Obtenção de U e L

for coluna in range(0,n):

    pivoteamento(A,coluna,b)

    for linha in range(coluna, n):

        A[linha][cont]=(A[linha][cont])/A[cont][cont]

#2.1)Obtenção de U

U=A

for i in range(0,n):

    for j in range(0,n):

        if i>j:

            U[i][j]=0

#2.2)Obtenção de L:

L=A

for i in range(0,n):

    for j in range(0,n):

        if i==j:

            L[i][j]=1

        if j>i:

            L[i][j]=0

#3)Solução dos subsistemas

y=resolveL(L,b)

x=resolveU(U,y)

#4)O polinômio interpolador de grau n

print(x)

5) Resultados e Conclusão

                Os coeficientes dos polinômios foram obtidos executando o algoritmo no Shell do Python:

Plotamos os polinômios obtidos no Microsoft Mathematics, podendo observar claramente o Fenômeno de Runge: os polinômios afastaram-se das extremidades à medida que o grau aumentava. No gráfico abaixo, a curva azul indica a Função de Runge, a curva verde indica P₄(x), a curva rosa indica o P₈(x) e, por vim a curva preta indica o P₁₀(x).