Notas de Aula Matematica Computacional, Notas de aula de Matemática Aplicada

Baixe Notas de Aula Matematica Computacional e outras Notas de aula em PDF para Matemática Aplicada, somente na Docsity!

Matemática Computacional

Adérito Luís Martins Araújo

Notas de apoio às aulas de Matemática Computacional do Mestrado

Integrado em Engenharia Electrotécnica e de Computadores, no ano

lectivo de 2013/2014.

Figura da capa: Documento jurídico sumério YBC 3879, do período Ur III (c. 2100 – c. 2000 a.C.), onde é descrito um algoritmo geométrico para determinar as soluções de equações quadráticas. Fonte: cdli.ucla.edu/pubs/cdlj/2009/cdlj2009_003.html.

Aritmética computacional 2

tempo; o facto de não ser possível realizar os cálculos podia passar de mero detalhe técnico. Afinal de contas, o problema dos dois corpos (isto é, dois corpos que interagem por via da força gravitacional, como a Terra e o Sol) tinha uma solução muito simples, que era estudada no primeiro ano das universidades. O facto é que a solução analítica deste problema é impossível de obter! Resta-nos assim recorrer à solução numérica.

O estabelecimento das várias leis da física permitiu aos matemáticos e aos físicos obter modelos para a mecânica dos sólidos e dos fluidos. As engenharias mecânica e civil usam esses modelos como sendo a base para os mais modernos trabalhos em dinânica dos fluidos e em estruturas sólidas, e a análise numérica tornou-se uma ferramenta essencial para todos aqueles que pretendem efectuar investigação nessas áreas da engenharia. Por exemplo, a construção de estruturas modernas faz uso do chamado método dos elementos finitos para resolver as equações com derivadas parciais associadas ao modelo; a dinâmica dos fluidos computacional é actualmente uma ferramenta fundamental para, por exemplo, desenhar aviões; a elaboração de novos materiais é outro assunto que recorre, de forma intensa, a algoritmos numéricos. A análise numérica é pois uma área que tem assumido crescente importância no contexto das ciências da engenharia.

No processo de resolução de um problema físico podemos distinguir várias fases.

1. Formulação de um modelo matemático que descreve uma situação real. Tal formula- ção pode ser feita recorrendo a (sistemas de) equações algébricas, transcendentes, integrais, diferenciais, etc. É necessário ter muito cuidado nesta fase uma vez que a grande complexidade dos problemas físicos pode-nos obrigar a fazer simplificações no modelo, simplificações essas que não devem alterar grandemente o comportamento da solução.
2. Obtenção de um método numérico que permite construir uma solução aproximada para o problema. Um método numérico que possa ser usado para resolver o problema é traduzido por algoritmo que não é mais do que um completo e não ambíguo conjunto de passos que conduzem à solução do problema. Esta fase constitui o cerne da análise numérica. Dado um determinado método numérico, temos necessidade de saber em que condições as soluções por ele obtidas convergem para a solução exacta; em que medida pequenos erros de arredondadmento (e outros) poderão afectar a solução final; qual o grau de precisão da solução aproximada obtida, etc.
3. Programação automática do algoritmo. Nesta fase teremos necessidade de recorrer a uma linguagem de programação como o Fortran, o Pascal, o C++, entre outras. Mais recentemente é usual o recurso a programas como o Mathematica ou o Matlab.

Os algoritmos numéricos são quase tão antigos quanto a civilização humana. Os babi- lónios, vinte séculos antes de Cristo, já possuiam tabelas de quadrados de todos os inteiros entre 1 e 60. Os egípcios, que já usavam fracções, inventaram o chamado método da falsa posição para aproximar as raízes de uma equação. Esse método encontra-se descrito no papiro de Rhind, cerca de 1650 anos antes da era cristã. Na Grécia antiga muitos foram os matemáticos que deram contributos para o impulso desta disciplina. Por exemplo, Arquimedes de Siracusa (278-212, a.C.) mostrou que

3

< π < 3

e apresentou o chamado método da exaustão para calcular comprimentos, áreas e volumes de figuras geométricas. Este método, quando usado como método para calcular aproximações,

Aritmética computacional 3

está muito próximo do que hoje se faz em análise numérica; por outro lado, foi também um importante precursor do desenvolvimento do cálculo integral por Newton e Gottfried Wilhelm von Leibniz (1646-1716).

√^ Heron de Alexandria^ (∼10-∼75), no século I, deduziu um procedimento para determinar a da forma (como deduzir este método?)

x(n+1)^ =

x(n)^ +

a x(n)

No ano 250, Diofanto de Alexandria (∼200-∼284) obteve um processo para a determinação das soluções de uma equação quadrática. Durante a Idade Média, os grandes contributos para o desenvolvimento da matemática algorítmica vieram, sobretudo, do médio oriente, Índia e China. O contributo maior foi, sem dúvida, a simplificação introduzida com a chamada numeração indo-árabe. O aparecimento do cálculo e a criação dos logaritmos, no século XVII, vieram dar um grande impulso ao desenvolvimento de procedimentos numéricos. Os novos modelos mate- máticos propostos não podiam ser resolvidos de forma explícita e assim tornava-se imperioso o desenvolvimento de métodos numéricos para obter soluções aproximadas. O próprio Newton criou vários métodos numéricos para a resolução de muitos problemas, métodos esses que possuem, hoje, o seu nome. Tal como Newton, muitos vultos da matemática dos séculos XVIII e XIX trabalharam na construção de métodos numéricos. De entre eles podemos destacar Leonhard Euler (1707-1783), Joseph-Louis Lagrange (1736-1813) e Johann Carl Friedrich Gauss (1777-1875). Foi, no entanto, o aparecimento, na década de 40 do século XX, dos primeiros compu- tadores que contribuiu decisivamente para o forte desenvolvimento da disciplina. Apesar de tanto Blaise Pascal (1623-1662) como Leibniz terem construído, já no séc. XVII, as primeiras máquinas de calcular e de Charles Babbage (1791-1871), milionário inglês, ter construído o que é considerado o primeiro computador (nunca funcionou!), foi apenas com o aparecimento do ENIAC, nos anos 40, que a ciência usufruiu, de facto, desses dispositivos de cálculo.

1.1 Erros absolutos e relativos

O processo de resolução numérica de um problema passa sempre pela substituição de um problema difícil para um problema mais simples que tenha a mesma solução ou, pelo me- nos, uma solução que se relaciona intimamente com a do problema original. As soluções numéricas são, por isso, aproximações dos resultados pretendidos. Os motivos para con- siderar essas aproximação podem ser de vária ordem e a precisão da solução final reflecte todo esse processo. Também há que ter em conta a incerteza nos dados iniciais, assim como as perturbações efectuadas ao longo dos cálculos, que podem ser amplificadas pelo algoritmo.

Exemplo 1.2 Pretende-se calcular a superície terrestre pela fórmula A = 4πr^2. Este processo involve uma aproximação. Primeiro, a Terra não é uma esfera. A esfera é uma idealização da sua verdadeira forma. Depois, o valor do raio da Terra é baseado em medições empíricas e cálculos anteriores. O valor de π requer a truncatura de um processo infinito. Finalmente, os valores dos dados iniciais e os resultados das operações aritméticas são arredondados durante os cálculos no computador.

O exemplo anterior mostra que a introdução de erros num determinado processo de cálculo pode ter várias causas. É nosso objectivo analisar quais são essas causas e estudar

Aritmética computacional 5

Definição 1.3 (Erro absoluto) Seja x ∈ Rn^ um vector cujas componentes são desco- nhecidas e x ∈ Rn^ um vector cujas componentes são aproximações para as componentes correspondentes de x. Chama-se erro absoluto de x à quantidade ‖e(x)‖.

Definição 1.4 (Erro relativo) Seja x ∈ Rn, x 6 = 0, um vector cujas componentes são desconhecidas e x ∈ Rn^ um vector cujas componentes são aproximações para as compo- nentes correspondentes de x. Chama-se erro relativo de x, e representa-se por r(x), à quantidade r(x) = ‖e(x)‖/‖x‖.

Como na definição de erro relativo o valor de x não é conhecido, é usual considerar a aproximação r(x) ≈ ‖e(x)‖/‖x‖. Melhor ainda, atendendo a que

‖x‖ ≥ |‖x‖ − ‖e(x)‖|,

podemos considerar o majorante

r(x) ≤ ‖e(x)‖ |‖x‖ − ‖e(x)‖|

O erro relativo, atendendo a que é uma quantidade adimensionada, é muitas vezes representado sob a forma de percentagem. Note-se também que o erro relativo nos dá uma maior informação quanto à precisão da aproximação que o erro absoluto. É com base nas definições de erro absoluto e erro relativo que iremos analisar os re- sultados numéricos que aparecerão como aproximações a valores que não conhecemos com exactidão. Suponhamos que pretendemos calcular o valor de uma função F : Rn^ −→ Rm para um dado argumento vectorial x ∈ Rn^ (valor exacto). O valor desejado é, pois, F (x). No entanto, o que realmente é calculado é o valor de F (x), onde x é o vector das aproxi- mações e F^ é a aproximação à função, o^ algoritmo. Assim, o erro total é dado por

F (x) − F (x) = F (x) − F (x) ︸ ︷︷ ︸ erro computacional

• F (x) − F (x) ︸ ︷︷ ︸ erro de propagação

O algoritmo não tem qualquer influência nos erros de propagação. Apenas influencia os erros computacionais, que são a soma dos erros de arredondamento com os erros de truncatura.

1.2 Erros de arredondamento e truncatura

Os dados de um determinado problema podem estar à partida afectados de imprecisões resultantes de medições incorrectas. Note-se que a escala de um instrumento de medição nos dá uma possibilidade de saber um limite superior para o erro com que vêm afectados os valores medidos. Por exemplo, com uma régua usual, a medição de uma distância de 2 mm pode vir afectada com um erro de 0 , 5 mm o que dá um erro relativo de 25%. Outra causa de erro resulta das simplificações impostas ao modelo matemático usado para descrever um determinado fenómeno físico. Por exemplo, é usual considerar que, para um dada problema, não há perdas de calor, o atrito é nulo, etc. Este tipo de erros fogem ao controlo do analista numérico e são muito difíceis de quantificar.

Aritmética computacional 6

O aspecto que vamos agora analisar tem a ver com os erros que resultam da forma como representamos os números reais. O conjunto dos números reais R não pode ser representado numa máquina de precisão finita. Numa máquina só é possível representar um seu subconjunto finito F. Os números desse conjunto F são chamados números de vírgula flutuante. Um número real x é geralmente truncado pela máquina dando origem a um novo número que (número de vírgula flutuante), que se designa por fl(x). Em geral, x 6 = fl(x). Além disso, podemos ter x 1 6 = x 2 e fl(xi) = fl(x 2 ). Usualmente, um computador guarda um número real na forma

x = (−1)s^ · (0,a 1 a 2... at) · βe^ = (−1)s^ · m · βe−t, a 1 6 = 0, (1.1)

onde s é 0 ou 1, conforme o sinal de x, β, inteiro positivo maior ou igual a 2, é a base adoptada pelo computador específico em que estamos a trabalhar, m é um inteiro chamado mantissa, cujo comprimento t é o número máximo de algarismos armazenados ai, com 0 ≤ ai ≤ β − 1 , e e um número inteiro chamado expoente. Os dígitos a 1 a 2... ap, com p ≤ t, são chamados os p primeiros algarismos significativos de x. O conjunto F fica completamente caracterizado à custa de 4 parâmetros: a base β, número de algarismos significativos t, e o intervalo de variação do expoente e, designado por ]L, U [, com L < 0 e U > 0. Escrevemos então F(β, t, L, U ). Em Matlab temos F = F(2, 53 , − 1021 , 1024). Note-se que, 53 algarismos significativos em base 2 correspondem a 15 algarismos significativos em base 10. Para perceber melhor o que está em causa, consideremos, por exemplo, o número x = 123, 9346. Este número não tem representação numa máquina de base decimal cuja mantissa só permita armazenar 6 dígitos. Temos assim necessidade de o aproximar por um outro que possa ser representado na referida máquina. Essa aproximação vai ser efectuada por um processo conhecido por arredondamento. A forma de arredondar um número real é a usual. Como tal

x = 123, 9346 ≈ 123 ,935 = x,

e este novo valor já tem representação na máquina que estamos a usar sob a forma 0 , 123935 × 102. Note-se que o arredondamento foi efectuado na terceira casa decimal e que

|e(x)| = |x − x| = 0, 0004 < 0 , 5 × 10 −^3 ,

r(x) = |e(x)| |x|

≈ 3 , 23 × 10 −^6 < 5 × 10 −^6.

Se o arredondamento tivesse sido efectuado na segunda casa decimal vinha

x = 123, 9346 ≈ 123 ,93 = x,

e assim |e(x)| = 0, 0045 < 0 , 5 × 10 −^2 ,

r(x) = |e(x)| |x|

≈ 3 , 63 × 10 −^5 < 5 × 10 −^5.

Daqui resultam as seguintes definições.

Aritmética computacional 8

para x < 0 , pode dar resultados catastróficos. Outro exemplo, pode ser dado no cálculo do desvio padrão. Como se sabe, a média e o desvio padrão de uma amostra {x 1 , x 2 ,... , xn} são dados por

x ¯ =

n

∑^ n

i=

xi, e σ =

n − 1

∑^ n

i=

(xi − ¯x)^2

A fórmula equivalente para o desvio padrão

σ =

n − 1

( (^) n ∑

i=

x^2 i − nx¯

é muitas vezes usada mas o cancelamento subtractivo é mais perigoso que na fórmula anterior.

Os erros de truncatura ou de discretização são, por definição, os erros que surgem quando se passa de um processo infinito para um processo finito ou quando se substitui um processo contínuo por um discreto. A título de exemplo, considere-se o conhecido Teorema do Valor Médio, estabelecido por Joseph Louis Lagrange (1736-1813).

Teorema 1.1 (Valor Médio de Lagrange) Se f for uma função contínua em [a, b] e diferenciável em ]a, b[ então existe pelo menos um ξ ∈]a, b[ tal que

f ′(ξ) = f (b) − f (a) b − a

-0.4 -0.2 0.2 0.4 0.6 0.8 1

2

3

4

a b

Figura 1.1: Teorema do Valor Médio.

Este resultado justifica o procedimento (muito comum) de substituir o cálculo da deri- vada de uma função definida num intervalo (pequeno) [a, b] pela diferença dividida

f [a, b] =

f (b) − f (a) b − a

isto é, para um valor de h pequeno

f ′(x) ≈

f (x + h) − f (x) h

O erro cometido nesta aproximação é um erro de truncatura.

Aritmética computacional 9

Também se comete um erro de truncatura quando se efectua a aproximação

e ≈

M

)M

Outro exemplo onde surgem este tipo de erros é dado pela chamada aproximação de Taylor que iremos considerar na próxima secção.

1.3 O polinómio de Taylor

Seja f uma função real definida num intervalo [a, b] ⊆ R. Um problema que frequentemente se coloca é o de determinar uma função g definida em [a, b] tal que |f (x) − g(x)| < ǫ, para todo o x ∈ [a, b], com ǫ > 0 uma tolerância dada. A existência de solução para tal problema é dada pelo Teorema de Weierstrass, devido a Karl Wilhelm Theodor Weierstrass (1815-1897).

Teorema 1.2 (Weierstrass) Seja f uma função contínua definida em [a, b]. Então para cada ε > 0 existe um polinómio p definido em [a, b] tal que

max x∈[a,b]

|f (x) − p(x)| < ε.

Notemos a grande importância deste resultado. De acordo com ele, podemos ter a certeza que dada uma função contínua f qualquer existe sempre um polinómio p que está tão próximo de f quanto se queira. Assim sendo, este resultado legitima a aproximação polinomial, isto é, a tarefa de, dada uma função, procurar um polinómio que a aproxime. No entanto, o teorema anterior não nos diz como podemos construir esse polinómio; ele apenas garante a existência. Consideremos agora o seguinte teorema, apresentado sem demonstração, devido a Brook Taylor (1685-1731).^1

Teorema 1.3 (Taylor) Se f admite derivadas contínuas até à ordem n (inclusivé) em [a, b], isto é, se f ∈ Cn([a, b]), e se f (n+1)^ existir em ]a, b[ então, para todo o x, x 0 ∈ [a, b],

f (x) = Pn(x; x 0 ) + Rn(x; x 0 ), (1.2)

onde

Pn(x; x 0 ) =

∑^ n

k=

f (k)(x 0 ) k!

(x − x 0 )k

e

Rn(x; x 0 ) = f (n+1)(ξ) (n + 1)!

(x − x 0 )n+1, ξ ∈ I{x, x 0 },

sendo I{x, x 0 } o intervalo aberto definido por x e x 0.

A (1.2) chamaremos fórmula de Taylor sendo Pn(x; x 0 ) o polinómio de Taylor de f em torno do ponto x 0 e Rn(x; x 0 ) o resto (de Lagrange) de ordem n (ou de grau n + 1). Se x 0 = 0 a (1.2) chamaremos fórmula de Maclaurin.^2

(^1) Taylor foi, entre outras coisas, o sucessor de Edmond Halley (1656-1742) como secretário da Royal Society. Publicou, em 1715, um livro intitulado Methodus Incrementorum Directa & Inversa no qual a sua expansão aparece descrita. O seu teorema foi enunciado em 1712. (^2) Colin Maclaurin (1698-1746) foi um menino prodígio sendo nomeado professor em Aberdeen com a idade de 19 anos. A sua expansão apareceu em 1742 no Treatise on Fluxions.

Aritmética computacional 11

Por tentativas... n = 9 ⇒

210 = 0, 254 × 10 −^2

n = 10 ⇒

211 = 0, 462 × 10 −^3.

Logo a aproximação pedida é

e^2 ≈

∑^10

k=

xk k!

1.4 Problemas

Exercício 1.3 Sejam x, y e z três quantidades exactas. Por arredondamento obtiveram-se as seguintes aproximações: x = 231, y = 2, 31 e z = 23, 147.

1. Conte o número de casas decimais correctas nas aproximações e calcule limites superiores para o erro absoluto em cada uma delas. Compare os resultados e comente.
2. Conte o número de algarismos significativos correctos nas aproximações e calcule limites superiores para o erro relativo em cada uma delas. Compare os resultados e comente.

Exercício 1.4 (Matlab) A adição em vírgula flutuante não verifica a propriedade associativa. Para comprovar esse facto, calcule, usando o computador, a soma x+y+z nas formas x+(y+z) e (x + y) + z para os casos em que:

1. x = 1, 0 ; y = − 5 , 0 ; z = 6, 0 ;
2. x = 1 × 1020 ; y = − 1 × 1020 ; z = 1, 0.

Explique os resultados obtidos.

Exercício 1.5 (Matlab) Escreva um programa para obter os primeiros n termos da sucessão gerada pela equação de diferenças

x(k+1)^ = 2, 25 x(k)^ − 0 , 5 x(k−1),

a partir dos valores iniciais x(1)^ = 13 e x(2)^ = 121.

1. Considere n = 60 e trace o gráfico semi-logarítmico dos resultados obtidos como função de k.
2. A solução exacta da equação de diferenças é dada por (prove)

x(k)^ =

41 −k 3

que decresce, de forma monótona, quando k cresce. Explique porque é que o gráfico obtido na alínea anterior não está de acordo com este comportamento teórico.

Exercício 1.6 (Matlab) Considere o desenvolvimento em série de Maclaurin

ex^ =

∑^ ∞

k=

xk k!

Calcule e−^12 usando este desenvolvimento. Repita novamente o cálculo mas usando o facto de e−^12 = (^) e^112. Compare os resultados e explique o sucedido.

Aritmética computacional 12

Exercício 1.7 O fluxo através de uma parte da camada fronteira num fluído viscoso é dado pelo integral definido (^) ∫ 0 , 8 0

1 ,4(1 − e−^4 x 2 )dx.

Usando a fórmula de Taylor na função integranda, aproxime o valor do integral com quatro casas decimais correctas.

Exercício 1.8 Consideremos uma viga uniforme de comprimento L, suspensa, sujeita a uma carga uniformemente distribuída, W , e a uma força compressiva, P , em cada extremo. A deflexão, D, no ponto médio é dada por

D =

W EI

P 2

(sec (0, 5 mL) − 1) −

W L^2

8 P

onde m^2 = P/EI, com E e I constantes. Usando o desenvolvimento em série de Maclaurin da função y = sec x, prove que, quando a força gravítica tende a anular-se, a deflexão, D, tende

para

5 W L^4

384 EI

Exercício 1.9 A lei dos gases perfeitos é dada por P V = nrT e relaciona a pressão, P , o volume, V , a temperatura, T , e o número de moles, n, de um gás ideal. O número r nesta equação depende apenas do sistema de medição a usar. Suponhamos que foram efectuadas as seguintes experiências para testar a veracidade da lei usando o mesmo gás.

1. Consideraram-se P = 1, 0 atmosferas, n = 0, 0042 moles, V = 0, 10 metros cúbicos e r = 0, 082. Usando a lei, a temperatura do gás foi prevista como sendo

T =

P V

nr

1 , 0 × 0 , 10

0 , 082 × 0 , 0042

= 290o^ Kelvin = 17o^ Celsius.

Quando medimos a temperatura do gás verificámos ser 17 o^ Celsius.

2. A experiência anterior foi repetida usando os mesmos valores de r e n mas aumentando o pressão quatro vezes enquanto se reduziu o volume na mesma proporção. Como P V é constante, a temperatura prevista é de 17 o^ Celsius mas agora, ao medir a temperatura do gás, encontrámos o valor 32 o^ Celsius.

Será que a lei não é válida nesta situação?

Sistemas de equações lineares 14

2.1 Classes de matrizes

Existem vários tipos de matrizes com relevância em aplicações práticas.

• Matrizes densas e matrizes esparsas Um matriz com muitos elementos nulos diz-se esparsa; caso contrário diz-se que a matriz é densa. Sistemas com matrizes esparsas modelam sobretudo problemas onde existem princípios de influência local. Note-se que, no caso da ponte, as equações em cada nodo apenas envolvem as barras que aí se encontram. O seu número é o mesmo quer a ponte tenha 50 metros e, digamos, 10 barras, ou 5 km e 1000 barras. Assim, para uma ponte grande, a maioria dos coeficientes da matriz são nulos. Como caso particular das matrizes esparsas temos as matrizes banda – uma matriz A = (aij ) é uma matriz banda se aij = 0 para todo |i−j| > β, onde β é o comprimento de banda de A – e, dentro dessa classe, as chamadas matrizes tridiagonais, isto é, as matrizes com comprimento de banda β = 3. Por outras palavras, as matrizes banda são aquelas cujos elementos não nulos se concentram apenas num conjunto de diago- nais paralelas à diagonal principal.
• Matrizes triangulares As matrizes triangulares são aquelas que possuem todos os seus elementos acima ou abaixo da diagonal principal iguais a zero. No primeiro caso, as matrizes dizem-se triangulares inferiores e no segundo triangulares superiores.
• Matrizes simétricas As matrizes simétricas são aquelas que coincidem com a sua transposta. Uma caracte- rística importante das matrizes simétricas é o facto de todos os seus valores próprios serem reais.
• Matrizes estritamente diagonal dominantes (EDD) Uma matriz A = (aij )ni,j=1 diz-se estritamente diagonal dominante por linhas se

|aii| >

∑^ n

j=1,j 6 =i

|aij |, i = 1, ..., n,

e estritamente diagonal dominante por colunas se

|ajj | >

∑^ n

i=1,i 6 =j

|aij | j = 1, ..., n.

A matriz diz-se estritamente diagonal dominante se for estritamente diagonal domi- nante por linhas ou estritamente diagonal dominante por colunas.

• Matrizes simétricas e positivas definidas (SPD) Uma matriz A, diz-se simétrica e positiva definida se for simétrica e se, para todo o vector x ∈ Rn, não nulo, se tem xT^ Ax > 0. Nesse caso diz-se que A ∈ SPD.

Sistemas de equações lineares 15

Exercício 2.1 Mostre que A ∈ SPD se e só se todos os valores próprios de A são reais e positivos.

Resolução: Por definição, λi é valor próprio de A ∈ Mn(R) se e só se

Axi = λixi, xi ∈ Rn, xi 6 = 0,

sendo xi o vector próprio de A associado a λi. Suponhamos que A ∈ SPD. Então xT^ Ax > 0 é válida para todo o vector não nulo x e, em particular, para os vectores próprios xi de A. Assim sendo,

0 < xTi Axi = xiλixi = λi‖x‖^22 ⇒ λi > 0.

Suponhamos agora que todos os valores próprios de A são positivos. Então, pelo que foi visto, para todos os valores próprios xi de A é válida a desigualdade xTi Axi > 0. Como A é uma matriz real simétrica, prova-se que existe uma base ortonormada de vectores próprios de A (teorema espectral). Nessa base, qualquer que seja o vector não nulo x, tem-se que existem constantes ci, i = 1, 2 , ..., n, onde n é a ordem da matriz A, tais que x = c 1 x 1 + c 2 x 2 + · · · + cnxn. Então Ax = c 1 Ax 1 + c 2 Ax 2 + · · · + cnAxn o que implica Ax = c 1 λ 1 + c 2 λ 2 + · · · + cnλn. Assim sendo, xT^ Ax = c^21 λ 1 + x^22 λ 2 + · · · + c^2 nλn > 0 , o que prova o pretendido.

2.2 Métodos directos: revisão

Consideremos, de novo, o problema de determinar o vector x ∈ Rn^ tal que Ax = b, sendo b ∈ Rn^ e A ∈ Mn(R). O primeiro tipo de métodos que iremos considerar resolver este problema são os chamados métodos directos. Estes métodos são aqueles que, supondo não haver erros de arredondamento ou quaisquer outros, nos permitem obter a solução exacta do problema num número finito de operações aritméticas. Os métodos directos baseiam-se no processo de eliminação de Gauss que consiste em transformar o sistema Ax = b num sistema equivalente U x = c, onde U é uma matriz triangular superior, atravéz de operações elementares efectuadas na matriz ampliada (lembrar estes conceitos dados na disciplina de álgebra linear). O sistema a resolver pode ser escrito na forma (^)   

 

u 11 x 1 + u 12 x 2 + · · · + u 1 nxn = c 1 u 22 x 2 + · · · + u 2 nxn = c 2

... .. .

unnxn = cn

e a sua resolução, caso uii 6 = 0, i = 1,... , n, é feita de acordo com o algoritmo seguinte.

Algoritmo 2.1 Resolução de um sistema triangular superior

Dados: ci, i = 1,... , n, e uij , i = 1,... , n, j = i,... , n xn := cn/unn Para i de n − 1 até 1 fazer xi :=

ci −

∑n j=i+1 uij^ xj

/uii Resultado: xi, i = 1,... , n

Sistemas de equações lineares 17

Algoritmo 2.3 Resolução de um sistema tridiagonal

Dados: matriz A e vector b Determinar as matrizes L e U (∗ Resolver Ly = b ∗) y 1 := b 1 /l 1 Para i de 2 até n fazer yi := (bi − αiyi− 1 )/li (∗ Resolver U x = y ∗) xn := yn Para i de n − 1 até 1 fazer xi := yi − uixi+ Resultado: xi, i = 1,... , n

Prova-se o seguinte resultado.

Teorema 2.1 Para uma matriz A ∈ Mn(R), a factorização A = LU existe e é única se e só se as submatrizes principais de A forem não singulares.

Quando, no processo da factorização de Gauss, encontramos um pivot nulo, podemos trocar de linhas por forma a evitar a divisão por zero. Uma possibilidade consiste em escolher, em cada passo da iteração, o pivot de módulo máximo. Obtém-se, assim, a factorização P A = LU,

com P a matriz de permutação da linha efectuada.

Note-se que, se A ∈ EDD ou A ∈ SPD então a matriz A verifica as hipóteses do teorema anterior e, como tal, a factorização LU de A existe e é única. Mais ainda, se A ∈ SPD pode obter-se a factorização de de A na forma A = RRT^ , com R uma matriz triangular inferior com elementos diagonais positivos. Esta factorização é chamada factorização de Cholesky, em homenagem a André-Louis Cholesky (1875-1918). No caso 2 × 2 , por exemplo, [ a 11 a 12 a 12 a 22

]

[

r 11 0 r 12 r 22

] [

r 11 r 12 0 r 22

]

implica que r 11 =

a 11 , r 12 = a 12 /r 11 , r 22 =

a 22 − r^212.

2.3 Normas de matrizes. Condicionamento

Seja Mn(R) o conjunto das matrizes quadradas de ordem n com coeficientes reais. Como este conjunto é um espaço vectorial podemos nele definir uma norma.

Definição 2.1 (Norma matricial) Seja ‖.‖ uma norma vectorial definida em Rn. A aplicação ‖ · ‖ : Mn(R) −→ R+ 0 tal que, para todo o A ∈ Mn(R),

‖A‖ = sup x 6 =

‖Ax‖ ‖x‖

é designada norma matricial subordinada à norma vectorial ‖.‖.

Sistemas de equações lineares 18

No entanto, a generalização para o caso rectangular Atendendo à definição anterior, podemos dizer que a norma de uma matriz mede o alongamento máximo que essa matriz provoca em qualquer vector não nulo em relação à correspondente norma vectorial. Prova-se que a função estabelecida na definição anterior é uma norma, isto é, que verifica as seguintes propriedades:

1. ∀A ∈ Mn(R), ‖A‖ = 0 ⇔ A = 0 (matriz nula),
2. ∀A ∈ Mn(R), ∀λ ∈ R, ‖λA‖ = |λ|‖A‖,
3. ∀A, B ∈ Mn(R), ‖A + B‖ ≤ ‖A‖ + ‖B‖,

Prova-se também que as seguintes aplicações verificam as propriedades que caracterizam um norma, com A ∈ Mn(R) a matriz de elemento genérico aij :

• ‖A‖ 1 = max j=1,...,n

∑n

i=

|aij |;

• ‖A‖ 2 =

ρ(AT^ A), com ρ(A) o raio espectral de A, definido por

ρ(A) = max i=1,...,n {|λi| : λi é valor próprio de A};

• ‖A‖∞ = max i=1,...,n

∑n

j=

|aij |.

• ‖A‖F =

√√∑^ n i,j=

a^2 ij.

A norma ‖ · ‖F é chamada norma de Frobenius, em homenagem a Ferdinand Georg Frobenius (1849–1917), e não é soburdinada a nenhuma norma vectorial.

Exercício 2.3 Considere a matriz

A =

Calcule ‖A‖ 1 , ‖A‖∞ e ‖A‖F.

Resolução: Vê-se facilmente que

‖A‖ 1 = max{| 2 | + | − 1 | + | − 2 |, | − 1 | + | − 2 | + | − 1 |, | 0 | + | 1 | + | 0 |} = 5,

‖A‖∞ = max{| 2 | + | − 1 | + | 0 |, | − 1 | + | − 2 | + | 1 |, | − 2 | + | − 1 | + | 0 |} = 4 e

‖A‖F =

22 + (−1)^2 + 0^2 + (−1)^2 + (−2)^2 + 1^2 + (−2)^2 + (−1)^2 + 0^2 = 4.