Map0151

| Comunidade no Stoa | Disciplina no Moodle|

Números no Computador

Muitos dos algoritmos de cálculo numérico, por uma questão prática, deverão ser executados numa máquina real. Nestas máquinas, no entanto, a capacidade de memória para representação dos números é finita. Vários números reais (infinitos, de fato) terão a mesma representação no computador (ou calculadora), daí originando-se os erros de arredondamentos. Vamos ver qual é a técnica usada atualmente para diminuir os erros de arredondamentos.

Representação de números inteiros numa base

Seja $\beta > 1$ um número natural. Então podemos representar qualquer número inteiro $k$ como a soma

$$k= \text{sgn}(k)*a_0\beta^0 + \cdots + a_s \beta^s$$ onde cada algarismo $a_j$ representa um número natural entre $0$ e $\beta -1$. Só $a_s$ é diferente de zero. Esta representação é única. A representação de $k$ na base $\beta$ é $$k=\pm [a_s\dots a_0]_{\beta}$$ Deixando o sinal de $k$ de fora, para simplificar. Podemos calcular os algarismos do número daddo na base $\beta$ aplicando repetidas vezes o algoritmo da divisão. Escrevendo $$k = a_0 + \beta(a_1 +\beta( a_2 + \cdots +\beta a_s)\cdots ))$$ vemos que $a_0$ é o resto da divisão $k//\beta$, $a_1$ é o resto da divisão $(k-a_0)//\beta$ e assim por diante.

Exemplo: $39 = [100111]_2$

Representação de números fracionários e decimais numa base

Seja novamente $\beta > 1$ um número natural.

Se $x \in (0,1)$ então

$$x= \frac{b_1}{\beta} + \cdots + \frac{b_k}{\beta^k} + \cdots$$ Diremos que $x=[0.b_1 b_2\dots]_{\beta}$ é a representação fracionária na base $\beta$ de $x$. Quando não houver dúvidas de que base se trata, omite-se a base da notação!

Exemplo: $0.9= [0.1110011001100...]_2$

Uma função para colocar um número decimal na forma binária

In[1]:

# -*- coding: utf-8 -*- """ Spyder Editor """ def binario(a):     # da a representacao binaria do numero a     ParteInteira = int(a)     ParteDecimal = a-int(a)     # representacao binaria da parte inteira     # a lista seguinte guarda os dígitos da parte inteira     ListaDigitos=[]     while (ParteInteira > 0):         ListaDigitos.append(ParteInteira%2)         ParteInteira=ParteInteira//2     # lista dos digitos depois da virgula     ListaResto=[]     k=1     while ((ParteDecimal!=0)&(k<50)):         ListaResto.append(int(2*ParteDecimal))         ParteDecimal=2*ParteDecimal - int(2*ParteDecimal)         k=k+1     # produz a string de representacao:             i=len(ListaDigitos)-1     p1=""     while (i>=0):         p1=p1+str(ListaDigitos[i])         i=i-1     # Depois disso p1 tem a parte inteira     l=0     p2=""     while (l<len(ListaResto)):         p2=p2+str(ListaResto[l])         l=l+1         return p1+"."+p2     print (binario(21.75))

Representação em ponto flutuante

Consideramos uma base fixa $\beta > 1$. Um número real $\alpha \in \mathbb{R}$ positivo pode ser escrito nesta base como:

$$\alpha = [a_k\cdots a_0.b_1b_2\cdots]_{\beta}$$ Isto significa que: $$\alpha = a_k\beta^k + \cdots + a_0 + \frac{b_1}{\beta} +\frac{b_2}{\beta^2} + \cdots$$ Na equação acima, colocando $\beta^{k+1}$ em evidência temos: $$\alpha = \left( \frac{a_k}{\beta} + \cdots + \frac{a_0}{\beta^{k+1}} + \frac{b_1}{\beta^{k+2}} +\frac{b_2}{\beta^{k+3}} + \cdots\right)\times \beta^{k+1}$$ ou ainda, lembrando da notação de um número numa base dada: $$\alpha = [0.a_k \cdots a_0 b_1 b_2\cdots]_{\beta} \times \beta^{k+1}$$ Esta última fórmula é importante. O número real $\alpha$ fica caracterizado por três dados:

• O número $m = [0.a_k \cdots a_0 b_1 b_2\cdots]_{\beta} \in [0.1,1)$ chamado de mantissa.
• O número $e = k+1$ chamado de expoente
• O sinal do número $\sigma$

Esta representação do número $\alpha$ como $\sigma m \times \beta^e$ chamaremos de representação normal em ponto flutuante na base $\beta$. Em geral a base fica clara pelo contexto!

Números de máquina

Continuamos com a base fixa ($\beta$), mas na representação normal em ponto flutuante vamos admitir apenas números com a mantissa limitada a $D$ dígitos e o expoente limitado, em módulo, por um número inteiro $M$. O conjunto

$$\mathbf{M} = \left\{ [0.d_1\cdots d_D]_{\beta}\times \beta^e: |e| \leq M \right\}$$ é um conjunto de números de máquina definido pelos números $(\beta, D, M)$. Este é um conjunto finito com $(\beta-1)(\beta^{(D-1)})(2M+1)$ elementos. Vamos agora definir a função arredondamento, que é uma função sobrejetora $\operatorname{rd}: \mathbb{R} \to \mathbf{M}$. Seja $\alpha$ um número real; se o expoente de $\alpha$ for maior que $M$ não representamos o número. Senão tomamos sua mantissa e se o dígito $b_{d+1}\geq \beta/2$ escolhemos o menor número de máquina maior que $\alpha$ para representá-lo. No caso de $b_{d+1} <\beta/2$ o representante será o maior número de máquina menor que $\alpha$.

Se denotarmos por $\underline{\alpha}$ o truncamento do número $\alpha$, isto é, o maior número de máquina menor que $\alpha$. E por $\overline{\alpha}$ o menor número de máquina maior que $\alpha$. Temos que:

$$\overline{\alpha} - \underline{\alpha} = \beta^{e-D}$$ Com isto temos uma avaliação do erro absoluto e erro relativo do arredondamento: $$|\alpha -\operatorname{rd}(\alpha)| \leq \frac{\beta^{e-D}}{2}$$ $$\frac{|\alpha -\operatorname{rd}(\alpha)|}{|\alpha|}\leq \frac{\beta^{1-D}}{2}$$

Aritmética de números em pontos flutuantes

A extensão das operações aritméticas de soma e multiplicação de números reais faz-se pelo arredondamento:

$$\alpha \oplus \beta = \operatorname{rd}( \operatorname{rd}(\alpha) + \operatorname{rd}(\beta))$$ $$\alpha \otimes \beta = \operatorname{rd}( \operatorname{rd}(\alpha) \times \operatorname{rd}(\beta))$$ Estas operações deixam de ser associativas.

Aproximações usando polinômios de Taylor

Para fazer as avaliações de funções de variável real, utilizaremos apenas as operações aritméticas. Basicamente, isto significa que só usaremos polinômios. Em primeiro lugar veremos o método de avaliação de Horner.

Método de Horner

É a forma de se avaliar o polinômio num ponto dado executando o menor número de operações aritméticas. O polinômio da forma

$$p(x) =a_nx^n+ a_{n-1}x^{n-1}+\cdots + a_0$$ no ponto $x_0$ é avaliado na ordem $$p(x)= [[[[a_nx_0 + a_{n-1}]x_0+a_{n-2}]x_0+\cdots]x_0 +a_0]$$ Isto é, calculamos a sequência

• $A_0 = a_n$
• $A_{k+1} = A_kx_0+ a_{n-k}$ até o $A_n$ que é a avaliação.

Teorema de Taylor

Teorema: Seja $f:(a,b) \to \mathbb{R}$ uma função de classe $\mathcal{C}^r$, talque exista a derivada $f^{(r+1)}(c), \forall c \in (a,b)$. Então para todo $x$ do intervalo, e $c$ fixado existe $\xi \in |c,x|$ tal que:

$$f(x) = f(c) + f^\prime(c)(x-c)+ \cdots + \frac{f^{(r)}(c)}{r!}(x-c)^{r} + \frac{f^{(r+1)}(\xi)}{(r+1)!}(x-c)^{r+1}$$

Zeros de funções

Dada uma função $f:(a,b)\to \mathbb{R}$, queremos encontrar um $r\in (a,b)$ tal que $f(r)=0$. Um tal número será chamado de zero da função $f$. Para assegurar a existência de zeros no intervalo, precisamos de alguma regularidade da função. Se $f$ é uma função contínua, uma condição suficiente para a existência do zero é que $f(a)*f(b) < 0$, e uma técnica para chegar ao zero é ir contraindo este intervalo, mantento a relação $f(a)*f(b) < 0$, até chegarmos suficientemente próximo. No método da bissecção vamos dividindo o intervalo ao meio. Se $(a_n,b_n)$ é o intervalo com $f(a_n)*f(b_n) < 0$, definimos $x_{n+1} = (a_n+b_n)/2$ e o novo intervalo como:

• $(a_{n+1},b_{n+1})=(x_{n+1},b_n)$ se $f(a_n)*f(x_{n+1}) < 0$
• $(a_{n+1},b_{n+1})=(a_n,x_{n+1})$ caso contrário

O tamanho do intervalo no n-ésimo passo é $|b-a|/2^n$. A seguir uma função em python que implementa a bissecção.

Algoritmo da bissecção

In[2]:

#define a funcao bisseccao def bissec(f, a, b, eps=0.0001):     ''' Uma funcao que executa a biseccao do intervalo (a,b), não testa se tem zero '''     x0=a     x1=b     k=0     print("|  k  | x0       | x1       | xmid     |f(x0)*f(xmid)|")     while (x1-x0)/2 > eps :         xmid = (x1+x0)/2         # incluir print para debug e o contador         print("| {0:<2}  | {1:1.6f} | {2:<1.6f} | {3:<1.6f} | {4:<+2.8f} |".format(k,x0,x1,xmid,f(xmid)*f(x0)))         if f(xmid)*f(x0)>0:             x0,x1=xmid,x1         else : x0,x1 = x0,xmid         k+=1     return (x0+x1)/2     #Rodando um teste f = lambda x: x**2-2 a,b = 0,2 bissec(f,a,b)

Out[2]:

| k | x0 | x1 | xmid |f(x0)*f(xmid)| | 0 | 0.000000 | 2.000000 | 1.000000 | +2.00000000 | | 1 | 1.000000 | 2.000000 | 1.500000 | -0.25000000 | | 2 | 1.000000 | 1.500000 | 1.250000 | +0.43750000 | | 3 | 1.250000 | 1.500000 | 1.375000 | +0.04785156 | | 4 | 1.375000 | 1.500000 | 1.437500 | -0.00726318 | | 5 | 1.375000 | 1.437500 | 1.406250 | +0.00245667 | | 6 | 1.406250 | 1.437500 | 1.421875 | -0.00048804 | | 7 | 1.406250 | 1.421875 | 1.414062 | +0.00000960 | | 8 | 1.414062 | 1.421875 | 1.417969 | -0.00000454 | | 9 | 1.414062 | 1.417969 | 1.416016 | -0.00000218 | | 10 | 1.414062 | 1.416016 | 1.415039 | -0.00000100 | | 11 | 1.414062 | 1.415039 | 1.414551 | -0.00000041 | | 12 | 1.414062 | 1.414551 | 1.414307 | -0.00000011 | | 13 | 1.414062 | 1.414307 | 1.414185 | +0.00000004 |

Algoritmo de Newton

O algoritmo de Newton-Raphson usa propriedades de diferenciabilidade da função $f$ para criar uma sequência recorrente que converge ao zero da função que atrairá um ponto $x_0$, que é uma estimativa inicial. A fórmula para gerar a sequêmcia é:

$$x_{k+1} = x_k - \frac{f(x_k)}{f^\prime(x_k)}$$ Abaixo o algoritmo em python:

In[3]:

# ilustra o processo de newton-raphson def newtonRaphson(f,flinha,x0,eps=0.0001):     '''executa o procedimento de newton'''     a=x0     b= a - (f(a)/flinha(a))     k=1 #quero contar as iterações     print("| k | x_k       | |x_k-xk-1|   |")     print("|{0:<2} | {1:<2.7f} | {2:<2.7f}".format(k,b,abs(b-a)))       while (abs(b-a)>eps):         a=b         b= a - (f(a)/flinha(a))         k+=1         print("|{0:<2} | {1:<2.7f} | {2:<2.7f}".format(k,b,abs(b-a)))             return b # função de teste f = lambda x: x**2-2 flinha = lambda x: 2*x x0=1.5 newtonRaphson(f,flinha,x0)

Out[3]:

| k | x_k | |x_k-xk-1| | |1 | 1.4166667 | 0.0833333 |2 | 1.4142157 | 0.0024510 |3 | 1.4142136 | 0.0000021

Algoritmo da secante

O método da secante pode ser visto como uma adaptação do método de Newton, onde a derivada da função $f$ é computada pela fórmula de aproximação

$$f^\prime(x_k) \sim \frac{f(x_k)-f(x_{k-1})}{(x_k - x_{k-1})}$$ definindo a sequência $$x_{k+1} = x_k - \frac{f(x_k)(x_k -x_{k-1})}{(f(x_k)-f(x_{k-1}))}$$

O código em Python pode ser o seguinte

In[4]:

# Uma nova definição da secante usando o pseudo-código Kincaid def secante(f,x0,x1,eps=0.0001):         '''Método da secante'''     a=x0     b=x1     u = f(a)     v = f(b)     if abs(v) > abs(u):         a,b=b,a         u,v=v,u     c = b - v*(b-a)/(v-u)     k=1     print("| k | x_k       | f(xk)   |")     print("|{0:<2} | {1:<2.7f} | {2:<2.7f}".format(k,a,u))     print("|{0:<2} | {1:<2.7f} | {2:<2.7f}".format(k,b,v))     while abs(c-b) > eps:         a=b         b=c         u=f(a)         v=f(b)         c = b - v*(b-a)/(v-u)         k+=1         print("|{0:<2} | {1:<2.7f} | {2:<2.7f}".format(k,b,v))     return c #Teste f= lambda x: x**2 - 2 secante(f,1,2)

Out[3]:

| k | x_k | f(xk) | |1 | 2.0000000 | 2.0000000 |1 | 1.0000000 | -1.0000000 |2 | 1.3333333 | -0.2222222 |3 | 1.4285714 | 0.0408163 |4 | 1.4137931 | -0.0011891 |5 | 1.4142114 | -0.0000060

Resolução de sistemas lineares

Do ponto de vista de algoritmos numéricos existem duas formas para se resolver um sistema de equações lineares da forma:

$$\begin{align} a_{11}x_1 + \cdots + a_{1n}x_n &=& b_1 \\ \vdots & & \\ a_{n1}x_1 + \cdots +a_{nn}x_n &= &b_n \end{align}$$

Os Métodos Diretos têm, em geral, duas fases. Numa primeira fase, um sistema genérico é transformado, numa sequência finita de passos, num sistema de uma classe de sistemas, do qual se conhece um algoritmo eficiente para a solução. Na segunda fase, resolvemos o sistema gerado. O método da eliminação de Gauss (MEG) é um método direto que descreveremos a seguir.

Os Métodos Iterativos são obtidos por uma sequência de estimativas da solução. Esta sequência deve convergir para a solução, e muito raramente se obtem a solução exata. Mas tem-se um bom controle de erros. Descreveremos o método de Gauss-Seidel que ilustra esta técnica.

Método da Eliminação de Gauss

Em primeiro lugar vamos descrever os algoritmos para a resolução de sistemas triangulares. Uma matriz quadrada $U = (u_{ij})$ é triangular superior quando tivermos a condição $u_{ij} = 0$ para $i> j$. Neste caso podemos resolver o sistema linear $Ux=b$,usando o algoritmo de baixo pra cima em $n$ passos:

• Passo 1 : Encontro $x_n = b_n/u_{nn}$
• Passo k : Encontro $x_{n-k+1} = \frac{1}{u_{n-k+1 n-k+1}}\left(b_{n-k+1} - \sum_{j=n-k+2}^{n} u_{n-k+1 j}x_j\right)$

A matriz $L=(l_{ij})$ é triangular inferior quando $l_{ij}=0 \text{ e } l_{ii}=1$ para $i < j$. Para resolver os sistema $Lx=b$ usando um algoritmo parecido com o anterior.

• Passo 1 : Encontro $x_1 = b_1/l_{11}$
• Passo k : Encontro $x_{k} = \frac{1}{u_{kk}}\left(b_{k} - \sum_{j=1}^{k-1} u_{k j}x_j\right)$

No método da eliminação de Gauss (MEG), vamos transformar um sistema linear numa sistema equivalente triangular superior, isto é, quando a matriz dos coeficientes $A$ é triangular superior. Faremos isso usando as operações elementares, que descrevemos a seguir:

• $\mathcal{O}(i,j)$ é a troca da linha $i$ com a linha $j$
• $\mathcal{O}(i,\alpha)$ é a multiplicação da linha $i$ por um número, não nulo, $\alpha$
• $\mathcal{O}(i,\alpha, j)$ é a substituição da linha $i$ pela soma desta com $\alpha$ multiplicado pela $i-$ésima linha.

O algoritmo é a aplicação sucessivas dos seguintes passos com $k$ variando de $1$ até $n-1$

 Passo $k$: Caso $a_{kk} \neq 0$, faça de $i=k+1$ até $i=n$ $$m_{ik} = a_{ik}/a_{kk}$$  e aplique ao sistema a operação elementar $\mathcal{O}(i,m_{ik}, k)$.

Ao fim do passo $n-1$ o sistema estará triangularizado. Para evitar a possibilidade de não se triangularizar o sistema, mesmo ele sendo possível e determinado, executamos a pivotação que é o seguinte procedimento: antes da execução do passo $k$ do MEG. Buscamos, na $k-$ésima coluna, a partir da $k-$ésima linha o índice $i$ tal que $|a_{ik}|$ seja o maior possível, e executamos $\mathcal{O}(i,k)$. Com isto garantimos a triangularização de qualquer sistema que tenha solução única.

A seguir o código em python para a eliminação de Gauss

In[5]:

# gauss elimination as seen in the classroom import numpy as np # A is type array from numpy # remember: A[0.0] is the first element # Elementary row operations: def  troca_linha(A,i,j):     '''Troca as linhas i e j da matriz A'''     buffer = A[i].copy()     A[i] = A[j]     A[j] = buffer     return A def mult_linha(A,i,alfa):     '''Multiplica por alfa a linha i'''     A[i] = alfa*A[i]     return A def subs_linha(A,i,k,alfa):     '''soma a linha i com um multiplo da linha k'''     A[i] = A[i] + alfa*A[k]     return A # step k of gauss elimination: def gauss_step(A,k):     '''O k-esimo passo na eliminacao de gauss'''     L=range(len(A))     for j in L[k+1:]:         m=-A[j,k]/A[k,k]         A=subs_linha(A,j,k,m)         A[j,k]=-m # Aproveitamos a matriz A pra guardar o multiplicador     return A # triangle form, without care and pivoting def triangle_form(A):     '''Retorna a matriz escalonada,     com os multiplicadores na parte triangular inferior'''     L=len(A)     for k in range(L-1):         p=pivo(A,k)         troca_linha(A,k,p)         gauss_step(A,k)           return A     def pivo(A,n):     ''' retorna o indice de pivotacao da matriz A na coluna n,     a partir da linha n'''     l=len(A) # numero de linhas de A     c=len(A[0]) # numero de colunas     if (l<=n) or (c<=n):         raise ValueError("n deve ser menor que dimensao de A")     p=n # inicio o pivo como n, e vou aumentando     for k in range(n,l):         if abs(A[k,n]) > abs(A[p,n]):             p=k     return p # roda um teste A=np.array([[1,2,3],          [2,-1,9],          [3.,0., 2.]]) B=triangle_form(A)

Versão sem pivotação (não sei se precisa!)

In[6]:

# gauss elimination as seen in the classroom import numpy as np # A is type array from numpy # remember: A[0.0] is the first element # Elementary row operations: def  troca_linha(A,i,j):     buffer = A[i]     A[i] = A[j]     A[j] = buffer     return A def mult_linha(A,i,alfa):     A[i] = alfa*A[i]     return A def subs_linha(A,i,k,alfa):     A[i] = A[i] + alfa*A[k]     return A # step k of gauss elimination: def gauss_step(A,k):     L=range(len(A))     for j in L[k+1:]:         m=-A[j,k]/A[k,k]         A=subs_linha(A,j,k,m)     return A # triangle form, without care and pivoting def triangle_form(A):     L=len(A)     for k in range(L):         if A[k,k] != 0:             A=gauss_step(A,k)             print A # for testing         else: print ' falhou no passo ', k     return A # roda um teste A=np.array([[1,2,3],          [2,-1,9],          [3.,0., 2.]]) B=triangle_form(A)

Decomposição LU

Diremos que uma matriz $L$ é triangular inferior se $L=(l_{ij})$ com

$$l_{ij}=0 \text{ se } i \lt j \\ l_{ii}=1$$

Por outro lado, a matriz $U$ será triangular superior se $U=(u_{ij})$ com

$$u_{ij}=0 \text{ se } j \lt i$$ Uma matriz quadrada $A$ tem uma decomposição LU se existem uma matriz triangular inferior, $L$ e uma triangular superior $U$ tal que temos a igualdade: $$LU=A$$

Uma matriz inversível tem decomposição LU se e somente se seus menores principais têm determinantes não nulos. Neste caso pode-se resolver o sistema

$$Ax=b$$ Resolvendo primeiro $$Lz=b$$ usando substituição de cima pra baixo, e depois, resolver $$Ux =z$$ usando substituição de baixo pra cima.

Métodos Iterativos

A idéia dos métodos iterativos é produzir uma sequência de estimativas vetoriais $\mathbf{x}^{(k)}$, que idealmente convirja para a solução $\mathbf{x}$ e cada passo é calculado a partir da estimativa anterior. Costumam ser bem mais rápidos no caso de matrizes muito grandes. Para construir recursivamente a sequência que dá origem ao método iterativo fazemos uma decomposição da matriz $A$ de coeficientes do sistema como

$$A = N - P$$ Onde $N$ é de uma classe de sistemas inversíveis e fáceis de resolver. Então resolver o sistema: $$A\mathbf{x} =b$$ é equivalente a resolver a equação: $$N\mathbf{x} = b+ P\mathbf{x} \text{ ou } \mathbf{x}=N^{-1}(b+ P\mathbf{x})$$ A recursão $$N\mathbf{x}^{(k+1)} = b+ P\mathbf{x}^{(k)}$$ é a que converge para a solução caso $||N^{-1}P|| < 1$ onde a norma de matrizes é induzida de uma norma em $\mathbb{R}^n$. Destacamos dois casos

• Método de Jacobi $$N = \text{diag}(a_{11}, \cdots, a_{nn})$$
• Método de Gauss-Seidel $$N = (n_{ij}) \text{ com } n_{ij} = a_{ij} \text{ se } i\geq j \\ n_{ij}=0 \text{ se } i\lt j$$

Método dos Mínimos Quadrados: Caso discreto

Consideremos uma tabela de pontos $T = \{(x_0,y_0),\cdots,(x_k,y_k)\}$. Podemos interpretar esta tabela de duas formas: 1. Pode ser um conjunto de informações aproximadas sobre uma função de determinado tipo. E neste caso procuraremos a função daquele tipo que melhor se ajuste à tabela. 2. Pode ser um conjunto de informações precisas sobre uma função diferenciável. Mas não sabemos mais o tipo da função.

No primeiro caso faremos um ajuste da função para a tabela pelo Método dos Mínimos quadrados. A função ajustada não respeita, necessariamente, a tabela dada. Para o segundo caso estudaremos a Interpolação polinomial

Comecemos com o método dos mínimos quadrados (MMQ). Lembremos que temos dadas uma tabela de pontos $T = \{(x_0,y_0),\cdots,(x_k,y_k)\}$ e uma família de funções $\mathcal{F}$ onde o domínio de cada uma das funções de $\mathcal{F}$ contém os pontos $x_i$ da tabela $T$. O resíduo quadrático de uma função $f$ com relação à tabela é:

$$r(f,T)= \sum_{i=0}^k \left(y_i -f(x_i)\right)^2$$ O problema de MMQ é encontrar uma função $f^* \in \mathcal{F}$ que minimize $r(f,T)$