64
______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 1 CALCULO NUMÉRICO NOTAS DE AULA Este trabalho não tem objetivo de publicação e foi desenvolvido apenas na tentativa de aproveitar melhor o tempo disponível para as aulas de cálculo numérico, devido a extensão dos programas da disciplina, e não deve ser usado em hipótese alguma como única fonte de pesquisa, pelo fato de não contem maiores detalhes referentes aos assuntos aqui abordados, que serão melhor trabalhados em sala, inclusive com a resolução dos exemplos. Para um melhor aprofundamento dos conteúdos são sugeridas as referencias bibliográficas usadas na elaboração desse. Possivelmente não contenha todo o conteúdo que eventualmente faça parte do programa, que serão inseridos dentro da sequencia lógica conforme surgir a necessidade. Prof. Sergio Marcussi Gaspechak

Apostila Cálculo Numérico e Exercícios

Embed Size (px)

DESCRIPTION

Apostila de Cálculo Numérico- sistemas numéricos- método iterativo linear- método de newton raphson - método das secantes- sistemas lineares- pivotamento- sistema de matrizes- exercícios

Citation preview

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 1

CALCULO NUMÉRICO

NOTAS DE AULA

Este trabalho não tem objetivo de publicação e foi desenvolvido apenas na tentativa

de aproveitar melhor o tempo disponível para as aulas de cálculo numérico, devido a

extensão dos programas da disciplina, e não deve ser usado em hipótese alguma como

única fonte de pesquisa, pelo fato de não contem maiores detalhes referentes aos assuntos

aqui abordados, que serão melhor trabalhados em sala, inclusive com a resolução dos

exemplos. Para um melhor aprofundamento dos conteúdos são sugeridas as referencias

bibliográficas usadas na elaboração desse.

Possivelmente não contenha todo o conteúdo que eventualmente faça parte do

programa, que serão inseridos dentro da sequencia lógica conforme surgir a necessidade.

Prof. Sergio Marcussi Gaspechak

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 2

CALCULO NUMÉRICO

CAPITULO 1

Conceitos básicos do calculo numérico

1 Introdução

As fases de resolução de um problema físico podem de modo geral ser assim

representados.

A partir do problema físico, com emprego de leis de conservação, de relações

constitutivas, modelos de turbulência, das condições de contorno, etc., chega-se a um

modelo matemático.

A área da matemática que trata da concepção de processos numéricos e sua

exeqüibilidade para encontrar solução do modelo matemático denominam-se análise

numérica.

Com o surgimento do computador na década de 40, notou-se a importância da

analise numérica, uma vez que, por meio de processos eletrônicos tornaram viáveis as

técnicas numéricas.

Num curso de graduação a disciplina de calculo numérico tem por objetivo propiciar

ao estudante o conhecimento de processos já concebidos pela análise numérica.

1.1 Problema Numérico

É todo problema que é resolvido por meio de calculo numérico, onde os dados de

entrada e saída, input e output, são dados por um conjunto finito de números.

Problema físico

Modelagem

Modelo Matemático

Resolução

Implementação Computacional

Solução

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 3

1.2 Métodos numéricos

É o conjunto de procedimentos utilizados para transformar o modelo matemático a

ser resolvido em um problema numérico.

Para resolver um problema numérico deve-se escolher o método mais eficiente

levando em consideração os seguintes aspectos:

i) Precisão desejada.

ii) Capacidade do método de conduzir a resultado esperado como, convergência

e ordem de convergência.

iii) Esforço computacional despendido ( tempo de processamento, economia de

memória necessária para resolução).

1.3 Algoritmo

É a descrição seqüencial, completa, dos passos que levam a solução do problema, e

caracterizam um método numérico. O algoritmo fornece descrição completa de operações

bem definidas por das quais o conjunto de dados de entrada é transformado em dados de

saída.

1.4 Iteração ou Aproximação Sucessiva

Iteração tem um sentido de repetição de um processo. É uma das idéias básicas do

calculo numérico, isto é, a maioria dos métodos numéricos são iterativos, mas temos

também outros processos não iterativos, com por exemplo, os métodos diretos e métodos

de passos múltiplos.

Dentro do método iterativo o objetivo é criar uma sequencia de recorrência

convergente, onde o valor da convergência é a solução do problema em questão, para isso

um fator importante, é temos o erro tendendo a zero.

Um método iterativo possui as seguintes características:

i) Tentativa inicial: primeira tentativa para obter a solução desejada.

ii) Equação de recorrência: equação por meio da qual será feita as iterações

com o objetivo de obter a solução desejada.

iii) Teste de parada: condição para que o processo iterativo seja finalizado.

1.5 Aproximação local

A idéia básica da aproximação local é substituição de uma função por outra, mais

fácil de trabalhar. Como por exemplo, subustituir uma função não linear por outra linear em

um determinado intervalo.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 4

1.6 Erros

Os erros merecem atenção especial, pois ao tentar encontrar a solução de um

problema através de calculo numérico, surgem varias fontes de erros, que podem levar a

resultados muitos distantes da solução ou que nada tenham a ver com a solução do

problema.

As principais fontes de erros são:

i) Erros nos dados de entrada.

ii) Erros na fase de modelagem.

iii) Erros de arredondamento.

iv) Erros de truncamento.

v) Erros humanos e de máquinas.

1.6.1 Tipos de Erros

Quando calculamos a solução de um problema por aproximação temos de ter

critérios de como delimitar o erro envolvido, neste caso usaremos os tipos de erros assim

definidos:

EA: erro Absoluto

aE = n

x x

ER: erro Relativo

rE = n

x x

x

O erro relativo é o erro percentual do resultado, pois um erro significativo para um

resultado, pode não significar nada em outro resultado se analizada sua proporcionalidade,

ou seja, leva em consideração a ordem de grandeza do valor calculado.

Temos que x é a solução exata do problema e n

x é a aproximação para a solução

na n-ésima iteração.

O problema em encontrar o erro é que não conhecemos x e apenas aproximações

para este, fornecida pelo processo iterativo.

Neste caso podemos estimar o erro da seguinte forma:

aE < 1n n

x x

no caso do erro absoluto

rE < 1

1

n n

n

x x

x

no caso do erro relativo

Onde nx e 1nx é a solução aproximada obtida nas duas ultimas iterações.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 5

Quando 0,5 10 t

aE então dizemos que a aproximação tem possui t casas

decimais corretas, ou t casas decimais de precisão.

OBSERVAÇÃO: TODO PROBLEMA RESOLVIDO ATRAVÉS DE MÉTODOS NUMÉRICOS

DEVE CONTER UMA RESPOSTA COM A PRECISÃO DO RESULTADO, EXCETO PARA OS

MÉTODOS EXATOS, SEM ISSO A SOLUÇÃO ENCONTRADA NÃO FARÁ NENHUM SENTIDO.

Pois se você encontrou o valor de 3.1415 esta errado!, mas se for dito que o

valor de 3.1415 com 3 casas decimais corretas, ai sim faz sentido sua resposta, pois

fica claro que temos um erro apartir da 4 casa decimal.

Ex 1) Vamos supor que durante o processo iterativo pra obter a raiz de uma equação

obteve-se os seguintes resultados:

10 2.34567x e 11 2.34538x assim o erro absoluto será

3

10 11 2.34567 2.34538 0,00029 0,5.10EA x x

note que o erro é inferior a 30,5.10

, assim 3t e a resposta ficara:

A raiz da função é 2.34538 com 3 casa decimais corretas.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 6

CAPITULO 2

SISTEMAS NUMERICOS

Todo numero N, pode ser representado em diferentes bases. O sistema numérico

que utilizamos é o sistema decimal, posicional, visto que todo número é escrito como soma

de produto de potencias de base 10, e cada algarismo, dígito, vale de acordo com a sua

posição.

Outras bases, e outros sistemas numéricos, já foram usados na antiguidade como a

base 60, base 12 e o sistema aditivo conhecido por algarismos romanos pro exemplo.

Vejamos Alguns exemplos:

Ex 1): 10(4362) = 3 2 1 04.10 3.10 6.10 2.10

2.1 Sistema Binário

Atualmente com o uso do computador uma base muito usada é a base 2, também

conhecido como sistema binário, pois na base dois todos os números podem ser escritos

como soma de potencias de base 2, onde temos apenas dois dígitos 0 e 1. 3 2 1 0

10 2(13) (1101) 1 2 1 2 0 2 1 2

2.2 0utras bases

Base 5

Ex 2: 3 2 1 0

10 5(176) 1 5 2 5 0 5 1 5 (1201)

Qualquer número natural > 0 pode ser usado como base, quanto estamos na base

dez omitimos o índice que indica a base usada.

2.3 Sistemas discretos de números no computador

Alguns números não possuem representação exata em determinadas bases, por esse

fato, alguns sistemas são ditos discretos, ou seja possuem falhas, ou não permitem a

representação exata de alguns valores e alguns casos nem mesmo uma representação

aproximada.

Inicialmente, vejamos como os números são representados num computador.

A representação de um número inteiro N não possui maiores dificuldades é dada por:

0 1 1

1 1 0 0 1 1

0

( ... ) ( ... )k k

k k k k

k

i

i

i

N d d d d d d d d

d

Onde i

d , 0, 1, 2,...i k são inteiros satisfazendo 0i

d e 0k

d .

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 7

3 2 1 03) ) 2356 2.10 3.10 5.10 6.10

2000 300 50 6

Ex i

3 2 1 0

2

10

) (1011) 1.2 0.2 1.2 1.2

8 0 2 1 (11)

ii

2.4 Representação de um número real.

Podemos representar num computador um número real de duas maneiras;

i) Ponto fixo: Usado durante muito tempo.

Se x é um número real não nulo então sua representação será dada por: n

i

ii k

x d

Onde k e n são inteiros satisfazendo k < n e usualmente k > 0 e n > 0 e os valores de

id satisfazem 0i

d .

Ex 4:

1352,48 2

3

i

i

i

d

3 2 1 0 1 2

3 2 1 0 1 2( )d d d d d d

3 2 1 0 1 21352,48 (1.10 3.10 5.10 2.10 4.10 8.10 )

ii) Ponto flutuante: Usado nos dias atuais.

Se x é um número real não nulo este será representado em ponto flutuante por:

. ex d onde é a base, e é o expoente, d é a mantissa, parte decimal de um

número, escrita como ponto fixo.

1

1( ... )

ni k k n

i k k ni k

d d d d d

Obs: 1

0d caracteriza um sistema de ponto flutuante normalizado.

O zero pode ser representado em qualquer sistema, com mantissa igual a zero e

expoente –m.

Ex 5: Se a base é 10 então:

2 1 0 1 2 3)235,231 2 10 3 10 5 10 2 10 2 10 1 10i

3= 0,235231 10

esta na forma normalizada.

ii) 312,63 0,01263 10 não esta na forma normalizada.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 8

2.5 Dígitos significativos

Definição: Seja a base do sistema de números em ponto flutuante. Os dígitos

significativos de um número x são todos os algarismos de 0 a 1 , desde que x esteja

escrito na forma de ponto flutuante normalizado.

Ex 6:

253,26 0,5326 10 possui 4 dígitos significativos.

30,000230 0,230 10 possui 3 dígitos significativos.

Para representarmos um sistema de números em ponto flutuante normalizado, na

base , com t dígitos significativos e com limites do expoente -m e M, usaremos a

notação. ( , , , )F t m M .

Ex 7: Represente os números abaixo no sistema F(10,3,2,2).

i) 0,25 = ...

ii) -5,149 = ...

iii) 0,0147 = ...

iv) 5478,23 = ... overflow

v) 0,0003 = ... underflow

Overflow: Termo corrente na computação para designar que um valor não pode ser representado pois seu valor excede a capacidade de armazenamento disponível. Underflow: Termo corrente na computação para designar que um valor não pode ser representado pois está contido entre 0 e o menor valor real normalizado representável.

Observe que os números dos itens iv) e v) não possuem representação no sistema

F(10,3,2,2). Como foi visto no exemplo anterior nem todos os números possuem

representação exata em um sistema ( , , , )F t m M , por isso é dito sistema discreto de

números.

Quando não podemos representar exatamente um número em um sistema

( , , , )F t m M , podemos tentar um arredondamento para uma melhor representação

dentro do sistema, isto é, uma aproximação com um erro mínimo possível, o que nem

sempre é possível.

Dentro de um sistema ( , , , )F t m M a representação ficara da forma:

exp

... ...

Sinal Sinal doValor damantissa Valor dodamantissa Expoente oente

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 9

Ex 7: A representação do número 10( 23) no sistema (2,10,10,10)F sera:

1 1 0 1 1 1 0 ... 0 0 0 0 0 0 1 0 1

Precisão Simples = 32 bits= 1 / 8 / 23 (Precisão Simples = 32 bits= 1 / 8 / 2Sinal / Expoente / Mantissa)

Precisã Dependendo da maquina e da linguagem utilizada é possível se trabalhar com

precisão simples, 32 bits, ou precisão dupla, 64 bits, para os dígitos da mantissa.

É importante observar que a precisão dupla aumenta significativamente o tempo de

execusão e a memória, logo deve ser usado somente em casos onde houver necessidade.

2.6 Regras de Arredondamento em Ponto Flutuante.

Definição 1) Arredondar um número x dentro de um sistema de numeração, por outro

número, com menor número de dígitos significativos, consiste em encontrar outro número

x que pertence ao sistema de numeração tal x x seja menor possível.

Dado um número x , seja x sua representação em ( , , , )F t m M através do

arredondamento, se x = 0 então x = x , caso contrario escolhemos S e e tais que:

ex S onde 1

1 1(1 ); (1 ) .

2 2

t tS

Se [ , ]e m M não temos condições de representar o número no sistema. Se

[ , ]e m M então calculamos:

1 2 3 10. ... ...

2

t

t ts d d d d d

e truncamos em t dígitos. Assim

1 2 3(0. ... ) e

tx x d d d d

Ex 8: Considere o sistema F(10,3,5,5), represente neste sistema os números

1 2 3 4 51234,56; 0,0005492; 0.9995; 123456,7; 0,0000001 x x x x x

Resolução:...

----------------------------------///----------------------------------

2.7 Operações Aritméticas com ponto Flutuante.

Nem todas as propriedades dos números reais são válidas dentro de um sistema de

ponto flutuante como, por exemplo, as propriedades comutativas para a adição e para a

multiplicação e divisão assim como as propriedades associativa da multiplicação em relação

a adição.

Ex 9: Considere o sistema de base 10 com 3 dígitos significativos. Efetue as operações

indicadas:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 10

i) (11,4 3,18) 5,05 e 11,4 (3,18 5,05)

ii) 3,18 11,4

5,05

e

3,1811,4

5,05

iii) 3,18 (5,05 11,4) e 3,18 5,05 3,18 11,4

2.8 Mal condicionamento

A maioria dos processos numéricos segue a seguinte linha geral:

Dados são fornecidos

Os dados são processados de acordo com um plano pré-estabelecido (algoritmo)

Resultados são produzidos

Ao tentar resolver um determinado problema podemos nos deparar com o que

chamamos de Mau condicionamento, ou problema mal condicionado. Neste tipo de

problemas temos que uma pequena perturbação nos dados de entrada fornece uma grande

variação nos resultados obtidos.

Ex 10: Considere o sistema

2

1,01 2,01

x y

x y

Graficamente temos:

Note que as retas são aparentemente coincidentes, o que não é verdadeiro.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 11

CAPITULO 3

ZEROS DE FUNÇÕES

3.1 Equações não lineares.

Um dos problemas que ocorre mais frequentemente em trabalhos científicos é

calcular as raízes de equações da forma:

( ) 0f x

Definição 1:

Se :[ , ]f a b R é uma função dada, um ponto [ , ]x a b é um zero de

( )f x se for uma raiz da equação ( ) 0f x .

Definição 2: Multiplicidade

Uma raiz x tem multiplicidade q se conseguirmos escrever

( ) ( ) ( )qf x x x g x com ( ) 0g x .Por conseqüência temos que se a raiz tem

multiplicidade q então: ( )( ) 0, '( ) 0, "( ) 0, '''( ) 0, ... , ( ) 0qf x f x f x f x f x

3.2 Teorema de Bolzano

Se uma função continua ( )f x assume valores de sinais opostos, isto é, muda de

sinal nos extremos do intervalo [ , ]a b teremos ( ). ( ) 0f a f b , então existe pelo

menos um valor x em [ , ]a b , tal que ( ) 0f x .

Obs: Se ( ). ( ) 0f a f b não podemos dizer nada sobre a existência ou não de raizes

em [ , ]a b .

Ex 1: Mostre que a função 1

( ) xf x ex

possui uma raiz no intervalo [0,5; 1].

Solução: ...

-------------------------------///----------------------------------

Usando o maple.

> f:=x->exp(x)-1/x;

> f(0.5);

> f(1.0);

---------------------------------///--------------------------------

Outra forma de obtermos intervalos que contenham raízes é rearranjando a equação

de ( ) 0f x conforme faremos no exemplo a seguir:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 12

Ex 2: Encontre os intervalos que contenha a(s) raízes das funções

i) 22 ( 2)( ) ( 1) 1xf x x e

Solução: ( ) 0f x assim 22 ( 2)( 1) 1 0xx e podemos rearranjar a equação

dada

por 22 ( 2 )( 1) xx e fazendo

2

1( 1)y x e

2( 2 )

2

xy e colocando as duas curvas

no mesmo gráfico, obtemos:

Note que temos uma raiz no intervalo (-2, -1) e outra no intervalo (0.5, 1.5).

Usando o maple:

> plot([(x+1)^2,exp(2-x^2)],x=-3..3);

Assim f possui raízes nos intervalos (-2,-1) e (1.5; 2.5).

ii) 3( ) 3 1f x x x

iii) 2( ) ( )f x x sen x

Solução:...

-----------------------------------------///-----------------------------------------------

O teorema de Bozano fornece o primeiro método para encontrar raízes de funções.

3.3 Método do Meio Intervalo (MMI), ou Bisseção.

O MMI consiste em fazer um refinamento no intervalo que contem a raiz, dividindo

sucessivamente o intervalo [a,b] que contem a raiz, mas descartando sempre os sub-

intervalos onde f não muda de sinal. O que pode ser executado através do algoritmo.

1 1

1 1

( , ) ( ) ( ) 0( , )

( , ) ( ) ( ) 0

k k k k

k k

k k k k

a m se f a f ma b

m b se f a f m

com

1 1

2

k kk

a bm

Onde k são ditos passos para encontrar a raiz com k = 1,2,3,...

Ex 3: Encontre a raiz da função 1

( ) xf x ex

através do MMI.

Solução:...

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 13

---------------------------------///------------------------------

3.4 Critério de parada

Como podemos ver no ex. anterior o processo pode ser realizado indefinidamente,

logo precisamos estabelecer um critério de parada.

O critério pode definido por um número finito de passos e/ou uma precisão

desejada, ou seja, satisfatória para o seu resultado, o que pode ser verificado calculando o

erro absoluto do resultado.

3.5 Método Iterativo Linear (MIL) ou Ponto Fixo.

Neste método construímos funções ( )x x a partir da equação ( ) 0f x o

que pode ser feita de infinitas formas, e realizamos iterações com a regra de recorrência

dada por 1

( )n n

x x , a sequência gerada pode convergir para a raiz, da função.

Ex 4: Construa regras de recorrências para tentar obter a raiz da função 2( ) 2f x x x e use duas delas para tentar encontrar a raiz com 3 casa decimais

corretas.

Solução ...

i) 2

12

n nx x

ii)

12

n nx x

Graficamente, esse dois casos se comportam da seguinte forma:

Como podemos ver seria um tanto quanto trabalhoso encontrar a raiz por tentativa,

mas para evitar isso podemos usar um critério na escolha da relação ( )x x .

3.6 Critério de convergência para MIL

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 14

Seja 0 1 . suponhamos que a sequência (n

x ) seja tal que

2 1 1n n n nx x x x

para todo n R . Afirmamos que (

nx ) é de Cauchy e,

portanto convergente.

Assim para a sequência gerada pelo MIL teremos:

2 1 1n n n n

x x x x

1 1

1

1

1

( ) ( )

( ) ( )1 '( ) 1 quando 0

n n n n

n n

n n n

n n

x x x x

x xx x x

x x

Com isso o critério para escolha da relação ( )x x pode ser feita de modo que

tenhamos '( ) 1n

x no intervalo que contenha a raiz da função, isto é, devemos

encontrar um intervalo que cotenha a raiz e depois verificar se '( ) 1n

x nesse

intervalo então a regra de recorrência 1

( )n n

x x ira gerar uma sequência de valores

que converve para raiz da função, caso contrario não teremos garantia de convergência e

devemos encontrar outra relação ( )x x

Para driblar a dificuldade de resolvermos a desigualdade '( ) 1n

x podemos

verificar se 0

'( ) 1x , mas para isso devemos estar certos de que 0

x x , ou seja, que

temos um bom valor inicial 0x , “chute inicial”.

Ex 5:

Dada a função ( ) 2 ln( 1) 2f x x x

i) Encontre duas regras de recorrência 1 ( )n nx x “método iterativo linear” para tentar

calcular uma raiz da função.

ii) Mostre quais das regras obtidas no item i) converge para as raízes de f .

iii) Encontre a maior raiz, através do método da iteração linear, com três casas decimais

corretas.

Um exemplo do MIL, usando o Maple:

Método Iterativo Linear

( Número de Iterações ou Erro maximo )

> restart:

> Digits:=30;

> n:=0:

> x[n]:=-0.2;

> phi:=x->(x^2-ln(x+2))/2;

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 15

> p:=diff(phi(x),x);

> g:=unapply(p,x);

> evalf(g(-0.2));

> ex:=1:

> while ex > 0.5*10^(-2) and n < 20 do

x[n+1]:=phi(x[n]):

ex:=abs(x[n]-x[n+1]):

n:=n+1:

end do;

> s:=x->(x^2-2*x-ln(x+2));

> s(-0.1);s(-0.3);

3.7 Ordem de convergência

A ordem de convergência de um método mede a velocidade com que as iterações

produzidas por esse método aproximam-se da solução exata. Assim quanto maior a ordem

de convergência do método, melhor o método.

Definição 3

Sejam k

x o resultado da aplicação de um método na iteração k e k k

e x x o

erro envolvido. Se existirem um número 1p e uma constante 0c tais que:

1lim

k

pk

k

ec

e

Então p é chamado de ordem de convergência do método.

Obs: A ordem de convergência do MIL é linear.

Como estimar o ordem de convergência a partir dos resultados.

Se supuser que k seja suficientemente grande e considerar os 3 últimos erros

calculados, então teremos:

1

1

ln

ln

k

k

k

k

ee

pe

e

Prova: ...

3.8 Método de Newton (Newton-Raphson) MNR

É um dos métodos mais populares para encontrar a raiz. O MNR é um método de

aproximações sucessivas dado por:

1

( ) 0,1,2,3,...

'( )

k

k k

k

f xx x k

f x

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 16

Existem varias formas de deduzir o método de Newton, vejamos uma interpretação

geométrica:

Vamos supor uma função f representada pelo gráfico abaixo, cujo a raiz é x .

O critério de convergência para o método de Newton é o mesmo do MIL, tomando

( )( )

'( )

f xx x

f x .

Prova: ...

Ex 6: Determinar através do MNR, com 4 casas decimais corretas a raiz da função:

( ) 4cos( ) xf x x e

Solução: ...

-----------------------///-------------------------

O MNR possui ordem de convergência quadrática, para estimar raízes simples, isto é,

raízes de multiplicidade um.

Nem sempre conhecemos a multiplicidade da raiz, os resultados sobre ordem de

convergência não são validos se a multiplicidade for maior que um. No caso do MNR se a

multiplicidade da raiz for maior que um então a ordem de convergência se torna linear.

Podemos no caso do MNR fazer a seguinte modificação, para devolver a ordem

quadrática ao método:

Se q “multiplicidade” for conhecido então:

1

( ) 0,1,2,3,...

'( )

k

k k

k

f xx x q k

f x

Se q for desconhecido, uma adaptação do método para “devolver” a ordem de convergência

do método é dada por:

( )( )

'( )

k

k

k

f xu x

f x e

1

( )

'( )

k

k k

k

u xx x

u x 0,1,2,3,...k

MNR usando o maple com q conhecido

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 17

Método de Newton-Raphson

> restart:

> Digits:=20; # Numero de digitos, 'Default': 10 digitos.

Digite a função da qual se procura a raiz

> f:=x->(exp(x)-x^6)^2;

" Função da qual se procura a raiz ";

Digite o número de iterações

> N:=20;

" Números de iterações ";

> p:=diff(f(x),x):

df:=unapply(p,x);

" Função derivada da f(x) ";

Verifique, geometricamente, a multipliciddade da raiz e o valor para o chute inicial.

> plot([f(x),df(x)],x=-1..-0.5,y=-5..5);

Digite a multiplicidade da raiz.

> m:=1;

"Multiplicidade da raiz";

Digite o valor do chute inicial

> x[0]:=0;

"chute inicial";

> "Algoritmo de NEWTON-RAPHSON"; > for n from 0 by 1 while n < N do

x[n+1]:=evalf(x[n]-m*f(x[n])/df(x[n]))

end do;

> for n from 1 by 1 while n <= N do

E[n]:=abs(x[n]-x[n-1])

end do;

> print(" A raiz é ", x[N] , " com erro máximo de ", E[N]);

3.9 Exercícios

1) Encontre através do MNR o valor de:

i) 5 21 ii) e iii)

Dica: Neste caso, basta encontrar uma função na qual o número procurado seja a raiz, ou seja, o zero da

função.

2) Mostre que a raiz da função

2

( ) ( )2

xf x sen x

tem multiplicidade 2 e encontre-

a desconsiderando a multiplicidade, considerando a multiplicidade q = 2, e considerando a

existência da multiplicidade, mas não sendo conhecida e compare os resultados:

3.10 Métodos das Secantes

Um inconveniente do MNR é ter que calcular o valor da derivada da função a cada

iteração. Podemos modificar o MNR da seguinte forma, se k for suficientemente grande

então:

1

1

( ) ( )'( ) k k

k

k k

f x f xf x

x x

pois

10

k kx x

substituindo em

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 18

1

( )

'( )

k

k k

k

f xx x

f x teremos

1 1

1

1

( ) ( )

( ) ( )

k k k k

k

k k

x f x x f xx

f x f x

conhecido por

Método das Secantes.

Ex 7: Encontre uma aproximação da raiz da função ( ) 5 xf x x e através do método

das secantes fazendo 4 iterações.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 19

CAPITULO 4

SISTEMAS LINEARES

4.1 Introdução

Dado um sistema linear de n equações e n incógnitas, do tipo.

11 1 12 2 1 1

21 1 22 2 2 2

1 1 2 2

...

...

...

n n

n n

n n nn n n

a x a x a x b

a x a x a x b

a x a x a x b

Podemos reescrevê-lo em uma notação matricial AX = B.

11 12 1 1 1

21 22 2 2 2

1 2

n

n

n n nn n n

a a a x b

a a a x b

a a a x b

Assim teremos:

11 12 1

21 22 2

1 2

n

n

n n nn

a a a

a a aA

a a a

é a matriz dos coeficientes do sistema.

1

2

n

x

xX

x

é a matriz das incógnitas, e

1

2

n

b

bB

b

é a matriz dos termos independentes.

A matriz ampliada do sistema será:

11 12 1 1

21 22 2 2

1 2

n

n

n n nn n

a a a b

a a a b

a a a b

Na matriz ampliada é como se cada linha fosse uma abreviação de cada linha do sistema.

4.2 Sistema triangulares

Um sistema é dito triangular superior se a matriz dos coeficientes do sistema for

triangular superior. Assim a sua matriz ampliada será dada por:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 20

11 12 1 1

22 2 20

0 0

0 0

n

n

nn n

a a a b

a a b

a b

4.3 Operações elementares

São ditas operações elementares toda operação que realizada nas linhas do sistema,

não altera a solução do sistema. Sendo elas:

Multiplicação da linha por um número real “escalar”

Substituição de uma linha, pela soma ou diferença dela com outra linha.

Permutação, troca da posição, de linhas.

4.4 Posto de uma matriz.

Definimos como posto de uma matriz o número de linhas não nulas, que tem pelo

menos um elemento diferente de zero, dessa matriz.

Assim definimos cP como posto da matriz dos coeficientes, aP como posto da matriz

dos coeficientes.

Assim se ao escalornarmos uma matriz pelo MEG tivermos:

a cP P então o sistema não terá solução, isto é, temos um sistema indeterminado.

a cP P n , onde n é o numero de incógnitas, então o sistema terá única solução.

a cP P n então o sistema terá infinitas soluções.

O número ( )cn P é dito nulidade do sistema, ou grau de liberdade do sistema e indica o

número de variáveis livres que teremos no sistema, isto é, o número de variáveis que

podemos atribuir qualquer valor e obter as demais soluções em função destas.

4.5 Método de Eliminação de Gauss MEG

O objetivo do MEG é transformar a matriz ampliada do sistema, em uma matriz

ampliada, cujo a matriz dos coeficientes seja triangular superior e seja equivalente a matriz

ampliada inicial. Para isso fazemos uso das operações elementares, que sistematicamente

resultara nas equações:

; 0

.

1, 2, 3, ... , (n-1)

( 1), 3, 4, ... ,

,

, 2, 3,..., (n+1)

ikik kk

kk

ij ij ik kj

am a

a

a a m a

Para cada k

i k n

e para cada i temos

j k

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 21

( 1)

1 ; ,( 1),...,1

n

i n ik k

k ii

ii

a a x

x i n na

Ex 1: Resolva o sistema abaixo, através do MEG.

3 2 11

2 2 9

2 5

x y z

x y z

x y z

Solução:...

---------------------------------------///------------------------------------------

Exercício

Resolva o sistema abaixo através do MEG.

1 2 3 4

1 2 3 4

1 2 3 4

1 2 3 4

10 5 2

2 10 2 26

2 10 2 20

3 2 10 25

x x x x

x x x x

x x x x

x x x x

4.6 Pivotamento

Um problema que pode surgir durante o processo do MEG é termos o pivô do

sistema igual a zero, ou muito próximo de zero, em algum passo k. O pivotamento consiste

em permutar linhas, no caso de pivotamento parcial e permutar linhas e colunas, no caso do

pivotamento total, com o objetivo de mudar o pivô do sistema.

i) Pivotamento Parcial

Permuta-se a linha, com outra linha de modo que o termo equivalente, mesma

coluna, possui o maior valor absoluto..

ii) Pivotamento Total

Neste caso identificamos o elemento de maior valor absoluto do sistema,

independente da coluna, e permutamos as colunas e em seguida permutamos as linhas. Um

cuidado que se deve ter, neste caso, é que ao permutarmos colunas as incógnitas também

são permutadas.

Ex 2 : Durante o processo de resolução de um sistema, no passo 2k , através do MEG

obteve-se o sistema:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 22

2 1 0 3 4

0 0 2 3 7

0 5 1 1 3

0 3 1 8 1

i) Realize um pivotamento parcial e

ii) Pivotamento total indicando a ordem das incógnitas.

Resolução...

---------------------------///---------------------------------

4.7 Decomposição LU

4.8 Método de Doolittle

Neste método veremos como escrever uma matriz como produto de duas matrizes,

uma triangular inferior “L” e outra triangular superior “U”, de modo que a diagonal de L seja

unitária, isto é, todos elementos iguais a um. Assim teremos:

11 12 1

21 22 2

1 2

n

n

n n nn

a a a

a a aA

a a a

=21

1 ( 1)

1 0 0 0

1 0 0

1 0

1n n n

l

l l

*

11 12 1

22 20

0 0

0 0 0

n

n

nn

u u u

u u

u

O algoritmo usado para obter LU neste método é dado por:

( 1)

1

( 1)

1

1 ; 1,2,...,

( . ) ; ,( 1),...,

( . )

; ( 1), (k+2) , ... ,

kk

k

kj kj kp pj

p

k

ik ip pj

p

ik

kk

se l k n

u a l u j k k n

a l u

l i k nu

Definição 1: ( ) 0; se b

n a

f n a b

Ex 1: Decompor a matriz A através do Método de Doolittle.

A =

5 2 1

3 1 4

1 1 3

Solução:...

---------------------------------------///------------------------------------------

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 23

Aplicação

Dado um sistema AX B podemos resolvê-lo através das matrizes LU.

Veja que se A LU LUX B , tomando UX Y termos

LY B encontramos o vetor Y e depois resolvemos UX Y encontrando assim a

solução X do sistema.

Ex 2: Resolva o sistema AX B , onde (0, 7, 5)tB e A é a matriz do exemplo

anterior.

4.9 Teorema LU

Se ( )ij n n

A a

e k

A o menor principal constituído pelas k primeiras e k primeira

colunas de A . Se det 0k

A para todo 1,2,3,..., 1k n .Então existe uma única

matriz L com 1kk

l triangular inferior e um única matriz U triangular superior tal que

A LU . Alem disso 11 22 33

det . . . ... .nn

A u u u u .

4.10 Matriz simétrica positiva definida

Definição: Uma matriz A real simétrica é positiva definida se para todos menores

principais k

A , vale det 0k

A com 1, 2, 3, ... , k n .

Corolário: Se A é simétrica positiva definida então A pode ser decomposta no produto tA GG , ondeG é uma matriz triangular inferior com elementos diagonais positivos.

4.11 Método de Cholesk

Se A é simétrica positiva definida então tA GG é tal que

ijG g então:

1

1 22

1

1

1

( ) 1,2,3,...,

1 , 2 ,...,

k

kk kk kpp

k

ik ip kpp

ik

kk

g a g k n

a g g

g i k k ng

Temos ainda que 2

det detA G .

Obs: O Método de Cholesk é utilizado também para verificar se a matriz é simétrica

definida positiva

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 24

Ex 1: Seja

4 2 4

2 10 4

4 4 9

A

i) Mostre que A satisfaz as condições para o Método de Cholesk.

ii) Decompor A em tA GG .

iii) Calcular o determinante de A através da decomposição do item anterior.

4.12 Inversão de Matrizes

Vamos considerar A LU com det 0A , isto é, A inversível. Então por

propriedade de matrizes temos 1 1 1A U L .

Se 1 1 Y L e Z U então a inversa de L e U pode ser encontrado através do

algoritmo: 1

1

( )

; , ( 1), ...,

( )

; , ( 1),...,1

1 se i j

0 se i j

i

ij ip pj

p j

ij

ii

j

ij ip pj

p i

ij

ii

ij

l y

y i j j nl

u z

z i j ju

Ex 1: Consideremos A LU dada por:

2 1 21 0 02 1 2

311 2 3 1 0 0 22 2

4 1 2 2 22 1 0 03 3

Encontre 1 1 1A U L .

Solução:...

4.13 Método de Eliminção de Gauss Compacto

O Método de Eliminação de Gauss pode ser utilizado para resolver equações

matriciais do tipo AX B , onde A e B são matrizes. Neste caso a única diferença é a

variação de j.

Assim o Método de Eliminação Compacto será:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 25

1, 2, 3, ... , (n-1);

; 0 2, 3, 4, ... , ;

1, 2, 3,..., 2

.

ikik kk

kk

ij ij ik kj

k para cada k faremosa

m a i n para cada i faremosa

j n

a a m a

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 26

CAPITULO 5

MÉTODOS ITERATIVOS PARA APROXIMAR SOLUÇÕES DE SISTEMAS LINEARES

5.1 Introdução

Em certos casos os métodos iterativos são melhores que os exatos, como por

exemplo, no caso de sistemas esparsos, isto é, a matriz dos coeficientes do sistema é

esparsa “possui um número grande elementos iguais a zero”, o que pode provocar um

número grande de pivotamento.

Os métodos iterativos ainda possuem a vantagem de que se um erro for cometido

eles se auto corrigem.

5.2 Método de Jacobi-Richardson (MJR)

Vamos considerar o sistema:

11 1 12 2 1 1

21 1 22 2 2 2

1 1 2 2

...

...

...

n n

n n

n n nn n n

a x a x a x b

a x a x a x b

a x a x a x b

Isolando i

x na i-ésima equação do sistema, teremos as relações do tipo

1 2 1 1( , ,... , ,..., )

i i i nx x x x x x

e a partir dessas criamos as regras de recorrência

semelhante ao que foi feito no MIL, assim teremos( 1) ( ) ( ) ( ) ( ) ( )

1 2 1 1( , ,... , , )

k k k k k k

i i i nx x x x x x

. Onde

os índices sobre-escritos indicam os passos das iterações.

Assim dado o sistema AX B teremos: ( 1) ( ) ( ) ( )

1 1 12 2 13 3 1

11

( 1) ( ) ( ) ( )

2 2 21 1 23 3 2

22

( 1) ( ) ( ) ( )

1 1 2 2 ( 1) 1

1[ ( ... )]

1[ ( ... )]

1[ ( ... )]

k k k k

n n

k k k k

n n

k k k k

n n n n n n n

nn

x b a x a x a xa

x b a x a x a xa

x b a x a x a xa

Neste caso o método será convergente quando o vetor erro tender ao vetor nulo,

isto é,

( 1) ( ) ( 1) ( ) ( 1) ( )

1 1 2 2, ,..., 0,0,...,0

k k k k k k

n nx x x x x xe

quando k tender a infinito. A

solução terá t casas decimais corretas quando o valor absoluto da componente máxima do

vetor erro for inferior a 0,5 10 t .

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 27

5.3_ Critérios de convergência

Quando isolamos i

x na i-ésima linha e criamos as regras de recorrência, não temos

garantia de convergência do método dado um chute inicial (0) (0) (0)

1 2( , ,..., )

nx x x , vejamos então

como garantir a convergência do método.

5.4 Critérios de linhas

Considere a matriz A dos coeficientes do sistema, se:

1

1max a , 1,2,3,..., 1

a

n

k ij

jkk

com j k e k n

então o MJR

converge independente do chute valor inicial.

Ex 1: Faça uma expanção do critério de linhas.

Ex 2: Encontre uma aprox. para a solução do sistema:

10 2 7

2 3 10 6

5 8

x y z

x y z

x y z

Através do MJR com k = 5 e diga qual a precisão do resultado.

Solução:...

------------------------------------///----------------------------

5.5 Método de Gauss-Siedel (MGS)

O MGS pode ser visto como MJR modificado, pois o que difere este agora é que na

iteração de cada equação já começamos a utilizar os resultados obtidos nas iterações

anteriores. Assim teremos:

( 1) ( ) ( ) ( ) ( )

1 1 12 2 13 3 14 4 1

11

( 1) ( 1) ( ) ( ) ( )

2 2 21 1 23 3 24 4 2

22

( 1) ( 1) ( 1) ( ) ( )

3 3 31 1 32 2 34 4 3

33

( 1)

1( ... )

1( ... )

1( ... )

k k k k k

n n

k k k k k

n n

k k k k k

n n

k

n

x b a x a x a x a xa

x b a x a x a x a xa

x b a x a x a x a xa

x

( 1) ( 1) ( 1) ( 1)

1 1 2 2 3 3 ( 1) 1

1( ... )

k k k k

n n n n n n n

nn

b a x a x a x a xa

No MGS se o critério de linhas não for satisfeito deveremos verificar outro critério

para a convergência.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 28

5.6 Critério de Sassenfeld

O MGS converge se tivermos o critério de linhas satisfeito ou se:

1

1 1

1max ( a j ) + a , 1,2,3,..., 1

a

i n

i ij ijj j i

ii

j i e i n

Obs 1: se o critério de linhas for satisfeito, podemos verificar que o critério de Sassenfeld

também será. ( Porque? ).

Obs 2: Caso os critérios não sejam satisfeitos podemos fazer permutações nas linhas e/ou

colunas do sistema com objetivo de satisfazer algum critério, mas a resolução devera ser

realizada a partir do sistema permutado e não do sistema original.

Ex 3: Faça uma expanção do Critério de sassenfel.

Ex 4: Resolva o sistema do ex anterior através do MGS com k = 5 e diga qual a precisão do

resultado.

5.7 Matriz estritamente diagonal dominante

Definição1

Uma matriz A e estritamente diagonal dominante se 1

1,2,...,n

ij iijj i

a a i n

.

Se tivermos 1

1,2,...,n

ij iijj i

a a i n

então A é dita diagonal dominante.

Obs1: se uma matriz A for estritamente diagonal dominante então o critério de

linhas e consequentemente o de Sassenfeld serão satisfeitos.

Obs 2: Se os critérios não forem satisfeitos, podemos modificá-los permutando

linhas e/ou colunas, o que não altera a solução, de modo que o sistema equivalente obtido

possua a matriz dos coeficientes estritamente diagonal dominante.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 29

CAPITULO 6

INTERPOLAÇÃO

Interpolar uma função consiste em encontrar uma função tal que a curva de sua

representação gráfica passe por um conjunto finito de pontos conhecidos, ou substituir uma

função, digamos complicada, por outra que seja mais simples de se trabalhar. A

interpolação é feita dentro de certo grau de precisão.

6.1 O problema da interpolação

Consideremos (n+1) pontos ( , ), 0,1, ... ,i i

x y i n chamados de pontos de

interpolação.

A interpolação tem por objetivo encontrar uma função ( )g x de tal modo que

0 0 1 1( ) ( ); ( ) ( ); ... ; ( ) ( )

n ng x f x g x f x g x f x .

6.2 Teorema de Weistrass (interpolação Polinomial)

Se f é uma função contínua em um intervalo [a,b] então, dado qualquer >0 existe

um polinômio P(x) tal que ( ) ( )P x f x .

6.3 Teorema

Dado um conjunto ( , ) 0,1,2,..., i i

F x y i n formado por (n+1) pontos,

então existe um único polinômio ( )n

P x de grau no máximo n tal que:

( ) ( ), 0,1,2,...,n i i

P x f x i n

6.4 Interpolação de Lagrange

Seja ( , ) 0,1,2,..., i i

f x y i n o conjunto de pontos de interpolação e

( )n

P x o polinômio de grau no máximo n que interpola f nesses pontos. Podemos

representar ( )n

P x da seguinte forma:

0 0 1 1 2 2( ) ( ) ( ) ( ) ( ) ( ) ( ) ... ( ) ( )

n n nP x L x f x L x f x L x f x L x f x

0

( ) ( ). ( )n

n k kk

P x L x f x

onde

1

0 0 1

( ) ( ) ( )( )

( ) ( ) ( )

n k nj j j

kj j j kj k k j k j k j

x x x x x xL x

x x x x x x

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 30

Ex 1: Considere a função tabelada abaixo:

x -2 0 1

y = f(x) 15 3 0

i) Encontre o polinômio ( )n

P x , através da int. de Lagrange, que interpole ( )f x .

ii) Estime o valor de (0.5)f .

Solução:...

--------------------------------///-------------------------------

Ex 2: Usando o Maple, para encontrar o polinomio interpolador que passa pelos (n+1)

primeiros pontos dos 10 pontos de uma função com valores tabelados. Obs: podemos

incluir mais pontos.

> restart:

> (x[0],x[1],x[2],x[3],x[4],x[5],x[6],x[7],x[8],x[9]):=(-

2,3,1,0,0,0,0,0,0,0);

> (y[0],y[1],y[2],y[3],y[4],y[5],y[6],y[7],y[8],y[9]):=(-

8,0,1,0,0,0,0,0,0,0);

> n:=3; #Grau maximo do pol. interpolador

> p:=expand(sum((product((x-x[j])/(x[k]-x[j]), j=0..k-

1)*product((x-x[j])/(x[k]-x[j]), j=k+1..n))*y[k],k=0..n));

> P:=unapply(p,x);

Ex 3: Usando o maple, dada a expressão da função a ser interpolada.

> restart:

> Digits:=20;

> f:=x->(sqrt((ln(x)+1)^(x+ln(x+2))));

> n:=4; # Grau maximo do polinomio interpolador #

> (x[0],x[1],x[2],x[3],x[4],x[5],x[6]):=(1.1,1.2,1.3,1.4,1.5,1.6,

1.7); # (n+1) valores de x do intervalo de interpolação #

> (y[0],y[1],y[2],y[3],y[4],y[5],y[6]):=(f(x[0]),f(x[1]),f(x[2]),

f(x[3]),f(x[4]),f(x[5]),f(x[6])); # (n+1) valores de y do

intervalo de interpolação #

> p:=expand(sum((product((x-x[j])/(x[k]-x[j]), j=0..k-

1)*product((x-x[j])/(x[k]-x[j]), j=k+1..n))*y[k],k=0..n));

> P:=unapply(p,x);# copie aqui a expressão obtida para o

polinimio#

> P(1.1); # digite os valores x'is para conferir #

_Podemos comparar o gráfico da função com o polinômio interpolador, no intervalo de interpolação.

> with(plots):

> g:=plot(p(x),x=1.1..1.5,y=0..2,color=blue):

> w:=plot(f(x),x=1.1..1.5,y=0..2,color=red ):

> display([g,w]);

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 31

Comparando o valor da integrais, teremos:

> int(f(x),x);

> int(P(x),x=1.1..1.5);

> int(f(x),x=1.1..1.5);

6.5 Erro na Interpolação

Quando estimamos um valor de ( )f x em um número x através da interpolação

temos como encontrar o erro máximo envolvido nessa estimativa, vejamos os resultados

abaixo.

6.6 Teorema

Se ( )f x for contínua em um intervalo [a, b], supondo que a derivada ( 1) ( )nf x

exista para todo x em [a, b].

Seja ( )n

P x o polinomio que interpola ( )f x sobre os (n+1) números distintos

0 1,...,

nx a x x b e suas respectivas imagens.

Então

0 1 ( 1)

0

( ).( ). ... .( )( ) ( )

( 1)!

n n

n n

x x x x x xE R x f onde x x

n

Obs: A existência desse valor é garantida, mas o problema é encontrar esse valor.

Para contornar isso temos o resultado.

6.7 Corolário

Seja ( ) ( ) ( )n n

R x P x f x se as derivadas até ordem (n+1) existem e são

contínuas em [a, b] então:

0 1 ( 1)

0

( ).( ). ... .( )( ) max ( )

( 1)!

n n

n n

x x x x x xE R x f t onde x t x

n

Ex 3: Calcule o erro máximo cometido ao calcular (0.52)f usando um 2( )P x ,

com ( ) 2 xf x xe onde foram usados os pontos de interpolação:

x 0,3 0,4 0,5 0,6 0,7 0,8

y= ( ) 2 xf x xe 0,81 1,2 1,65 2,19 2,82 3,56

6.8 Diferenças Divididas

Para o próximo método de interpolação usaremos as Diferenças Divididas, definidas

por:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 32

Se 0 1, , ... ,

nx x x são (n+1) números arbitrariamente espaçados e

0 1( ), ( ), ... , ( )

nf x f x f x suas imagens definimos:

[ ] ( ) 0,1,2, ... ,i i

f x f x i n

1 1 0

1 0

0

[ ,..., ] [ ,..., ][ ,..., , ] n n

n

n

f x x f x xf x x x

x x

Por diferenças divididas de ordem n da função f sobre os (n+1) pontos. Para facilitar

a visualização dos cálculos montamos uma tabela conforme ilustrado a seguir:

Tabela 01 - Diferenças Divididas

nx (0)ordem (1)ordem (2)ordem (3)ordem

0x

0( )f x

1 0

[ , ]f x x

1x

1( )f x

2 1 0[ , , ]f x x x

2 1

[ , ]f x x 3 2 1 0

[ , , , ]f x x x x

2x

2( )f x

3 2 1[ , , ]f x x x

3 2

[ , ]f x x 4 3 2 1

[ , , , ]f x x x x

3x

3( )f x

4 3 2[ , , ]f x x x

4 3

[ , ]f x x

4x

4( )f x

Ex 1:

0[ ]f x

1[ ]f x

1 2[ , ]f x x

2 1 0[ , , ]f x x x

E assim por diante.

6.9 Forma de Newton para o polinômio interpolador

Teorema 1

A forma de Newton para o polinômio ( )n

P x que interpola uma função ( )f x em

(n+1) pontos é dada

por:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 33

0 1 0 2 0 1 0 1 1( ) .( ) .( ).( ) ... .( ).( ). ... .( )

n n nP x d d x x d x x x x d x x x x x x

Onde 1 1 0

[ , ,..., , ] 0,1,2,...,k k k

d f x x x x com k n

Ex 1: Dada a função expressa pelos pares ordenados abaixo, construa a tabela de diferenças

divididas e o polinômio que interpola a função.

X -1 0 2 3

Y 4 1 -1 -2

Solução: ...

---------------------------///------------------------------

Teorema 2

Para [ , ], , 0,1,2,...,k

x a b x x k n então:

( 1)

1 0 0

( )[ ,..., , , ] ; ( , )

( 1)!

n

n n

ff x x x x x x

n

Em particular, a diferença dividida de ordem n de um polinômio ( )n

P x é independente do

ponto x e igual a n

a (coeficiente de grau n). Assim ao examinarmos uma tabela de

diferenças divididas de uma função, se as diferenças divididas de ordem k são praticamente

constantes, então as diferenças divididas de ordem k+1 são muito próximas de zero.

Podemos então usar um polinômio de grau k para interpolar tal função.

6.10 Erro na interpolação de Newton

Temos que o erro pode ser calculado por:

0 1 ( 1)

0

( ).( ). ... .( )( ) ( )

( 1)!

n n

n n

x x x x x xE R x f onde x x

n

Como não conhecemos então estimamos o erro fazendo:

0 1 ( 1)

0

( ).( ). ... .( )( ) max ( )

( 1)!

n n

n n

x x x x x xE R x f t onde x t x

n

Que também pode ser estimado, pelo teorema anterior, por:

0 1 1 1 0( ) ( ).( ). ... .( ) max [ ,..., , ]

n n nE R x x x x x x x f x x x

Ex 2: dada a função

X 2 3 4 5 6 8

Y 0.13 0.19 0.27 0.38 0.51 0.67

i) Determine qual o polinômio de grau mais adequado para interpolar valores de f.

ii) Calcular f(4.5) e estime o erro cometido.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 34

Quando os pontos bases 'ix s são igualmente espaçados, podemos aproveitar a

estrutura de dados para simplificar o processo de interpolação, manipulação e

implementação computacional.

6.11 Diferenças Ordinárias (ou Finitas)

Se os pontos as abssisas dos pontos de interpolação forem igualmente espaçados

podemos usar a estrutura anterior e simplificar os calculos.

Definição 1

Sejam 0 1, ,...

nx x x números igualmente espaçados, ou seja,

0ix x ih onde

1k kh x x

com 0,1,...,k n , do intervalo [a,b] e sejam também

0 1( ), ( ), ... , ( )

nf x f x f x suas respectivas imagens, defini-se:

0 ( 1) ( 1)( ) ( ) ( ) ( ) ( )r r r

k k k k kf x f x e f x f x h f x

Por diferenças ordinárias (ou finitas) de ordem r em k

x x . Com 0,1,...,k n

Ex 3: 0 0

0 0 1 1

1 0 0

0 0 0

1 0 0

0 1 0

1

0 1 0

( ) ( ); ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

( ) ( ) ( )

f x f x f x f x

f x f x h f x

f x f x f x

f x f x f x

E assim por diante.

Tabela 02 - Diferenças Ordinárias

nx (0)ordem (1)ordem (2)ordem (3)ordem

0x

0

0( )f x

1

0( )f x

1x

0

1( )f x 2

0( )f x

1

1( )f x 3

0( )f x

2x

0

2( )f x 2

1( )f x

1

2( )f x 3

1( )f x

3x

0

3( )f x 2

1( )f x

1

3( )f x

4x

0

4( )f x

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 35

Teorema 3

Se 0i

x x ih , isto é, se tivermos pontos igualmente espaçados, então:

0

1 0

( )[ ,..., , ]

. !

n

n n

f xf x x x

h n

.

Demonstração:

Usando indução finita ou matemática temos, para n = 1 1

0 1 0 0 01 0 1

1 0

( ) ( ) ( ) ( ) ( )[ , ]

.1!

f x f x f x h f x f xf x x

x x h h

Logo para n = 1 a afirmação do teorema é verdadeira.

Vamos supor que para um 1k n , 1n , qualquer ( 1)

01 1 0 ( 1)

( )[ ,..., , ]

.( 1)!

n

n n

f xf x x x

h n

Seja verdadeiro, logo devemos mostrar que para 1k n que a afirmação feita pelo

teorema também é verdadeira.

De fato:

2 1 1 1 01 0

0

[ ,..., , ] [ ,..., , ][ ,..., , ]

( )

n nn

n

f x x x f x x xf x x x

x x

, por definição.

( 1)( 1)

01

( 1) ( 1)

( )( )

.( 1)! .( 1)!

nn

n n

f xf x

h n h n

nh

, pela hipótese anterior.

( 1) ( 1)

0 0

( 1)

( ) ( )

. . .( 1)!

n n

n

f x h f x

n h h n

,

0( )

!

n

n

f x

h n

, simplificando e pela definição.

Logo para 1k n a afirmação também é verdadeira, logo pelo principio da

indução finita fica provado o Teorema.

6.12 Forma de interpolação de Newton-Gregory.

A forma de Newton para o polinômio interpolador, pelo teorema anterior sera:

1 20 0 0

0 0 0 11 2

00 1 1

( ) ( )( ) ( ) ( ) ( ).( ). ...

.1! .2!

( ) ( ).( ).....( ).

. !

n

n

n n

f x f xP x f x x x x x x x

h h

f xx x x x x x

h n

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 36

Ex 1: Considere a função tabelada abaixo:

X 0,3 0,4 0,5 0,6 0,7

Y 0.5466 0.8902 1.3591 1.9921 2.3386

Encontre através da formula de Newto-Gregory um polinômio quadrático para

estimar o valor de (0,52)f e determine o erro máximo cometido.

Solução:...

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 37

CAPITULO 7

DIFERENCIAÇÃO NUMÉRICA

7.1 Diferenciação por Diferenças Finitas

Se f é uma função contínua em ( , )a b temos por definição que a derivada dessa

função é dada por:

0

( ) ( )limh

f x h f x

h

com ( , )x a b

Assim se h for pequeno sulficiente podemos estimar a derivada de f em ponto

( , )x a b por:

( ) ( )'( )

f x h f xf x

h

Ex 1: Dada a função ( ) ln( )f x x ,estime o valor de '(1.8)f com 0.01h e

0.0001h e compare com o valor obtido na calculadora.

Solução:...

-----------------------------------///-------------------------------------

Vamos desenvolver o método através da interpolação para estimar o valor das

derivadas.

Tomando 3 pontos e encontrando o polinômio interpolador de Lagrange teremos:

2

20

( ) ( ). ( )k k

k

P x L x f x

onde 2

0

( )( )

( )

k

kj

k jj k

x xL x

x x

e assim

2

0

' ( ) ' ( ). ( )n k k

k

P x L x f x

e assim:

1 2 2 10 0

0 1 0 2 0 1 0 2

0 1 0 21 1

1 0 1 2 1 0 1 2

0 1 0 12 2

2 0 2 1 2 0 2 1

( ).( ) 2( ) ' ( )

( ).( ) ( ).( )

( ).( ) 2( ) ' ( )

( ).( ) ( ).( )

( ).( ) 2( ) ' ( )

( ).( ) ( ).( )

x x x x x x xL x L x

x x x x x x x x

x x x x x x xL x L x

x x x x x x x x

x x x x x x xL x L x

x x x x x x x x

Resultara em:

2 1 0 2 0 10 1 2

0 1 0 2 1 0 1 2 2 0 2 1

2 2 2'( ) ( ) ( ) ( ) (*)

( ).( ) ( ).( ) ( ).( )

k k kk

x x x x x x x x xf x f x f x f x

x x x x x x x x x x x x

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 38

Vamos considerar agora 5 pontos igualmente espaçados

2 , , , , 2 .k k k k kx h x h x x h x h conforme gráfico:

Veja que podemos tomar 3 pontos de 3 diferentes formas, com o objetivo de calcular

'( )kf x :

7.2 Derivadas por diferenças finitas retroativas

Essa derivada em um número kx , é obtida se tomarmos os 3 primeiros pontos, dos

5, assim:

0 1 2( 2 ), ( ), ( )k k kx x h x x h x x e também:

0 1 0 2 1 0 1 2

2 0 2 1

( ) , ( ) 2 , ( ) , ( ) ,

( ) 2 , ( )

x x h x x h x x h x x h

x x h x x h

Substituindo em (*) teremos:

1

'( ) ( 2 ) 4 ( ) 3 ( )2

k k k kf x f x h f x h f xh

7.3 Derivadas por diferenças divididas centrais.

Se tomarmos os 3 pontos centrais, de forma equivalente teremos:

1

'( ) ( ) ( )2

k k kf x f x h f x hh

7.4 Derivadas por diferenças finitas progressivas.

E tomando os 3 últimos pontos, de forma equivalente teremos

1

'( ) ( 2 ) 4 ( ) 3 ( )2

k k k kf x f x h f x h f xh

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 39

CAPITULO 8

APROXIMAÇÃO DE FUNÇÕES

8.1 Método dos Mínimos Quadrados ( MMQ )

O objetivo do MMQ é substituir uma função ( )f x por outra ( )x , ou encontrar

uma função que melhor se aproxime de um conjunto de pontos, mas, levando em

consideração o comportamento destes, e não o total de pontos como na interpolação,

médoto esse muito usado no processo de modelagem.

Neste caso escreveremos a função procurada como combinação linear de funções

conhecidas. Assim teremos:0 0 1 1( ) ( ) ( ) ( ) ... ( )n nf x x a g x a g x a g x onde

as funções ( ) 0,1,...,ig x i n são conhecidas. A substituição devera ser feita de modo

que ( ) ( )f x x seja mínima. Isto é devemos minimizar a soma dos erros dado por:

0

( ) ( )n

i i

i

E f x x

A substituição da f é indicada quando o uso da função oferece alguns

inconvenientes tais como:

i) f é dada por processos não finitos, como soma de séries, integrais etc.

ii) f é dada por um conjunto de pares ordenados, obtidos experimentalmente e

queremos encontrar uma função que melhor modele o fenômeno estudado.

8.2 Aproximação Polinomial caso discreto

Vamos considerar um caso discreto, onde queremos aproximar um conjunto de

pares ordenados ( , ) 0,1,...,i ix y i m por um polinômio, que podemos escrever

como: 1 2

0 1 2

0

( ) ...

=

n

n n

nk

k

k

P x a a x a x a x

a x

de grau n < m, através do MMQ.

Neste caso devemos minimizar

0

| ( ) |n

i n i

i

E y P x

o problema que surge, que essa

soma sempre será nula, sem que função passe exatamente pelos pontos.

Para contornar esse problema vamos então minimizar a soma do quadrado desses

erros, que nunca se anularão, e se o quadrado for mínimo, a soma dos erros também serão,

o que justifica o nome de Métodos dos Mínimos Quadrados.

Assim Teremos:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 40

2

0

[ ( )]n

i n i

i

E y P x

Para encontrar o polinômio, basta encontrar os seus coeficientes , 0,1,...,ja j n , logo

vamos a busca dos números ja tal que 0j

E

a

. O que nos dara:

0 0 0

0 0 0 0 0

0 0 0 0 0 0

0

( )2( ( )). 0 ( ).( ) 0

0 0

0

m m nk ji

i n i i k i i

i i kj

m n m m nj k j j k j

i i k i i i k i

i k i i k

m m n m n mj k j k j j

i i k i k i i i

i i k i k i

nk

k i

k

P xy P x y a x x

a

y x a x y x a x

y x a x a x y x

a x

0 0

; 0,1,2,...,m m

j j

i i

i i

y x j n

0 0

0 0 0

0 1 2 0

0 1 2

0 0 0 0 0

0

:

0 :

;

... (1)

1 :

n m mk

k i i i

k i i

m m m m mn

i i i n i i i

i i i i i

nk

k i

k

Atribuindo valores para j e Expandindo as somatórias teremos

Para j temos a equação

a x y x

a x a x a x a x y x

Para j temos a equação

a x

1 1

0 0

1 2 3 1 2

0 1 2

0 0 0 0 0

1 2 2

0 1 2

0 0 0 0 0

;

... (2)

:

... ( )

m m

i i

i i

m m m m mn

i i i n i i i

i i i i i

m m m m mn n n n n

i i i n i i i

i i i i i

y x

a x a x a x a x y x

Para j n teremos

a x a x a x a x y x n

Note que temos n equações, chamadas de equações normais, e n incógnitas, logo

para descobrir o valor dos coeficientes do polinômio, basta resolver o sistema escrito na

forma matricial abaixo;

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 41

0 1 00

0 0 0 0

1 2 1 1

1

0 0 0 0

1 2

0 0 0 0

...

...

m m m mn

i i i i i

i i i i

m m m mn

i i i i i

i i i i

m m m mn n n n

i i i i ini i i i

ax x x x y

x x x a x y

x x x x ya

Cuja solução será os coeficientes procurados dos polinômios. Note que a matriz dos

coeficientes, somatórias, desse sistema é simétrica.

Neste caso vamos considerar 0 1x para todo x , o que ocorre na verdade apenas

se x não for nulo.

Ex 1: Considere o conjunto de pares ordenados abaixo:

x 0 0,25 0,5 0,75 1

y 1 1,2840 1,6487 2,1170 2,7183

i) determine através do MMQ o polinômio quadrático que melhor se aproxime da função

tabelada.

ii) Estime o valor de (0.3)f .

8.4 Aproximação Polinomial caso continuo

Vamos supor que tenhamos uma função f explicita por uma sentença de x , e não

apenas por um conjunto de pares ordenados e vamos supor que queremos substituir f por

outra função , em um intervalo [a,b], que tenha uma definição mais fácil de se trabalhar.

Graficamente teríamos:

Neste caso tevemos encontrar

0 0 1 1 2 2( ) ( ) ( ) ( ) ... ( )n nx a g x a g x a g x a g x de modo que a área

destacada seja mínima, o que é equivalente a minizar:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 42

2

( ) ( )b

aE f x x dx

fazendo

0j

E

a

; j = 0, 1, 2, ... , n

2 ( ) ( ) . ( ) 0b

ja

j

Ef x x g x dx

a

( ) ( ) . ( ) 0b

ja

f x x g x dx

Substituindo 0 0 1 1 2 2( ) ( ) ( ) ( ) ... ( )n nx a g x a g x a g x a g x , na equação

anterior teremos:

0 0 1 1 2 2( ) ( ( ) ( ) ( ) ... ( )) . ( ) 0b

n n ja

f x a g x a g x a g x a g x g x dx

Aplicando a propriedade distributiva e as propriedades de integral teremos

0 0

0 0

( ) ( ) ( ) ( ) ... ( ) ( ) 0

( ) ( ) ... ( ) ( ) ( ) ( )

b b b

j j n n ja a a

b b b

j n n j ja a a

f x g x dx a g x g x dx a g x g x dx

a g x g x dx a g x g x dx f x g x dx

Para cada valor de j = 0, 1, ... ,n, assim como no caso discreto teremos uma equação,

resultando no sistema (nxn) escrito na forma matricial abaixo:

0 0 1 0 0 0

0

10 1 1 1 1 1

0 1

b b b b

na a a a

b b b b

na a a a

b b b bn

n n n n na a a a

g g dx g g dx g g dx f g dxa

ag g dx g g dx g g dx f g dx

ag g dx g g dx g g dx f g dx

Obs: Para simplificar a notação estamos considerando ( )g g x .

Considerando um espaço vetorial euclidiano, e o produto vetorial dado por:

. ( ) ( )b

af g f x g x dx

podemos escrever o sistema acima como:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 43

0 0 1 0 0 0 0

0 1 1 1 1 1 1

0 1

, , , .

, , , .

, , , .

n

n

n n n n n n

g g g g g g a f g

g g g g g g a f g

g g g g g g a f g

, ( ).

b

a

f g f x g x dx

A solução do sistema serão os coeficientes do polinômio.

Ex 1: Seja 4( ) 5f x x x aproximar ( )f x por um polinômio de grau dois,o itervalo [-

1,1] através do MMQ.

8.5 Aproximações não Lineares

Em casos onde a função na qual se queira obter a aproximação não seja do tipo:

0 0 1 1( ) ( ) ( ) ( ) ... ( )n nf x x a g x a g x a g x

Podemos fazer uma adaptação para cair em um caso como este.

Ex: vamos supor que se queira obter um modelo do tipo xy e a partir de um conjunto

de pares ordenados. Note que

ln( ) ln( ) ln( )

ln( ) ln( )

ln( ) (1)

x x

x

y e y e se z y

z e

z x

Chamando 0ln( ) a e

1a teremos, 0 1z a a x que neste caso é uma

linearização do problema inicial.

Note que 0 1z a a x é a equação de uma reta e temos

0 ( ) 1g x , 1( )g x x e

caímos em um caso linear visto anteriormente. Uma vez encontrada a equação

0 1z a a x basta lembra que:

0 1

0 1

0 1

0 1

ln( ) .

. (2)

a a xz

a a xy

a a x

z a a x e e

e e e

y e e

Comparando (1) e (2) temos:

0 1.a a x xy e e e

Logo temos que 0ae e

1a

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 44

8.6 Teste de Alinhamento

Uma vez escolhida uma função não linear para ajustar uma função linear dada, uma

forma de verificarmos se a escolha feita foi razoável é aplicarmos o teste do alinhamento,

que consiste em:

i) fazer a linearização da função não linear escolhida;

ii) fazer o diagrama de dispersão dos novos dados;

iii) se os pontos do diagrama ii) estiverem alinhados, isto significará que a função não

linear escolhida foi uma boa escolha.

Obs: devidos a erros de observação e cálculos aproximados consideremos

satisfatório o diagrama de dispersão onde os pontos se distribuem aleatoriamente em torno

de uma reta.

Ex 2: Faça um teste de alinhamento e encontre uma função do tipo xy e para

modelar a função tabelada abaixo:

x -1 -0.7 -0.4 -0.1 0.2 0.5 0.8 1

y 36.547 17.264 8.155 3.852 1.820 0.860 0.406 0.246

Solução:...

Dica: faça uma pesquisa sobre outros casos não lineares na apostila da Neide Maria

Berttoldi Franco.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 45

CAPITULO 9

INTEGRAÇÃO NUMÉRICA

Formulas fechadas de Newton-Cotes

9.1 Regra dos trapézios

Vamos supor uma função ( )f x contínua e posítiva em um intervalo [ , ]a b .

Na regra dos trapézios, ao invés de dividirmos a área abaixo da curva ( )y f x ,

acima do eixo x , e entre as retas 0 e nx a x b em retângulos, dividimos-a em

trapézios.

Dividindo o intervalo [ , ]a b em n subintervalos cujo os extremos serão:

0 1 2, , , ... , nx a x x x b , isto é, em subintervalos do tipo

1[ , ] com 1,2,3,...,i ix x i n . Conforme ilustra o gráfico abaixo:

Considerando um trapézio qualquer com as alturas medindo 1( ) e ( )i if x f x

e

base 1i ix x x sua área pode ser calculada por: 1( ) ( )

2

i if x f xA x onde

1i ix x pode ser calculado da seguinte forma b a

xn

onde n é o número de

subintervalos de [ , ]a b . Obs: as vezes usamos h ao invés de x

Para obtermos uma aproximação para área total, isto é, uma aproximação para o

valor da integral neste caso, somamos as áreas dos trapézios. Assim teremos:

0 1 2 31 2

2 1 1

( ) ( ) ( ) ( )( ) ( )( ) ...

2 2 2

( ) ( ) ( ) ( )

2 2

b

a

n n n n

f x f x f x f xf x f xf x dx x x x

f x f x f x f xx x

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 46

0 1 1 2 2 1 1

0 1 2 2 1 1

( ) ( ) ( ) ( ) ( ) ... ( ) ( ) ( ) ( )2

( ) ( ) 2 ( ) 2 ( ) ... 2 ( ) 2 ( ) 2 ( ) ( )2

b

n n n na

b

n n n na

xf x dx f x f x f x f x f x f x f x f x

xf x dx f x f x f x f x f x f x f x

1

0

1

( ) ( ) ( ) 2 ( )2

nb

n ia

i

xf x dx f x f x f x

É possível mostrar que o erro máximo cometido na regra dos trapézios é dado por:

3

0

( ).max ''( ) ;

12n

n xE f t x t x

Ex 1: Estime o valor da integral:

4

21 16

dx

x através da regra dos trapézio tomando n = 6.

9.2 Regra 1/3 de Simpsom

Teorema 1

Se 0 0 0 1 1 1 2 2 2( , ); ( , ); ( , )P x y P x y P x y forem 3 pontos distintos sobre a parábola

com equação 2y Ax Bx C come

0 1 2, , 0y y y e

1 0 2 0; 2x x h x x h , ou seja, igualmente espaçados.

Então a medida da área da região limitada pela parábola, pelo eixo x e pelas retas

0 2; x x x x será da por:

0 1 243

hA y y y

Assim se ( )y f x for continua e positiva em um intervalo [ , ]a b podemos dividir

[ , ]a b em n subintervalos conforme feito na regra dos trapézios e tomar arcos de

parábola, e usando idéia semelhante teremos:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 47

Grafico

0 1 2 2 3 4

2 1

0 1 2 3 4

2 1

( ) ( ( ) 4 ( ) ( )) ( ( ) 4 ( ) ( )) ...3 3

( ( ) 4 ( ) ( ))3

( ) ( ( ) 4 ( ) 2 ( ) 4 ( ) 2 ( ) ...3

2 ( ) 4 ( ) ( ))

b

a

n n n

b

a

n n n

h hf x dx f x f x f x f x f x f x

hf x f x f x

hf x dx f x f x f x f x f x

f x f x f x

12 2

0 2 1 2

1 1

( ) ( ) ( ) 4 ( ) 2 ( )3

n nb

n i ia

i i

hf x dx f x f x f x f x

O erro cometido pela regra 1/3 de Simpson pode ser calculado por:

5(4)

0.max ( ) ; 90

n

nhe f t x t x .

Ex 1: Estime o valor da integral:

4

21 16

dx

x através da regra da regra 1/3 de Simpson com n = 4.

9.3 Polinômios ortogonais e Quadratura de Gauss

Definição 1

Vamos considerar o produto escalar:

( ), ( ) ( ) ( ) ( )b

af x g x w x f x g x dx

Onde ( ) 0w x contínua ( , )a b , e ( ) w x é dita função peso.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 48

Definição 2

Seja 0 1 2, , ,... uma família de polinômios de graus 0,1,2...

Se

, 0 se

, 0 se 0

i j

i i i

i j

Então os polinômios 0 1 2, , ,... se dizem ortogonais.

Teorema 1

Sejam os polinômios 0 1 2, , ,... definidos por:

0 ( ) 1x

0 0

1

0 0

. ( ), ( )

( ), ( )

x x xx

x x

e para 1,2,3,...k

1 1( ) ( ) ( )k k k kx x x x

1 1

. ( ), ( ) ( ), ( );

( ), ( ) ( ), ( )

k k k k

k k

k k k k

x x x x x

x x x x

Os polinômios 0 1 2, , ,... , assim definidos são dois a dois ortogonais.

9.4 Principais Polinômios Ortogonais

i) Polinômio de Legendre

Os polinômios 0 1 2, , ,...P P P , são obtidos segundo o produto escalar:

1

1( ) ( )f x g x dx

Neste caso temos ( ) 1w x e [ , ] [ 1, 1]a b .

ii) Polinômios de Tchebyshev

O produto escalar usado para obter os polinômios de Tchebyshev 0 1 2, , ,...T T T é

dado por:

1

21

1( ) ( )

1f x g x dx

x ou seja neste caso temos

2

1( )

1w x

x

com

[ , ] [ 1, 1]a b

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 49

iii) Polinômios de Laguerre

Para obter os polinômios de Laguerre usamos o produto interno:

0

( ) ( )xe f x g x dx

assim temos ( ) xw x e e [ , ] [0, )a b .

iv) Polinômios de Hermite

O produto interno neste caso será:

2

( ) ( )xe f x g x dx

assim temos 2

( ) xw x e e [ , ] ( , )a b .

Ex 1: Obter os 3 primeiros polinômios de Legendre.

Solução:...

---------------------------------///---------------------------------

9.5 Propriedade dos Polinômios Ortogonais

i) Todo polinômio de grau menor ou igual a n, pode ser escrito como combinação

linear de polinômios ortogonais 0 1 2, , ,..., n .

ii) Sejam 0 1 2, , ,..., n polinômios ortogonais.

Então n é ortogonal a qualquer polinômio ( )Q x de grau menor ou igual a n.

iii) Sejam 0 1 2, , ,..., n polinômios ortogonais segundo o produto escalar:

( ), ( ) ( ) ( ) ( )b

af x g x w x f x g x dx com ( ) 0w x e continua em ( , )a b . Então

n possui n raízes reais, distintas em ( , )a b .

iv) Sejam 0 1 2, , ,..., n nas condições da prop. 3).

Sejam 0 1 2, , ,..., nx x x x as raízes de 1n

. Se é um polinômio de grau menor ou igual a

2 1n então:

0

( ) ( ) ( )nb

k ka

k

w x f x dx A f x

onde ( ) ( )b

k ka

A w x L x dx com ( )kL x da

interpolação de Lagrange.

9.6 Fórmula da Quadratura de Gauss

São as formulas usadas para calcular:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 50

( ) ( )b

aw x f x dx valendo-se da propriedade 4. calculamos o valor aproximado da integral

usando:

0

( ) ( ) ( )nb

k ka

k

w x f x dx A f x

onde ( ) ( )b

k ka

A w x L x dx com ( )kL x

da interpolação de Lagrange, sobre os pontos 0 1 2, , ,..., nx x x x que são as raízes do

polinômio ortogonal 1n

.

Assim o procedimento para se calcular uma integral usando quadratura de Gauss, é o

seguinte:

a) Determinar o polinômio ortogonal 1n

,segundo o produto escalar conveniente,

isto é, a função peso ( )w x e no intervalo [ , ]a b .

b) Calcular as raízes 0 1 2, , ,..., nx x x x de

1n .

c) Determinar os polinômios de Lagrange ( )kL x usando os pontos obtidos em b).

d) Calcular ( ) ( )b

k ka

A w x L x dx , 0,1,2,...k n .

e) Calcular o valor de ( ); 0,1,2,...kf x k n .

f) Calcular finalmente

0

( ) ( ) ( )nb

k ka

k

w x f x dx A f x

.

Veja que é um processo trabalhoso, mas para facilitar, a partir do momento em que

se define o polinômio ortogonal a ser usado, podemos usar valores já calculados e

tabelados de suas raízes e os valores de kA também já calculados e tabelados, assim só

será necessário realizar os itens e) e f).

9.7 Fórmula de Gauss-Legendre

Para usarmos a formula de Gauss-Legendre a integral a ser calculada deve ter função

peso ( ) 1w x e os limites de integração devera ser o interval [-1, 1], caso os limites

de integração não sejam estes, estão deve-se fazer uma mudança de variável para realizar a

integração numérica.

Ex 1: Estime o valor da integral, usando quadratura de Gauss indicada.

i)1

3 2

1( 5 )x x dx

Gauss-Legendre.

ii) 2

1

1dx

x Gauss-Legendre.

9.8 Fórmula de Gauss-Tchebyshev

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 51

Neste caso a integral a ser calculada deverá ter a função peso

2

1( )

1w x

x

com [ , ] [ 1,1]a b .caso o intervalo de integração seja outro, deve-

se fazer uma mudança de variável.

Ex 2: use a formula de Gauss-Thebyshev para calcular:

3 22

22

2

4 4

t tdt

t

9.9 Fórmula de Gauss-Laguerre

Aqui a integral, imprópria, a ser calculada deverá ter função peso ( ) xw x e e

intervalo de integração [ , ] [0, )a b . Caso o intervalo seja outro deveremos tentar

mudar de variável.

Ex 3: Use a formula de Gauss-Laguerre para calcular a integral:

33

1

1dx

x

9.10 Fórmula de Gauss-Hermite

Neste deveremos ter a integral imprópria com função peso 2

( ) xw x e e intervalo

de integração ( , ) .

Ex 4: Calcule:

2( 1)xe dx

9.11 Erros nas Fórmulas de Gauss

Quando a função a ser integrada é um polinômio sabemos que as formulas de

quadraturas fornecem o valor exato da integral.

Caso contrario o método fornece um valor aproximado, onde os erros podem ser

estimados resolvendo:

9.12 Erro da formula de Gauss-Legendre

2 3 4(2 2)

3

2 [( 1)!]max ( ) ;

(2 3).[(2 2)!]

n

nnE f t a t b

n n

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 52

9.13 Erro da formula de Gauss-Tchebyshev

(2 2)

(2 2)

2max ( ) ;

(2 2)!

n

nE f t a t b

e n

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 53

CAPITULO 10

MÉTODOS NUMÉRICOS PARA APROXIMAR SOLUÇÕES DE EQUAÇÕES

DIFERENCIAIS ORDINARIAS

Como sabemos para encontrarmos a solução analítica de uma EDO as vezes o

processo é um tanto trabalhoso e até impossível em certos casos, assim veremos aqui

métodos numéricos que permite descobrir valores numéricos aprox. da solução, mesmo

sem que se conheça a solução analítica da EDO.

Sempre partiremos de um valor conhecido, que é a condição inicial, neste caso

temos um “p.v.i” problema de valor inicial.

Ex 1:

'

(ln(2)) 1

y y

y

Solução:...

De modo mais geral, vamos trabalhar com problemas do tipo:

0 0

' ( , )

( )

y f x y

y x y

Onde os valores 0 0 e x y serão valores conhecidos, no caso, a condição inicial.

CUIDADO: 0x nem sempre é zero, pode ser um valor qualquer.

10.1 Método de Taylor de ordem q

Dividindo um intervalo [ , ]a b em subintervalos, com extremos

0 1 2, , ,..., nx a x x x b igualmente espaçados, a uma distancia h , e expandindo a

função 1( ) ( )n ny x y x h em torno do ponto

nx teremos:

2 ( 1) ( 1)

1

'( ) ''( ) ( ) ( )( ) ( ) ...

1! 2! ! ( 1)!

q q q q

n n nn n

hy x h y x h y x h yy x y x

q q

(*)

Onde o ultimo termo é o erro de truncamento local, com 1n nx x .

Se f for suficientemente diferençiável, as sua derivadas podem ser obtidas,

considerando ' ( , )y f x y como uma função só x de pois ( )y y x , o que nos

permite pensar na derivada exata de ( , )f x y em termos de x . Pela regra da cadeia

teremos:

' ( , ) ( , )dy

y f x y f x ydx

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 54

( , ) ( , )'' '( , )

f x y dx f x y dyy f x y

x dx y dx

'( , ) ( , ) ( , ) ( , )x yf x y f x y f x y f x y

'( , ) '( , )''' ''( , )

f x y dx f x y dyy f x y

x dx y dx

''( , ) ' ( , ) ' ( , ) ( , )x yf x y f x y f x y f x y

E assim por diante.

Para simplificar a notação considerar daqui por diante:

( )n ny x y , ( ) ( )( , ) 0,1,2,...,j j

n n nf x y f j n onde (0)

n nf f .

Truncando a equação (*) em que q termos, e substituindo as derivadas de f já na

nova notação teremos: 2 ( 1)

1

'...

1! 2! !

q q

n n nn n

hf h f h fy y

q

conhecido como Método de Taylor de

ordem q

Ex 1: Resolver o p.v.i através do método de Taylor de ordem 3. 2' 2 1

(0) 1

y y x

y

No intervalo [0; 3] 0.1com h

10.2 Métodos dos Passos Múltiplos

10.3 Método de Euller

Se tomarmos o método de Taylor de ordem 1, teremos o que também é chamado de

Método de Euller, Mas a dedução do método pode ser feita de varias formas.

Se ' ( , )y f x y como ( )y y x podemos pensar em f apenas como uma função

( )f x o que nos permite algumas deduções.

Vejamos uma dedução geométrica do método:

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 55

Avaliando as áreas. Teremos:

1 1

1' ( ) ( )n n

nn

x x

n nxxI y dx y y x y x

por outro lado, neste caso, temos

1( ). ( ) . ( )n n n nI x x f x h f x o que nos dara

1 1( ) ( ) . ( ) ( ) ( ) ( )n n n n n ny x y x h f x y x y x hf x

então se tivermos

1 0n nx x h

lembrando que ( ) ( , ( ))f x f x y x assim podemos ver que

1( ) ( ) ( , )n n n ny x y x hf x y

que é o Método de Euller.

Escrevendo na notação adotada teremos:

O erro do método de Euller pode ser calculado por

(2) '( )

2!

h f com

1n nx x

Note ser h não for muito pequeno então o erro será grande. Para contornar isso é

feito um aprimoramento do método de Euller, usando a idéia de áreas, aplicando a regra

dos trapézios teremos:

10.4 Método de Euller Aprimorado

Da mesma forma temos:

1n n ny y hf

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 56

1 1

1' ( ) ( )n n

nn

x x

n nxxI y dx y y x y x

por outro lado, neste caso, teremos:

11 1

[ ( ) ( )]( ). .[ ( ) ( )]

2 2

n nn n n n

f x f x hI x x f x f x

que resultara em

1 1 1 1( ) ( ) .[ ( ) ( )] ( ) ( ) .[ ( ) ( )]2 2

n n n n n n n n

h hy x y x f x f x y x y x f x f x

Se tivermos, como no exemplo anterior 1 0n nx x h então

1 1( ) ( ) .[ ( ) ( )]2

n n n n

hy x y x f x f x como estamos considerando

( ) ( , )f x f x y então: 1 1 1( ) ( ) .[ ( , ) ( , )]

2n n n n n n

hy x y x f x y f x y . O

problema aqui é que temos 1( )ny x

em ambos os lados da equação, assim deveríamos

isolar 1( )ny x

o que nem sempre é fácil, as vezes impossível, tomaremos então

1( ) ( ) ( , )n n n ny x y x hf x y do método de Euller. O que resultara

1 1( ) ( ) . ( ) , ( )2

n n n n n n

hy x y x f x f x y hf x que é chamado de Método

de Euller Aprimorado. Podendo ser escrito como:

1 1.2

n n n n

hy y f f

Ex 1:

Encontre valores para 1y e

2y através do método de Euller aprimorado com 0.025h

Onde:

' 1 4 (0) 1y x y y

Solução:

----------------------------------------///--------------------------------------

Usando o Maple podemos fazer os cálculos da seguinte forma:

Método de Euller Aprimorado. > restart: Digits:=6;

f:=(x,y)->3*x^2-7*x+3-3*y;

h:=0.1;

x[0]:=1; y[0]:=2; n:=0;

> while n<2 do

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 57

x[n+1]:=x[0]+(n+1)*h;

("x"=x[n], "y"=y[n]);

k[1]:=evalf(f(x[n],y[n]));

("x"=x[n]+h, "y"=y[n]+h*k[1]);

k[2]:=evalf(f(x[n]+h,y[n]+h*k[1]));

y[n+1]:=evalf(y[n]+h/2*(k[1]+k[2]));

print("----------------(((o)))------------------");

n:=n+1:

end do;

10.5 Erros nos método de Euller e Euller aprimorado

Como o Método de Euller pode ser interpretado como método de Taylor de ordem q

erro pode ser estimado por:

2

1

'( , ( ))max ;

2n n

h f t y tE x t x

Note que deduzimos o Método de Euller Aprimorado, através da regra dos trapézios,

assim o seu erro pode ser estimado por:

3

1

''( , ( ))max ;

12n n

h f t y tE x t x

Podemos verificar que no Método de Euller tomarmos 2

hh então o erro será:,

2

2'( , ( )) '( , ( ))max max

2 4 . 2

h f t y t h f t y tE ;

1n nx t x assim com

2

hh o erro foi reduzido em 4 partes.

Já no Método de Euller Aprimorado, podemos verificar facilmente que se tomarmos

2

hh então o erro será reduzido em 8 vezes.

O que torna o tamanho do passo h importante para controlar o erro dentro de

parâmetros aceitáveis.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 58

10.6 Método de Runge-Kutta

Definição 1

O método geral de Runge-Kutta de R estágios é definido por:

1 ( , , )n n n ny y h x y h onde

1

( , , )R

n n r r

r

x y h c k

1 ( , )k f x y

1

1

( , ) 2,3,4,...,r

r r rs s

s

k f x a h y h b k r R

1

1

; 2,3,4,...,r

r rs

s

a b r R

10.7 Método de Runge-Kutta de ordem 2

Nesta caso temos então R = 2, assim:

2

1 1 2 2

1

( , , )n n r r

r

x y h c k c k c k

1 ( , )k f x y

1

2 2

1

( , ) rs s

s

k f x a h y h b k

2 2 21 1 2 2 1( , ) ( , ) k f x a h y hb k f x a h y a hk pois

2 21a b

Para determinar as constantes 1 2 2, , c c a , que omitirei os cálculos aqui, mas estes

podem ser encontrados nas referências citadas no final desse trabalho, desenvolvendo

2k em em série de Taylor torno do ponto ( , )x y e desenvolvendo ( , , )n nx y h em série

Taylor, e comparando os resultados obtidos teremos o sistema:

1 2

2 2

1

1

2

c c

c a

Que possui infinitas soluções, logo atribuindo valores para uma das incógnitas,

encontramos os valores das demais em função desta.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 59

Um método de Runge-Kutta de ordem 2 e 2 estágios mais conhecido é obtido

atribuindo os valores:

1 2 2

1 1 1

2 2c c e a

1 1 2

1 1

2 2n ny y hk hk

1 1 2

1

2n ny y h k k

Onde: 1 ( , )k f x y e

2 1( , ) k f x h y hk conhecido também como

Método de Euller Aprimorado.

10.8 Método de Runge-Kutta de 4 estágios e ordem 4

É o mais utilizado na pratica e é dado por:

1 1 2 3 42 26

n n

hy y k k k k

Onde:

1

2 1

3 2

4 3

( , )

( , )2 2

( , )2 2

( , )

n n

n n

n n

n n

k f x y

h hk f x y k

h hk f x y k

k f x h y hk

Ex 45: Encontre as soluções do p.v.i no intervalo [0,0.4]com 0.2h e (0) 1y

' 1 4y x y

Solução;...

-------------------------------///----------------------------------

A mesma resolução usando o Maple. Método de Runge-Kutta >

> restart:

> f:=(x,y)->1-x-4*y

> h:=0.2;# tamanho do passo

> x[0]:=0; y[0]:=1; # cond. iniciais

> n:=0;

> while n < 2 do

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 60

x[n+1]:=x[0]+(n+1)*h;

k[1]:=f(x[n],y[n]);

k[2]:=f(x[n]+1/2*h,y[n]+h/2*k[1]);

k[3]:=f(x[n]+1/2*h,y[n]+h/2*k[2]);

k[4]:=f(x[n]+h,y[n]+h*k[3]);

print("---------//-----------");

y[n+1]:=y[n]+h/6*(k[1]+2*k[2]+2*k[3]+k[4]);

print("---------//-----------");

n:=n+1:

end do;

10.9 Sistemas de Equações Diferenciais Ordinárias

Consideremos um sistema de n equações diferenciais de primeira ordem:

1 1 1 2

2 2 1 2

1 2

' ( , , ,..., )

' ( , , ,..., )

' ( , , ,..., )

n

n

n n n

y f x y y y

y f x y y y

y f x y y y

Para que possamos encontrar as soluções, deveremos ter as condições iniciais:

1 0 1 0 2 0 2 0 3 0 3 0 0 0( ) ( ) , ( ) ( ) , ( ) ( ) , ... , ( ) ( )n ny x y y x y y x y y x y

Todos os métodos vistos até aqui, podem ser usados para resolver esses sistemas.

Ex1 : Vejamos um caso de um sistema de 2 equações.

0 0

0 0

' ( , , )

' ( , , )

( )

( )

y f x y z

z g x y z

y x y

z x z

10.10 Usando método de Euller:

1

1

. ( , , )

. ( , , )

n n

n n

y y h f x y z

z z h g x y z

10.11 Usando método de Euller Aprimorado

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 61

1 1 2.( )2

n n

hy y k k

1 1 2.( )2

n n

hz z l l onde

1 1

2 1 1 2 1 1

( , , ) e ( , , )

( , . , . ) e ( , . , . )

n n n n n n

n n n n n n

k f x y z l g x y z

k f x h y h k z h l l g x h y h k z h l

10.12 Usando método de Runge-Kutta

1 1 2 3 4

1 1 2 3 4

1 1

2 1 1 2 1 1

3 2 2 3

2 26

2 2 6

( , , ) e ( , , )

( , , ) e ( , , )2 2 2 2 2 2

( , , ) e ( ,2 2 2 2

n n

n n

n n n n n n

n n n n n n

n n n n n

hy y k k k k

hz z l l l l onde

k f x y z l g x y z

h h h h h hk f x y k z l l g x y k z l

h h h hk f x y k z l l g x y

2 2

4 3 3 4 3 3

, )2 2

( , , ) e ( , , )

n

n n n n n n

h hk z l

k f x h y hk z hl l g x h y hk z hl

Ex 1: Considere o sistema de equações diferenciais de ordem um:

2

2

' 2

' ( )

(0) 1

(0) 1

y y yz

z xy y sen z

y

z

Encontre, através do método de Euller Aprimorado, 1 2 y e y com 0.1h .

Solução:

----------------------------------///----------------------------

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 62

CAPITULO 11

EQUAÇÕES DIFERENCIAIS DE ORDEM SUPERIOR

11.1 Equações Diferenciais de Ordem Superior

Consideremos a equação diferencial de ordem n:

( ) ( 1)( , , ',..., )n ny F x y y y com as condições iniciais

( 1) ( 1)

0 0 0 0 0 0 0 0( ) ; '( ) ' ; ''( ) '' ; ... ; ( )n ny x y y x y y x y y x y

Todos os métodos vistos até agora, sevem para encontrar soluções de p.v.i de ordem

um., para resolver os p.v.i de ordem n, vamos reduzi-los a um sistema de n equações

diferencias de ordem um e usar os métodos para resolver sistemas já vistos.

Sem perda de generalidade, vamos tomar um exemplo de ordem 3.

Ex 1:

0 0

0 0

0 0

"' ( , , ', '')

( )

'( ) '

''( ) ''

y F x y y y

y x y

y x y

y x y

Tomando

' '' ' ( , , , ) ' '' ' ( , , , )y z y z g x y z w z w z w u x y z w termos o

sistema:

0 0 0 0 0 0

' ( , , , )

' ( , , , )

' ( , , , )

( ) ; ( ) ; ( )

y f x y z w z

z g x y z w

w u x y z w

y x y z x z w x w

Ex 2: Encontre, através do método de Euller Aprimorado as soluções do p.v.i dado por:

''

(1) 1

'(1) 0

xy y e

y

y

[1; 1.2] 0.1x com h

Solução:...

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 63

------------------------------------///--------------------------------------

Ex 3: Reescreva o p.v.i na forma de um sistema de equações dif. de ordem um. 2 (4)( 1) ''' ' cos( )

(1) 0; '(1) 2; ''(1) 1; '''(1) 2

xx y e y y x

y y y y

Solução:...

----------------------------------///-------------------------------------

Ex 4: Encontre através do método de Euller,com h = 0.15, 1 2 y e y do p.v.i dado por:

23 ''' '' ' 1

(2) 0; '(2) 1; ''(2) 2

y xy y y x

y y y

Solução;...

-----------------------------------------///---------------------------------------

Ex 5: Escreva como ficaria o método de Runge-kutta para encontrar as soluções de um p.v.i

de ordem 4.

______________________________________________Cálculo Numérico-Prof.Sérgio Marcussi Gaspechak 64

CAPITULO 12

EQUAÇÕES DIFERENCIAIS PARCIAIS

(...)