Fisica_Comp(1)

Prévia do material em texto

Física Computacional
Professor Klaus Cozzolino
Universidade Federal do Pará
Instituto de Ciências Exatas e Naturais
Faculdade de Física
3 de junho de 2014
2
Sumário
1 Introdução 5
1.1 Ementa da disciplina . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Objetivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Conceitos iniciais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.4 Fluxogramas e Algoritmos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.5 Comandos e estruturas no ambiente MATLAB ou OCTAVE . . . . . . . . . . . . . . . 8
2 Integração e Diferenciação Numéricas 15
2.1 Integrais definidas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.1 Regra do Trapézio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.1.2 Regra de Simpson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.2 Quadratura Guassiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.1 Quadratura Gauss-Legendre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.3 Integrais Impróprias . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3.1 Quadratura Gauss-Laguerre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.4 Diferenciação Numérica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.4.1 Aproximações de primeira e segunda ordem para f ′(x) . . . . . . . . . . . . . . 23
2.4.2 Aproximações de primeira e segunda ordem para f ′′(x) . . . . . . . . . . . . . . 24
2.4.3 Aproximações de f ′(x), f ′′(x) e f ′′′(x) cúbicas e quárticas . . . . . . . . . . . . 24
2.5 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3 Raízes de Funções de uma Variável 29
3.1 Método da Iteração Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.1.1 Convergência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2 Método de Newton-Raphson . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.2.1 Convergência . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.3 Método da Bissecção . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.4 Método da Secante (Régula Falsi) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.5 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3
4 SUMÁRIO
4 Ajuste de Curvas 35
4.1 Mínimos Quadrados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4.1.1 Exemplos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
4.1.2 Mínimos Quadrados Ponderados . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.2 Ajuste por uma função qualquer (Newton e Gauss-Newton) . . . . . . . . . . . . . . . 40
5 Sistema de Equações Lineares 47
5.1 Definições . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.2 Método de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.3 Método de Gauss-Jordan . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.4 Decomposição LU (ou Fatoração LU) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.5 Método Iterativo de Jacobi-Richardson . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.6 Método de Gauss-Sidel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.7 Convergência dos métodos Iterativos . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.8 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
6 Solução de Equações Diferenciais Ordinárias 57
6.1 Decaimento Radioativo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
6.2 Equação de Poisson 1-D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
6.2.1 Algoritmo de Thomas - Sistemas tridiagonais . . . . . . . . . . . . . . . . . . . 60
6.3 Oscilações Simples, Amortecidas e Forçadas . . . . . . . . . . . . . . . . . . . . . . . . 62
6.3.1 Oscilador Harmônico Simples . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
6.3.2 Oscilador Harmônico Amortecido . . . . . . . . . . . . . . . . . . . . . . . . . . 64
6.3.3 Oscilações Forçadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6.4 Algoritmo Verlet (Velocity-Verlet) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
6.5 Sistemas de Equações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6.6 Método de Runge-Kutta 4ª Ordem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
7 Solução de Equações Diferenciais Parciais 75
7.1 Equação da Difusão 1-D (caso transiente) . . . . . . . . . . . . . . . . . . . . . . . . . 75
7.2 Soluções Numéricas da Equação da Difusão . . . . . . . . . . . . . . . . . . . . . . . . 76
7.2.1 Forward no tempo e centrada no espaço (1-D) . . . . . . . . . . . . . . . . . . . 76
7.2.2 Esquema de Crank-Nicolson (1-D) . . . . . . . . . . . . . . . . . . . . . . . . . 76
7.3 Solução da Equação de Schrödinger 1-D transiente (pacote Gaussiano) . . . . . . . . . 77
7.4 Equação da Onda . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
7.5 Equação de Laplace e Poisson 2-D (caso estacionário) . . . . . . . . . . . . . . . . . . . 82
7.6 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
Referências Bibliográficas 89
A Obtendo Pontos e Pesos para Quadraturas 91
A.1 Quadratura Gaussiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
A.2 Quadratura Gauss-Legendre . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92
Capítulo 1
Introdução
1.1 Ementa da disciplina
Introdução a linguagem científica de programação a ser usada no curso: FORTRAN (preferencial-
mente), MATHEMATICA, MAPLE, MATLAB (opcionais). Introdução aos métodos básicos de cálculo
numérico: zeros de função, ajuste de dados, integração, diferenciação e solução de sistemas algébricos.
Solução de equações diferenciais ordinárias que modelem sistemas físicos por métodos numéricos.
Resolução numérica de equações diferenciais parciais em uma e duas dimensões, Transformada de
Fourier, Processos Estocásticos (opcional), Dinâmica Molecular (opcional).
1.2 Objetivos
Fornecer meios para que o estudante receba, por uma abordagem moderna, as informações básicas,
mas necessárias para o uso das facilidades computacionais existentes na maioria das Faculdades e
Universidades.
1.3 Conceitos iniciais
Linguagem de programação: comunicação pessoa ⇐⇒ computador.
Como exemplo de linguagens de programação temos:
• ASSEMBLY - Esta linguagem de programação foi criada no meio da década de 50 para substituir
a linguagem de máquina (linguagem binária que a Unidade Central de Processamento ou CPU
entende).
• ALGOL - ALGOrithmic Language. Criada em 1958 é constituída por uma mistura de linguagens
de programação. Considerada a primeira linguagem de programação estruturada.
• COBOL - COmmon Business Oriented Language. Criada em 1960 com objetivo de manipular
grandes arquivos. Utilizado por bancos e empresas para manutenção de contas correntes, folhas
de pagamento, controle de estoque e outras aplicações.
• BASIC - Beginner’s All-purpose Symbolic Instruction Code. Foi criada em 1964 no Dartmouth
College e muito utilizada nos primeiros computadores pessoais (PCs).
• VBASIC - Visual Basic. Evolução da Linguagem Basic, muito usada na construçãode aplicativos
para sistemas Windows como MSOffice e outros aplicativos gráficos.
5
6 CAPÍTULO 1. INTRODUÇÃO
• PASCAL - Criada em 1970 por Niklaus Wirth, tinha por objetivo a substituição do BASIC por
ser uma linguagem estruturada e robusta. As liguagens de programação Turbo Pascal e Delphi
são evoluções do PASCAL.
• C - Foi lançada em 1972 para o desenvolvimento de sistemas operacionais UNIX sua evolução
resultou na linguagem C++ muito utilizada para programação científica.
• FORTRAN - FORmula TRANslation. A primeira versão padronizada foi lançada em 1966 tendo
o nome de FORTRAN66. Aprimoramentos do compilador culminaram nas versões FORTRAN
77, 90 e 95.
• JAVATM - muito utilizada atualmente em ambientes de internet, nos celulares, em vizualizações
3D e aplicativos gráficos. Desenvolvida pela companhia Sun Microsystems é uma linguagem de
programação multiplataforma.
Qualquer que seja a linguagem de programação usada precisamos seguir algumas importantes etapas
para que tenhamos sucesso. São elas
1. planeje os projetos antecipadamente - identifique as principais variáveis a serem usadas nos
cálculos, o método ou métodos numéricos a serem empregados e como os resultados devem ser
apresentados;
2. desenvolva-os por estágios - início, corpo e saídas;
3. modularize-os - se possível estruture seu programa em módulos estanques;
4. mantenha-o simples - códigos complexos dificultam a depuração, ou seja, a identificação e remoção
dos erros tipográficos, lógicos e metodológicos;
5. sempre teste cada estágio - a saída de um estágio (dados, números, decisões) é normalmente
usada no próximo módulo. Se um estágio estiver com erros os resultados do próximo também
estarão errados;
6. documente todo o programa - insira comentários para facilitar a compreensão das etapas sempre
que julgar necessário;
7. obtenha e interprete seus resultados
denominadas de REGRAS de OURO da PROGRAMAÇÃO RINO (2008) ou golden rules, na língua
inglesa .
1.4 Fluxogramas e Algoritmos
Todo programa - ou código fonte - é escrito com um editor de texto normalmente muito simples como o
Notepad (do Windows), o Vi (do Linux), o Edit (no DOS) mas é possível editar também em ambientes
de programação próprios à linguagem em uso como o MSDeveloperStudio para Fortran ou C++ ou
Visual Studio C++.
Os programas mais simples assemelham-se à uma redação possuindo início, meio e fim (ou conclusão).
O início de um programa é normalmente denominado de etapa de inicialização e/ou entrada de dados:
contém informações importantes que serão utilizadas tais como constantes, dimensões de matrizes e
de vetores bem como valores iniciais. O corpo do programa contém as instruções principais a serem
passadas à Unidade de Processamento Central (CPU) da máquina para que a mesma realize as tarefas.
É nesta etapa que os cálculos são realizados através de linhas de comando/instruções específicas à cada
linguagem de programação, i. e.
a = 2 + b, atribua ao a o resultado da soma de 2 com b e
c = a ∗ exp(−b ∗ t), atribua ao c o resultado das operações.
1.4. FLUXOGRAMAS E ALGORITMOS 7
Na conclusão do código fonte encontramos instruções como saída de resultados para arquivos, saí-
das gráficas ou simplismente uma menssagem de finalização da execução. Uma simples saída como
“Aplique na Ação X” ou “NÃO Aplique na Ação Y” pode ser o resultado de inúmeras linhas de ins-
truções de um código executável.
O Algoritmo: é um conjunto resumido de instruções (ou camandos) a serem realizadas pelo compu-
tador.
O Fluxograma: apresenta as instruções de um algoritmo na forma gráfica através de figuras geomé-
tricas com as instruções principais em seu interior.
Exemplos de algoritmos, fluxogramas introdução a programação podem ser encontrados em GUSMAN;
VASCONCELLLOS (1985); MEDINA; CRISTINA (2005).
Exemplo 1: desejamos construir uma tabela com 1000 linhas que contenha em suas colunas o tempo,
a velocidade e a posição de um objeto que se desloca com aceleração constante. Para isso usamos
as equações da Física pertinentes à cinemática inserindo, via teclado, alguns dados iniciais como
posição e velocidade em t = to. Abaixo temos um exemplo de fluxograma e um do algoritmo
desenvolvido para tal fim.
.
1. Inicialize as variáveis: N = 0, dt = 0.1, a, vo,
xo, to
2. Calcule:
novo tempo t = to + N*dt;
nova velocidade v = vo + a*t;
nova posição x = xo +vo*t +0.5*a*t^2
3. Escreva em TAB.DAT: t , v e x
4. Incremente o valor de N: N = N + 1
5. Se N < 1000 vá para o passo 2 caso contrário
6
6. Fim
Figura 1.1: Exemplo de Fluxograma (à esquerda) e Algoritmo (à direita) para geração de uma tabela
de dados (arquivo TAB.DAT) com 1000 linhas contendo em cada linha tempo (t), velocidade (v) e
posição (x).
Observe e reflita: no fluxograma encontramos uma figura de decisão (losângulo) com a instrução “
N é maior ou igual a 1000? “. Porém, no item 5 do algoritmo tem-se, “ Se N for menor do que 1000
vá ...”. Estas instruções são claramente diferentes mas realizam o mesmo propósito: a repetição das
instruções necessárias a construção de uma tabela contendo o tempo, a velocidade e a posição de um
móvel em movimento retilíneo uniformemente acelerado - MRUA.
8 CAPÍTULO 1. INTRODUÇÃO
1.5 Comandos e estruturas no ambiente MATLAB ou OCTAVE
Como em qualquer ambiente de programação é importante sabermos passar para o computador um
conjunto de comandos para que o mesmo compreenda e execute. A Tabela 1.1 abaixo contém alguns
comandos e estruturas encontrados na maioria dos programas científicos. Os dois primeiros são es-
truturas de laço usados para executar repetidamente uma instrução ou um conjunto de instruções ou
operações com diferentes dados a cada laço. O terceiro caso é usado em situações que exijam decisão,
por exemplo, se o relógio estiver parado dê corda e acerte a hora, caso contrário não faça nada.
Informações adicionais das funções pré-definidas no ambiente Matlab são encontradas em JALÓN;
RODRÍGUEZ; VIDAL (2001); JALÓN; RODRÍGUEZ; VIDAL (2005) bem como nos arquivos de
ajuda (help) que vêm junto com o aplicativo THE MATHWORKS (2002a,b,c)
Exemplo 2: para criar um vetor linha com quatro elementos escrevemos no “prompt” do ambiente
Matlab (>‌>):
>> v = [ 0 . 0 0 .5 0 .9 1 . 2 ] ;
e teclamos enter. Os elementos do vetor devem estar separados por espaço ou vírgulas. Um vetor
com elementos igualmente espaçados de 0.1, que inicie com 0 e termine com 0.5, pode ser criado
usando a linha de comando:
>> vi = 0 : 0 .1 : 0 . 5 ;
o que resultará em um vetor com seis elementos vi = [0.0 0.1 0.2 0.3 0.4 0.5]. Se desejamos definir
uma matriz 4x4 usamos
>> M = [ 0 . 5 3 .0 2 .5 3 .6 ; 1 . 0 6 .0 3 .0 2 .8 ; 2 . 0 5 . 0 3 .0 4 .0 ; . . .
1 . 0 8 . 0 4 .0 2 . 0 ] ;
onde a quebra da linha da matriz é indicada pela inserção do ’ ; ’. Multiplicar a matriz M por
um vetor coluna vt ou seja o vetor transposto de v fazemos:
>> v1 = M ∗ v ’ ;
Já uma matriz inversa pode ser obtida com o comando inv :
>> IM = inv (M) ;
Para obtermos o quadrado dos valores de um vetor usamos:
>> v2 = v .∗ v ;
cujos elementos serão v2 = [02 0.52 0.92 1.22] = [0 0.25 0.81 1.44], porém se fizermos
>> re s = v ∗ v ’ ;
o resultado será o produto interno de ~v · ~v = 02 + 0.52 + 0.92 + 1.22= 2.5. A sutil diferença - o
ponto antes do sinal de multiplicação - instrui a CPU que: multiplique o primeiro elemento de v
(antes de .* ) com o primeiro elemento de v (depois de .*) assim sucessivamente até os últimos
elementos dos vetores, armazenando cada resultado no vetor v2. Se usarmos
>> v2 = v .^2 ;
também obtemos um vetor resultante v2 idêntico ao anterior onde cada elemento de v2 é o
quadrado do respectivo elemento em v. Se desejamos criar uma matriz 3x3 com todos os elementos
nulos podemos usar o comando zeros na forma
>> MZ = ze ro s ( 3 , 3 ) ;
e atribuímos valores a cada elemento da matriz usando por exemplo
1.5. COMANDOSE ESTRUTURAS NO AMBIENTE MATLAB OU OCTAVE 9
Tabela 1.1: Estruturas e comandos iniciais em MATLAB/OCTAVE.
Comando/ Estrutura Exemplo Comentário
for ... end
for i = 1 : N;
comando 1 ;
comando 2; ...
end;
estrutura de laço - para i
variando de 1 até N repita
comando 1 e comando 2
while ... condição ...
end;
while condição;
comando1 ;
comando2 ; ...
end;
estrutura de laço - enquanto
condição for verdadeira
repita as instruções comando
1 e comando 2
if ... else ... end;
if condição;
comando 1;
else;
comando 2;
end;
estrutura de decisão - se
condição for verdadeira
execute comando 1 caso
contrário exec comando2
+ e - y = z + w - I soma e subtração
* e / x = vo + a * t; y = z / x ; multiplicação e divisão
^ y = yo + a * t^2 potenciação ou radiciação
cos , acos , sin , asin,
tan , atan , sec , asec ,
...
a = cos(x) , arc = asin(M) operações trigonométricas esuas inversas (radianos)
whos whos lista as variáveis e bem comosuas dimensões
clear clear ; clear x y z
suprime todas as variáveis do
espaço de trabalho ativo e
suprime apenas x, y e z
length( ) length(x) ; length(M)
no. de elementos do vetor x e
max. dimensão da matriz M
zeros( ) , ones( ) M=zeros(3, 2) ; N=ones(3, 2)
inicializa uma matriz 3x2 com
todos elementos nulos ou uns
( ... ) y = 2*pi / (3 * y) indica preferência na operação
help comando ou
função help plot
lista uma ajuda sobre o
comando plot
load load nome.ext -ascii
carrega na memória os dados
contidos no arquivo nome.ext
e cria uma matriz no ambiente
Matlab
save save nome.ext variável -ascii
salva no diretório corrente
os dados armazenados em
variável no arquivo
nome.ext
’ vt = v’ ; Mt = M’ ;
transposição do vetor/matriz
(valores reais) ou
transposição do complexo
conjugado
10 CAPÍTULO 1. INTRODUÇÃO
>> MZ(1 ,1 ) = 5 . 0 ; % o elemento de MZ l inha 1 coluna 1 passa a s e r 5 . 0
>> MZ(2 ,3 ) = 1 . 0 ; % o elemento de MZ l inha 2 coluna 3 passa a s e r 1 . 0
Exemplo 3: O código fonte depende da linguagem de programação.
Código fonte do Exemplo 1 em FORTRAN
program tabela
! cria TAB.DAT c/ t v e x
integer :: N, Nlin=1000
real :: dt=0.1
real :: t, v, x, to, vo, xo
open(unit=10, file=’TAB.DAT’,&
status=UNKNOWN)
print*,’entre c/ a,to,vo,xo’
read*,a,to,vo,xo
N = 0
do while (N < Nlin)
t = to + N * dt
v = vo + a * t
x = xo + vo *t + 0.5* a * dt
write(10,*) t, v, x
N = N + 1
end do
close(10)
write(6,*) ’FIM DA EXECUÇÃO’
stop
end
Código fonte do Exemplo 1 em MATLAB/OC-
TAVE
% tabela.m
% cria TAB.DAT c/ t v e x
Nlin = 1000;
dt = 0.1;
a =input(’Entre c/ a :’);
to=input(’Entre c/ to:’);
vo=input(’Entre c/ vo:’);
xo=input(’Entre c/ xo:’);
fid = fopen(’TAB.DAT’,’w’);
N = 0;
while N < Nlin;
t = to + N * dt;
v = vo + a * t;
x = xo + vo *t + 0.5* a * dt;
fprintf(fid,’%e %e %e \n’,t,v,x);
N = N +1;
end;
fclose(fid);
disp(’FIM DA EXECUÇÃO’)
Figura 1.2: Códigos fonte em FORTRAN (à esquerda) e Matlab ou OCTAVE (à direita).
No código para Matlab acima as duas primeiras linhas são comentários e portanto não são
executadas: começam com “ % ”. Na terceira e quarta linhas são inicializadas as variáveis
número de linhas da tabela (Nlin) e incremento do tempo (dt). Da quinta a oitava linha vemos
as instruções de entrada de dados via teclado (comando input). O arquivo “ TAB.DAT “ é
aberto para escrita na nona linha. Na décima linha inicializa-se o contador “ N ” necessário à
estrutura de laço “ while ” compreendido entre as linhas 11 à 18. Dentro do laço as instruções
serão executadas enquanto a condição “ N < Nlin “ for verdadeira. Os valores de t, v e x são
alterados a cada laço e o comando fprintf instrui a máquina a escrever os valores das variáveis
em uma linha no formato científico (%e) 1.2345678e00 que equivale a 1.2345678×1000. O “ \n “
instrui o computador a pular de linha após escrever os valores numéricos de t, v e x no arquivo,
isto para cada valor do contador N . Na sequência o valor de N é incrementado de uma unidade
para um nova repetição das instruções dentro do laço. Quando N for igual a 1000 o laço é
interrompido e o computador irá executar a linha posterior ao “ end “ ou seja irá fechar o arquivo
“ TAB.DAT ” para então escrever na janela de comando do ambiente Matlab o fim da execução
do script “ tabela.m “.
Deve-se ressaltar que em ambos os códigos existe um erro no cálculo de x(t): o último termo da
equação deveria ter 0.5*a*t*t. Este tipo de erro não é incomum, ou o programador desconhece a
equação da cinemática ou ocorreu erro de digitação. Em ambos os casos os resultados da posição
da partícula no arquivo TAB.DAT estarão errados!
Exemplo 4: A série de Taylor da funçao f(x) = senx é dada por
f(x) = senx = x− x
3
3!
+
x5
5!
− x
7
7!
+ · · · =
∞∑
k=1
(−1)k+1x2k−1
(2k − 1)!
1.5. COMANDOS E ESTRUTURAS NO AMBIENTE MATLAB OU OCTAVE 11
Para responder quantos termos são necesários para que o resultado da série em x = pi/3 tenha
precisão de 4 casas decimais foi desenvolvido o seguinte algoritmo
1- limpar memória
2- definir soma como zero
3- definir precisão desejada
4- inicializar o contador k ← 0
5- definir ponto x
6- enquanto |termo| > = precisao faça
k ← k + 1
termo ← (-1)^(k+1) * x^(2*k -1) / (2*k -1) !
soma ← soma + termo
fim do laço 6
7- escreva o valor do contador k
8- escreva o resultado do seno de x
Observações: no ambiente Matlab/Octave o fatorial de um número (num) é obtido usando
a instrução factorial(num). Um erro frequente ao codificar o algoritmo acima para uma lingua-
gem de programação pode ocorrer: o valor da variável termo não foi definido antes do teste na
instrução de número 6. O valor absoluto de um número é obtido com a função abs(num).
Exemplo 5: Um número denominado de precisão da máquina (eps) tem valor aproximado de 2, 22×
10−16. Ao somar ou subtrair este valor de um escalar o resultado permanece inalterado. Com o
intúito de calcular o eps o seguinte algoritmo foi desenvolvido
1- limpar memória
2- definir um como 1
3- definir eps como um
4- enquanto eps + um > um faça
eps ← eps/2
fim do laço 4
5- escreva o valor encontrado de eps
Exemplo 6: Segundo Butkov (1988, p.362) a função de Bessel de primeiro tipo e ordem zero pode
ser expressa segundo
J0(x) =
∞∑
k=0
(−1)k
(k!)2
(x
2
)2k
.
Em virtude de limitações computacionais esta série deverá ser truncada em N termos antes de
ser utilizada em um código (programa). Estamos interessados em determinar quantos termos
devemos usar para atingir o limite de precisão do computador.
Para simplificar a estimativa do número de termos da série cosideremos x = 2. O fator principal
que determinará a magnitude do termo da série será o denominador (k!)2 logo, ao igualar este
termo ao eps teremos
2, 22× 10−16 = 1
(N !)2
→ N ! =
√
1016
2, 22
→ N ! ≈ 68× 106.
Usando divisões sucessivas teremos: 68 × 106/2/3/4/5/6/7/8/9/10/11 = 1, 7 ⇒ N ≈ 11, ou
seja, são necessários aproximadamente onze termos da série para alcançar a precisão da máquina.
Para testar a estimativa acima o seguinte algoritmo foi desenvolvido
1- limpar memória
2- definir soma como zero
3- definir precisão desejada (= eps)
4- inicializar o contador k ← 0
5- definir ponto x← 2
12 CAPÍTULO 1. INTRODUÇÃO
6- definir termo (deve ser maior que o valor em precisão)
7- enquanto |termo| > = precisao faça
termo ← (-1)^k * (x/2)^(2*k) / (k !)^2
soma ← soma + termo
k ← k + 1
fim do laço 7
8- escreva o valor do contador k
Exercícios
1.1- Edite um arquivo de nome “tabela.m” copiando o script (código fonte) para Matlab/Octave
apresentado na Figura1.2. No ambiente Matlab execute o código e verifique a presença do arquivo
TAB.DAT no diretório de trabalho, após a execução do script. Abra o arquivo TAB.DAT com o
editor e verifique se número de linhas está correto.
1.2- Ao final do código (script) “tabela.m” adicionea linha load TAB.DAT e abaixo outra linha
plot(TAB(:,1), TAB(:,3),’b–’); grid; salve, execute o código e descreva o que aconteceu.
Mude ’b–’ no comando gráfico plot para ’ro’, salve, execute o script e descreva o que ocorreu.
Digite help plot no prompt (>‌>) do ambiente MATLAB, leia o texto e tire suas conclusões.
1.3- (a) Crie um algoritmo com intuito de gerar um arquivo de dados (mhs_a.dat) que contenha
100 linhas de modo que cada linha possua o tempo (t), a deformação da mola (xs) do MHS e
a deformação da mola (xa) do movimento harmônico amortecido. (b) Transcreva o algoritmo
para um fluxograma. (Dica: reveja o conteúdo de oscilações para obter as equações horárias de
x(t) para cada sistema físico).
1.4- Reveja o exemplo 5 neste capítulo. Codifique o algoritmo acima em um script para Matlab/Octave
e verifique o valor encontrado do eps.
1.5- Crie um script para o ambiente Matlab/Octave que gere um vetor de elementos igualmente
espaçados de 1 a 1001, com incremento de 10. Em seguida use a estrutura de laço
for i = inicio : incremento : fim;
comandos;
end;
para somar os elementos do vetor. ( Obs.: no Matlab obtemos o valor numérico do elemento
i=2 de um vetor usando v(i). Isto é equivalente a escrever valor=v(2), o que atribui a variável
valor o número armazenado na posição 2 do vetor v. Dica: reveja o script na Figura 1.2 na
página 10 e preste atenção no contador N. Esta estrutura também é chamada de acumulador e
pode ser ligeiramente modificada para acumular as sucessivas somas dos elementos do vetor.)
1.6- O vetor que posiciona o centro de massa (CM) de um sistema discreto de massas é calculado por
~rCM = xCM iˆ+ yCM jˆ + zCM kˆ com xCM =
∑
mixi/
∑
mi
e assim sucessivamente. Considere um sistema de seis massas, m = [1.0, 3.0, 5.0, 2.0, 4.0, 3.0]
no plano xy. Suas posições no eixo x e y são
x = [0.0, 0.5, 1.0, 2.0, 3.0, −5.0] e y = [−1.0, 0.0, 2.0, 4.0, 6.5, 8.0],
respectivamente. Faça um script para Matlab/OCTAVE que a partir dos dados de m, x e y
calcule a posição do centro de massa ~rCM . (Use estrutura de laço, um somador para xCM e
outro somador para yCM ). Adicione a instrução plot(x, y, ’ro’, xcm, ycm, ’bs’) e execute
novamente o script. Descreva o resultado gráfico.
1.5. COMANDOS E ESTRUTURAS NO AMBIENTE MATLAB OU OCTAVE 13
1.7- Considere n ∈ [1, 1000]. Escreva um código usando instruções na linguagem Matlab/Octave que
obtenha o resultado da soma dos números pares bem como dos números ímpares entre 1 e 1000.
Para tanto use uma estrutura de decisão
Se condição for verdadeira faça
somapar ← somapar + n
caso contrário faça
somaimp ← somaimp + n
fim do Se
que teste se o número é par (condição) usando o resto da divisão entre dois números mod(x, y)
ou rem(x, y). Veja na ajuda a documentação sobre a instrução mod: use help mod.
1.8- Uma carga teste (próton) encontra-se no vácuo sob ação
de um campo elétrico uniforme de 10−5V/m e sentido
como indicado na Figura 1.3. Em t = 0s sua posição é
de 1 m a direita da origem e possui velocidade de 53,6
m/s da direita para esquerda. Assumindo a massa de
1, 67× 10−27Kg e carga de 1, 6× 10−19C
(a) escreva a equação horária de sua posição e
(b) faça um script (Matlab/Octave) para encontrar
os dois instantes em que o próton cruza a origem do sis-
tema de coordenadas: defina as quantidades pertinentes;
com elas obtenha a aceleração da carga teste; (o script
deve resolver a equação do segundo grau em t);
(c) crie um vetor com valores de tempo no intervalo
[0, 0,15]s e incremento de 10−4s; em seguida usando a
equação horária x(t) do item (a) calcule as posições do
próton nestes tempos e faça um gráfico usando a instru-
ção plot(t, x);
(d) a partir do gráfico estime a máxima distância ne-
gativa alcançada pelo próton e
(e) em que instante ele estará na posição x = 2m.
Figura 1.3: Próton em campo elétrico
uniforme. Esquema do exercício 1.8.
Observação: neste exercício deconsideramos as pertubações geradas pelo campo da carga teste
sobre o campo elétrico uniforme bem como o campo irradiado por ela em virtude de sua pequena
aceleração, como preconiza o eletromagnetismo clássico.
1.9- Considere os dados do exercício anterior e a expressão da massa relativística M dada por
M =
m√
1−
(v
c
)2
onde v é a velocidade instantânea do próton e c é a velocidade da luz.
(a) Escreva um código para gerar um gráfico que ilustre a dependência da massa do próton
com o tempo (o gráfico é um dos produtos finais) ;
(b) execute o código considerando o intervalo de tempo do exercício anterior e
(c) altere a linha do código que gera o tempo: redefina o tempo final para 1000 s e incremento
de 1 s; execute-o e analizando o gráfico responda
(d) é necessário considerar o efeito relativístico da massa no problema ?
(e) Quais são os valores de x e v em t = 1000 s ? (Use x(1000) e v(1000) seguidos de <enter>
no prompt da janela de comando).
14 CAPÍTULO 1. INTRODUÇÃO
Capítulo 2
Integração e Diferenciação Numéricas
Cálculo de integrais em Física são frequentes: as Leis de Maxwell na forma integral, o valor do trabalho
de uma força que depende da posição ou o impulso provocado por determinada força que depende do
tempo são alguns exemplos. As vezes a função a ser integrada leva a manipulações algébricas complexas
e o importante em determinadas situações é o resultado numérico para interpretação.
Este capítulo têm por finalidade ilustrar alguns aspectos básicos dos processos de integração e diferen-
ciação numérica aplicadas a funções reais de uma variável.
2.1 Integrais definidas
2.1.1 Regra do Trapézio
Consideremos uma função f(x) conhecida, contínua e definida para qualquer valor de x ∈ [a, b].
Desejamos calcular o resultado da integral definida
I =
ˆ b
a
f(x) dx,
sendo f(x) finita no intervalo considerado. Uma possibilidade seria usar a definição da Integral de
Riemann (SPIEGEL, 1975, p.102)
I =
ˆ b
a
f(x) dx = lim
∆x→0
∑
i
f(ci)∆xi , (2.1)
onde ci é um valor compreendido no intervalo [xi, xi + ∆xi]. Fazendo ∆x pequeno, desconsideramos
o limite e assim obtemos a Regra dos Retângulos.
15
16 CAPÍTULO 2. INTEGRAÇÃO E DIFERENCIAÇÃO NUMÉRICAS
Outro método mais preciso é a Regra dos Trapézios:
discretizamos o intervalo [a, b] em N − 1 segmentos
de largura ∆x onde N é o número de avaliações da
função. Dado um conjunto de valores de xi a saber
x = [a, (a + ∆x), (a + 2∆x) + . . . . . . + (a + (N −
1)∆x), b] = [x1, x2, . . . , xN ] a área genérica Areai de
um trapézio é dada por (ver Figura (2.1))
Areai =
1
2
[f(xi) + f(xi+1)] ∆x
e, somando para todos os intervalos de 1 a N −1 vêm
I =
[
1
2
f(x1) + f(x2) + . . .+ f(xN−1) +
1
2
f(xN )
]
∆x.
Figura 2.1: Elemento genérico de área do mé-
todo dos trapézios.
Como f(x1) = f(a) e f(xN ) = f(b) podemos escrever
Iab =
{
N∑
i=1
f(xi)− 1
2
[f(a) + f(b)]
}
∆x =
{
N∑
i=1
fi − 1
2
[f(a) + f(b)]
}
∆x. (2.2)
Exemplo 1: um ’script ’ simples para implementar o método acima encontra-se a seguir.
% Reg_Trap.m
% use para dados com ’dx = cte.’
clear;
N1=input(’Entre com o no. de intervalos (N-1): ’);
N =N1+1;
a =input(’Entre com o limite inferior a: ’);
b =input(’Entre com o limite superior b: ’);
dx=(b-a)/N1;
soma=0;
for i = 1:N ; % inicio do laço para soma
x(i)= a+(i-1)*dx;
f(i)= função( x(i) ) % expressão da função a ser integrada nesta linha
soma = soma + f(i);
end; % fim do laço
disp(’O resultado da integral eh:’)
I = dx*( soma - 0.5*( f(1) + f(N) ) )
2.1.2 Regra de Simpson
Uma forma alternativa e normalmente mais precisa do que a Regra dos Trapézios é a Regra de
Simpson, onde se faz uso de três avaliações da função para definir a área elementar. Neste caso
podemos escrever de forma compacta a expressão desta regra a saber
Iab =
∆x
3
(N−1)/2∑
i=1
[f(x2i−1) + 4f(x2i) + f(x2i+1)]
Iab =
∆x
3
(N−1)/2∑
i=1[f2i−1 + 4f2i + f2i+1] . (2.3)
2.1. INTEGRAIS DEFINIDAS 17
Devemos ter em mente que neste método o número de intervalos (N − 1 ) entre a e b deve ser par.
Ao deduzir sua fórmula Simpson aproxima a função por uma parábola em cada conjunto de 3 pontos
como ilustrado na Fig. 2.2 sendo xi = a + i∆x, xi−1 = a + (i − 1)∆x e xi+1 = a + (i + 1)∆x para
i = 1, 1, 2,...,N-1.
x x xi-1 i i+1
y
x
yi
yi+1
yi-1
y=f(x)
Figura 2.2: Regra de Simpson 3 pontos: aproximação da área elementar abaixo da parábola.
A área elementar é obtida integrando uma função parabólica y(x) = cox2 + c1x+ c2 que se ajusta aos
três pontos de avaliação da função. O valor aproximado da integral será
Ai =
ˆ xi+1
xi−1
(
cox
2 + c1x+ c2
)
dx
=
[co
3
x3 +
c1
2
x2 + c2x
]xi+1
xi−1
=
[co
3
(
x3i+1 − x3i−1
)
+
c1
2
(
x2i+1 − x2i−1
)
+ c2 (xi+1 − xi−1)
]
Ai =
(xi+1 − xi−1)
6
[
2co
(
x2i+1 + xi+1xi−1 + x
2
i−1
)
+ 3c1 (xi+1 + xi−1) + 6c2
]
escrevendo em termos da função y(x) temos
Ai =
(xi+1 − xi−1)
6

y(xi+1)︷ ︸︸ ︷
cox
2
i+1 + c1xi+1 + c2 +
y(xi−1)︷ ︸︸ ︷
cox
2
i−1 + c1xi−1 + c2 +
+co
(
x2i+1 + 2xi+1xi−1 + x
2
i−1
)
+ 2c1 (xi+1 + xi−1) + 4c2
]
=
(xi+1 − xi−1)
6
{
y(xi−1) + y(xi+1) + 4
[
c0
(
xi+1 + xi−1
2
)2
+ c1
(
xi+1 + xi−1
2
)
+ c2
]}
Ai =
(xi+1 − xi−1)
6
[
y(xi−1) + y(xi+1) + 4y
(
xi+1 + xi−1
2
)]
.
Usando as igualdades xi+1 − xi−1 = 2∆x e (xi+1 + xi−1)/2 = xi escrevemos a área elementar abaixo
da parábola na forma
Ai =
∆x
3
[y(xi−1) + 4y (xi) + y(xi+1)] ,
o que leva à Regra de Simpson (Eq. 2.3) ao somarmos para todo o intervalo de integração.
18 CAPÍTULO 2. INTEGRAÇÃO E DIFERENCIAÇÃO NUMÉRICAS
2.2 Quadratura Guassiana
No último exemplo ficou claro a necessidade de se resolver uma integral definida no intervalo [a, b]. A
Quadratura Gaussiana é um método de integração que basicamente substitui a integral por uma soma
do produto de uma função peso w(xi) pela função f(x) a ser integrada, i.e.
ˆ b
a
f(x)dx ≈
n∑
i=1
w(xi)f(xi) (2.4)
onde os valores de xi, pontos de avaliação da função, são definidos em termos do número de pontos
n da soma. Como regra básica se a função f(x) é um polinômio de grau 2n − 1 serão necessários n
pares (xi, wi) para que o resultado da integral seja exato. Para a aplicação deste método é necessário
o conhecimento da função f(x).
Exemplo 2: o ’script ’ a seguir foi adaptado de uma subrotina em Fortran (PRESS et al., 1989). Os
resultados serão precisos caso f(x) seja um polinômio de grau menor ou igual à 19.
% Integ_gauss.m
% realiza a integral de f(x) por quadratura gaussiana
% preciso para polinômios de grau <= 19
clear;
% define pontos da quadratura xi
xi = [.1488743389,.4333953941,.6794095682,.8650633666,.9739065285];
%
% define os pesos wi
wi = [.2955242247,.2692667193,.2190863625,.1494513491,.0666713443];
%
% entrada dos limites a e b
a = input(’Entre c/ limite inferior: ’);
b = input(’Entre c/ limite superior: ’);
%
% define a função a ser integrada
ff = inline(’cos(8*x).*exp(-x)’,’x’)
%
%
xm = 0.5*(b+a);
xr = 0.5*(b-a);
soma = 0;
for j = 1:5;
dx = xr*xi(j);
soma = soma + wi(j)*( ff(xm+dx) + ff(xm-dx) );
end;
disp(’O resultado da integral eh: ’)
soma = soma*xr
Exemplo 3: o código acima pode ser reescrito usando o formalismo matricial a saber:
% Integ_gauss1.m
% realiza a integral de f(x) por quadratura gaussiana
clear;
% define pontos da quadratura xi
2.2. QUADRATURA GUASSIANA 19
xi = [.1488743389,.4333953941,.6794095682,.8650633666,.9739065285];
%
% define os pesos wi
wi = [.2955242247,.2692667193,.2190863625,.1494513491,.0666713443];
%
% entrada dos limites a e b
a = input(’Entre c/ limite inferior: ’);
b = input(’Entre c/ limite superior: ’);
%
% define a funcao a ser integrada
ff = inline(’cos(8*x).*exp(-x)’,’x’)
%
%
xm = 0.5*(b+a);
xr = 0.5*(b-a);
dx = xr*xi;
soma = wi*( ff(xm+dx) + ff(xm-dx) )’;
disp(’O resultado da integral eh: ’)
soma = soma*xr % desnormaliza o resultado
Exemplo 4: existe no Matlab/OCTAVE algumas funções escritas para este fim. Ao final do ’script ’
abaixo existem duas linhas de comando com as funções predefinidas no ambiente Matlab: quad e
quadl. A forma mais básica de usá-las para o cálculo das integrais é o mesmo: quad(nome_função,
lim_inf, lim_sup, tolerância).
% Integ_gauss2.m
% realiza a integral de f(x) por quadratura gaussiana
clear;
%
% entrada dos limites a e b
a = input(’Entre c/ limite inferior: ’);
b = input(’Entre c/ limite superior: ’);
%
% define a funcao a ser integrada
ff = inline(’cos(8*x).*exp(-x)’,’x’)
%
% realiza a integral - quadratura adaptativa de Simpson
disp(’O resultado da integral eh: ’)
soma = quad(ff,a,b)
%
% realiza a integral - quadratura adaptativa Gauss-Lobato
disp(’O resultado da integral eh: ’)
soma = quadl(ff,a,b)
2.2.1 Quadratura Gauss-Legendre
Quando os limites de integração de uma função são finitos é comum a aplicação de uma substituição
de variáveis
x = xR z + xM , onde
xR =
b− a
2
,
xM =
b+ a
2
.
20 CAPÍTULO 2. INTEGRAÇÃO E DIFERENCIAÇÃO NUMÉRICAS
com o intúito de normalizar os mesmos para o intrervalo [-1, 1]. Assim a integral definida no intervalo
x ∈ [a, b] é reescrita na forma de um somatório a saber
I =
ˆ b
a
f(x) dx =
ˆ 1
−1
f
(
b− a
2
z +
b+ a
2
)[
b− a
2
dz
]
=
b− a
2
ˆ 1
−1
f
(
b− a
2
z +
b+ a
2
)
dz
I ≈ b− a
2
n∑
i=1
Wi f
(
b− a
2
zi +
b+ a
2
)
(2.5)
I ≈ xR
n∑
i=1
Wi f (xR zi + xM ) .
onde zi são as raízes do Polinômio de Legendre Pn(z) com pesos calculados por
Wi =
2(
1− z2i
)
[P ′n(zi)]
2 .
Os polinômios de Legendre são obtidos recursivamente pela equação (ARFKEN; WEBER, 2001, p.749)
Pn+1(z) = 2zPn(z)− Pn−1(z)− zPn(z)− Pn−1(z)
(n+ 1)
.
Para obtê-los iniciamos com n = 0, sabendo que P0(z) = 1 e, recursivamente, derivamos os polinômios
de graus mais elevados. Como exemplo os quatro primeiros polinômios de Legendre são
P1(z) = z,
P2(z) =
1
2
(
3z2 − 1) , P3(z) = 18
(
5z3 − 3z) ,
P4(z) =
1
8
(
35z4 − 30z2 + 3) .
Abaixo um script em Matlab/OCTAVE (TREFETHEN, 2008) para o cálculo dos pontos e pesos
necessários a aplicação do método
% esta funcao deve estar salva no diretório de trabalho
% com o nome GaussLegendre.m
function [x, w] = GaussLegendre(n)
% x - pontos da quadratura (saida)
% w - pesos da quadratura (saida)
% n - numero de pontos da quadratura (entrada)
i = 1:n-1;
a = i./sqrt(4*i.^2-1);
CM = diag(a,1) + diag(a,-1);
[V L] = eig(CM);
[x ind] = sort(diag(L));
V = V(:,ind)’;
w = 2 * V(:,1).^2;
Com esta função salva no diretório de trabalho podemos obter os pontos e respectivos pesos da qua-
dratura Gauss-Legendre usando uma instrução como
>‌> [x , w] = GaussLegendre(3)
x =
-0.77459666924148 0 0.77459666924148
w =
0.55555555555556 0.88888888888889 0.55555555555556
o que atribui a x as raízes do polinômio de Legendre de terceiro grau (pontos da quadratura) e à w os
respectivos pesos. Neste caso é possível integrar com precisão um polinômio, f(x), de grau menor ou
igual à 2n− 1 = 5.
2.3. INTEGRAIS IMPRÓPRIAS 21
2.3 Integrais Impróprias
Algumas integrais possuem o limite inferior, o superior ou ambos igual a ±∞. Neste caso as integrais
são denominadas de impróprias existindo algumas técnicas capazes de simplifica-las numericamente
por substituição de variáveis. Vejamos o exemplo a seguir.
Exemplo 5: calcular numericamente a integral I =
´∞
0
dx
1+x2
. Uma tabela de integrais indefinidas
fornece o resultado como sendo I = tan−1(x) e, substituindo os limites, obtemos I = pi/2.
Agébricamente podemos realizar a seguinte substituição de variáveis: x = 1/y =⇒ dx = −dy/y2.
Ao fazê-lo temos
I =
ˆ ∞
0
dx
1+ x2
=
ˆ 1
0
dx
1 + x2
+
ˆ 0
1
− dy
y2 (1 + 1/y2)
,
=
ˆ 1
0
dx
1 + x2
−
ˆ 0
1
dy
1 + y2
,
I = 2
ˆ 1
0
dx
1 + x2
.
Neste exemplo vimos que a integral no intervalo [0, ∞] foi substituída por outra integral, agora
definida no intervalo [0, 1], sendo facilmente resolvida por um dos métodos descritos nas seções
anteriores. A Tabela 2.1 abaixo ilustra alguns casos de integrais impróprias e suas respectivas
substituições.
Tabela 2.1: Exemplos de integrais impróprias e suas respectivas substituições (PRESS et al., 1989,
cap. 4).
´ b
a f(x)dx
´ 1/a
1/b
1
t2
f(1t )dt , ab>0 Usada quando b→∞ e a > 0
ou quando a→ −∞ e b < 0´ b
a f(x)dx
1
1−γ
´ (b−a)1−γ
0 t
γ
1−γ f(t
1
1−γ + a)dt , b>a No caso do integrando
divergir em (x− a)−γ
próximo de x = a para
0 ≤ γ < 1´ b
a f(x)dx
1
1−γ
´ (b−a)1−γ
0 t
γ
1−γ f(b− t 11−γ )dt , b>a No caso do integrando
divergir em (x− b)−γ próximo
de x = b para 0 ≤ γ < 1´ b
a f(x)dx
´ √b−a
0 2tf(a+ t
2)dt , b>a Funções do tipo 1/[x− a]1/2
com singularidade em x = a´ b
a f(x)dx
´ √b−a
0 2tf(b− t2)dt , b>a Funções do tipo 1/[x− b]1/2
com singularidade em x = b´∞
a f(x)dx
´ exp(−a)
0
1
t f(− log t)dt Integrandos que decaem
exponencialmente exp(−x)
2.3.1 Quadratura Gauss-Laguerre
Similar ao método de Gauss-Legendre porém aplicado à integrais com um dos limites finito e outro
infinito, i.e.
I =
ˆ ∞
0
f(x) dx =
ˆ ∞
0
xα exp(−x)g(x) dx =
ˆ ∞
0
xα exp(−x)
[
exp(x)
xα
f(x)
]
dx
I ≈
n∑
i=1
Wi
[
exp(xi)
xαi
f(xi)
]
(2.6)
22 CAPÍTULO 2. INTEGRAÇÃO E DIFERENCIAÇÃO NUMÉRICAS
onde α > −1, xi são as raízes do Polinômio de Laguerre Ln(z) com pesos Wi calculados por
Wi =
Γ(n+ α+ 1) zi
n!
[
(n+ 1)Lαn+1(zi)
]2 .
com
Lαn(z) =
n∑
m=0
(n+ α)!
(n−m)!(α+m)!
(−z)m
m!
.
Usando a propriedades básicas da funçãoΓ, Γ(n+ 1) = n! obtemos para α = 0
Wi =
zi
(n+ 1)2 [Ln+1(zi)]
2 .
Os polinômios de Laguerre são recurssivamente obtidos pela equação (ARFKEN; WEBER, 2001, p.831)
Ln+1(z) = 2Ln(z)− Ln−1(z)− (1 + z)Ln(z)− Ln−1(z)
(n+ 1)
.
Os cálculos iniciam-se com n = 0 sabendo que L0(z) = 1 fornecendo os seguintes resultados para os
próximos quatro polinômios
L1(z) = −z + 1,
L2(z) =
1
2
(
z2 − 4z + 2) , L3(z) = 16
(−z3 + 9z2 − 18z + 6) ,
L4(z) =
1
24
(
z4 − 16z3 + 72z2 − 96z + 24) .
Abaixo encontra-se um script em Matlab/Octave (WILSON; TURCOTTE, 1998; RECKTENWALD,
2000) para o cálculo dos pontos e pesos da quadratura necessários a aplicação da Quadratura Gauss-
Laguerre.
% esta funcao deve estar salva no diretório de trabalho
% com o nome GaussLaguerre.m
function [x, w] = GaussLaguerre(n, alpha)
% x - pontos da quadratura (saida)
% w - pesos da quadratura (saida)
% n - número de pontos da quadratura (entrada)
% α - potência de x no integrando (entrada)
i = 1:n;
a = (2*i-1) + alpha;
b = sqrt( i(1:n-1) .* ((1:n-1) + alpha) );
CM = diag(a) + diag(b,1) + diag(b,-1);
[V L] = eig(CM);
[x ind] = sort(diag(L));
V = V(:,ind)’;
w = gamma(alpha+1) .* V(:,1).^2;
Observação: o maior real expresso em um sistema operacional 32 ou 64 bits é aproximadamente
1.79× 10308. Usando a instrução >‌> [x, w] = GaussLaguerre(185, 0), ou seja, 185 pontos
e α = 0, obtemos o maior ponto de avaliação da função x185 ≈ 708. Ao inserir este valor na
exponencial temos exp(708) ≈ 3.02× 10307 resultado este muito próximo do limite de expressão
de um número real. Concluí-se assim que a quadratura Gauss-Laguerre deve ser usada com
cuidado pois o número de pontos (e pesos) escolhido pode incorrer no resultado “Inf ”: o que
significa um real acima do valor 1.79× 10308.
2.4. DIFERENCIAÇÃO NUMÉRICA 23
2.4 Diferenciação Numérica
Em alguns ramos da Física, como no Eletromagnetismo, encontramos relações vetoriais entre o campo
elétrico e o campo magnético, com fontes e sumidouros. Um exemplo é o campo elétrico produzido
por uma densidade de carga ρ no vácuo. Assim entendemos da teoria eletromagnética que o campo
elétrico relaciona-se com a densidade de carga elétrica segundo a equação (lei de Gauss)
∇ · (�o−→E ) = ρ
onde �o é a permitividade elétrica do meio (vácuo) e ∇ é o operador Nabla a saber
∇ = ∂∂x iˆ+ ∂∂y jˆ + ∂∂z kˆ coord. cartesianas
∇ = ∂∂r rˆ + 1r ∂∂φ φˆ+ ∂∂z kˆ coord. cilíndricas
∇ = ∂∂r rˆ + 1r sin θ ∂∂φ φˆ+ 1r ∂∂θ θˆ coord. esféricas.
O operador gradiente também está presente na teoria da condução do calor. A Lei de Fourier estabelece
uma relação entre o fluxo de calor e o produto da condutividade térmica do material (k) pelo gradiente
de temperatura em certo ponto do meio, i.e.
−→
H = −k∇T
e que expressa a quantidade de calor por unidade de tempo por unidade de área sendo propagada
através do sólido.
2.4.1 Aproximações de primeira e segunda ordem para f ′(x)
Consideremos o gráfico ao lado de f(x). Desejamos
obter uma aproximação numérica da derivada pri-
meira da função f(x) no ponto x0. Considerando o
triângulo 1 podemos escrever a seguinte equação
f ′(xi) = tan θ =
∆y
∆x
≈ f(xi)− f(xi −∆x)
xi − (xi −∆x)
f ′(xi) ≈ f(xi)− f(xi −∆x)
∆x
sendo conhecida como aproximação atrasada de pri-
meira ordem da derivada primeira (aproximação
backward). De modo análogo, mas usando o tri-
ângulo 2 temos
f ′(xi) ≈ f(xi + ∆x)− f(xi)
∆x
sendo esta equação a derivada numérica de primeira
ordem (erro ≈ ∆x) conhecida por aproximação adi-
antada (forward) para f ′(x) em x = xi.
Figura 2.3: Representações gráficas da derivada
primeira: atrasada, adiantada e centrada.
Outra aproximação da derivada primeira de uma função é a aproximação centrada que pode ser obtida
pela média entre as aproximações backward e forward, ou ainda através do triângulo 3 na Fig. 2.3
resultando em
f ′(xi) ≈ f(xi + ∆x)− f(xi −∆x)
2∆x
sendo esta uma aproximação de segunda ordem (erro ≈ ∆x2), ou seja, mais precisa.
24 CAPÍTULO 2. INTEGRAÇÃO E DIFERENCIAÇÃO NUMÉRICAS
Observação: as derivadas numéricas também podem ser deduzidas através da Série de Taylor da
função f(x), i.e.
f(x+ ∆x) = f(x) + f ′(x) [x+ ∆x− x] + 1
2!
f ′′(x) (x+ ∆x− x)2 + · · ·
que pode ser simplificado para
f(x+ ∆x) ≈ f(x) + f ′(x)∆x+ 1
2
f ′′(x)∆x2 + · · ·
e, retendo apenas os dois primeiros termos da série e manipulanto temos
f ′(x) ≈ f(x+ ∆x)− f(x)
∆x
que representa a aproximação numérica de primeira ordem forward para f ′(x). Se os valores de
x forem discretos como em uma lista ou tabela, i.e., xi = a + (i − 1)∆x para i = 1, 2, 3, . . . por
substituição na equação anterior obtemos a forma compacta a saber
f ′(xi) ≈ f (a+ (i− 1)∆x+ ∆x)− f (a+ (i− 1)∆x)
∆x
,
f ′(xi) ≈ f (a+ i∆x)− f (a+ (i− 1)∆x)
∆x
,
f ′i ≈
fi+1 − fi
∆x
,
onde fi e fi+1 representam os valores da função f(x) obtidos nos pontos xi+1e xi+1, respectiva-
mente.
2.4.2 Aproximações de primeira e segunda ordem para f ′′(x)
A aproximação numérica da segunda derivada de f(x) pode ser obtida facilmente através da formulação
vista na seção anterior porém aplicada a derivada primeira, i.e.
f ′′(xi) ≈ f
′(xi)forward − f ′(xi)backward
∆x
≈
[
f(xi+∆x)−f(xi)
∆x
]
−
[
f(xi)−f(xi−∆x)
∆x
]
∆x
f ′′(xi) ≈ f(xi + ∆x)− 2f(xi) + f(xi −∆x)
∆x2
ou f ′i =
fi+1 − 2fi + fi−1
∆x2
sendo esta a aproximação centrada de segunda ordem (erro ≈ ∆x2) para a derivada segunda de f(x).
Outras formulações numéricas para f ′′(x) são
f ′′(xi) ≈ f(xi+2∆x)−2f(xi+∆x)+f(xi)∆x2 , aproximação (forward) f ′i =
fi+2 − 2fi+1 + fi
∆x2
f ′′(xi) ≈ f(xi)−2f(xi−∆x)+f(xi−2∆x)∆x2 , aproximação (backward) f ′i =
fi − 2fi−1 + fi−2
∆x2
ambas de primeira ordem (erro ≈ ∆x), ou seja, menos precisas.
2.4.3 Aproximações de f ′(x), f ′′(x) e f ′′′(x) cúbicas e quárticas
Alguns cálculos podem exigir aproximações numéricas das derivadas com uma precisão maior. A tabela
abaixo resumealgumas formulações de terceira e quarta ordem para as derivadas numéricas.
2.4. DIFERENCIAÇÃO NUMÉRICA 25
Tabela 2.3: Derivada Primeira f ′ em x0: expressões numéricas, erros e interpolação
Expressão Erro Polinômio Tipo
[−11f(x0) + 18f(x1)− 9f(x2) + 2f(x3)]/6∆x (∆x)3 cúbico Forward
[−2f(x−3) + 9f(x−2)− 18f(x−1) + 11f(x0)]/6∆x (∆x)3 cúbico Backward
[−25f(x0) + 48f(x1)− 36f(x2) + 16f(x3)− 3f(x4)]/12∆x (∆x)4 quarta Forward
[3f(x−4)− 16f(x−3) + 36f(x−2)− 48f(x−1) + 25f(x0)]/12∆x (∆x)4 quarta Backward
[f(x−2)− 8f(x−1) + 8f(x1)− f(x2)]/12∆x (∆x)4 quarta Central
Tabela 2.4: Derivada Segunda f ′′ em x0: expressões numéricas, erros e interpolação
Expressão Erro Polinômio Tipo
[2f(x0)− 5f(x1) + 4f(x2)− f(x3)]/∆x2 (∆x)2 cúbico Forward
[−f(x−3) + 4f(x−2)− 5f(x−1) + 2f(x0)]/∆x2 (∆x)2 cúbico Backward
[35f(x0)− 104f(x1) + 114f(x2)− 56f(x3) + 11f(x4)]/12∆x2 (∆x)4 quarta Forward
[11f(x−4)− 56f(x−3) + 114f(x−2)− 104f(x−1) + 35f(x0)]/12∆x2 (∆x)3 quarta Backward
[−f(x−2) + 16f(x−1)− 30f(x0) + 16f(x1)− f(x2)]/12∆x2 (∆x)3 quarta Central
Tabela 2.5: Derivada Terceira f ′′′ em x0: expressões numéricas, erros e interpolação
Expressão Erro Polinômio Tipo
f ′′′(x0)
[−f(x0) + 3f(x1)− 3f(x2) + f(x3)]/∆x3
[−f(x−1) + 3f(x0)− 3f(x1) + f(x2)]/∆x3 (∆x) cúbico Forward
f ′′′(x0)
[−f(x−3) + 3f(x−2)− 3f(x−1) + f(x0)]/∆x3
[−f(x−2) + 3f(x−1)− 3f(x0) + f(x1)]/∆x3 (∆x) cúbico Backward
f ′′′(x0) [−f(x−2) + 2f(x−1)− 2f(x1) + f(x2)]/2∆x3 (∆x)2 quarta Central
Deve ser observado nas tabelas acima a seguinte notação f−i = f(x−i) = f(x0 − i∆x) e fi = f(xi) =
f(x0 + i∆x) para efeito de cálculos. Também devemos ter em mente ao realizarmos aproximações
numéricas das derivadas a dificuldade em calculá-las nas extremidades dos intervalos. Se temos um
conjunto de dados do tipo (xj , fj) onde j = 1, 2, ..., 11 ou seja, possuímos 11 pares de dados, a
aproximação de segunda ordem da derivada primeira conterá apenas 9 pares de dados pois necessitamos,
para o cálculo da derivada na segunda amostra (j = 2), de (x1, f1) e de (x3, f3). Analogamente, para
derivada em j = 10 são necessários os pares (x9, f9) e de (x11, f11). Uma metodologia muito empregada
nestes extremos dos dados (j = 1 e j = 11) é a utilização das aproximações de primeira ordem, forward
e backward, nos extremos inicial e final respectivamente.
Exemplo 6: Um sistema amortecido oscila segundo a função g(t) = 10sen(0, 5t) exp(−0, 2t). A deri-
vada primeira (analítica) é g′(t) = 5 cos(0, 5t) exp(−0, 2t)−2sen(0, 5t) exp(−0, 2t). Considerando
o intervalo para t de 0 a 10 e um incremento dt = 0.5, determine:
1. A derivada primeira (numérica) da função;
2. O erro entre a aproximação numérica e a expressão analítica de g′(x)
Neste exemplo criamos o seguinte script:
% deriv1.m
clear
dt= 0.5;
t = 0:dt:10;
N = length(t);
fg= inline(’10*sin(0.5*t).*exp(-0.2*t)’);
26 CAPÍTULO 2. INTEGRAÇÃO E DIFERENCIAÇÃO NUMÉRICAS
gl_an = 5*cos(0.5*t).*exp(-0.2*t)-2*sin(0.5*t).*exp(-0.2*t);
%
% usando backward e forward para os ptos. inicial e final
gl_nu(1) = ( fg(t(2)) - fg(t(1)) )/dt; % avançada
gl_nu(N) = ( fg(t(N)) - fg(t(N-1)) )/dt; % atrasada
%
% usando a formulação centrada para pontos intermediários
for i = 2:N-1;
gl_nu(i)= 0.5*( fg(t(i+1)) - fg(t(i-1)) )/dt;
err(i)=gl_an(i)-gl_nu(i);
end;
err(1) = gl_an(1)-gl_nu(1);
err(N) = gl_an(N)-gl_nu(N);
% resultados gráficos
plot(t,gl_nu,’b.’ , t,gl_an); pause
plot(err)
2.5 Exercícios
2.1- Admita a função f(x) = cos(8x) exp(−x2). Realize a integral de f(x) para diferentes discretiza-
ções ∆x (número de intervalos) e verifique para que ∆x o resultado da integral não difere mais
usando
a) regra dos Trapézios
b) regra de Simpson
(Sugestão: escolha um intervalo [0, pi/2] e rode o código para (N-1) = 30, 40, 50, ... ; o
valor de ∆x é obtido escrevendo >‌> dx <enter> no prompt de comando.)
2.2- Para obtermos a posição do C.M. de uma distribuição contínua de massa usamos ~rCM =
(´
~r dm
)
/M
em que xCM =
(´
x dm
)
/M , e assim sucessivamente com M sendo a massa do corpo. Quando
a distribuição de massa não é uniforme a massa específica depende da posição, ou seja, ρ =
ρ(x, y, z). Considere uma vareta de comprimento L = 0, 9m, área da seção reta circular a =
6, 0mm2 e com massa específica dada por ρ(x) = 7+2x−x2 Kg/m3. Use a Regra de Simpson
para estimar o valor numérico de xCM .
2.3- Considere o problema do Exercício 2.2 uma vareta de comprimento L = 0, 9m, área da seção
reta circular a = 6, 0mm2 e com massa específica dada por ρ(x) = 7 + 2x − x2 Kg/m3. Use a
Quadratura Gaussiana (via exemplo 3) para estimar o valor numérico de xCM .
2.4- Considere novamente o problema do Exercício 2.2. Uma vareta de comprimento L = 0, 9m, área
da seção reta circular a = 6, 0mm2 e com massa específica dada por ρ(x) = 7 + 2x− x2 Kg/m3.
Use Quadratura Gauss-Legendre para estimar o valor numérico de xCM .
2.5- Em um laboratório foram realizadas medidas experimentais da posição de uma partícula em
relação ao tempo e os resultados encontram-se registrados na tabela abaixo.Construa um scrip
para Matlab definindo dois vetores t e x com os dados ao lado. Em seguida
2.5. EXERCÍCIOS 27
t (s) x (m) t (s) x (m)
0.5000 0.1947 3.0000 2.0082
0.9000 0.5165 3.4000 2.1118
1.4000 0.9733 4.0000 2.1654
1.9000 1.3961 4.5000 2.1343
2.5000 1.7907 5.0000 2.0521
1. use a aproximação numérica centrada da derivada
primeira para estimar o vetor velocidade a partir
dos dados fornecidos.
2. Aplique o mesmo procedimento para o vetor veloci-
dade com o intuito de estimar a aceleração da par-
tícula.
3. Faça um gráfico de x(t) , outro de v(t) e finalmente
um para a(t).
2.6- As velocidades instantâneas de um foguete foram transmitidas, via rádio, para a estação de lan-
çamento e os valores registrados encontram-se na tabela abaixo.
t (s) v (m/s) t (s) v (m/s)
0 0 96 747.0
12 31.5 108 987.3
24 72.2 120 1301.8
36 138.9 132 1481.6
48 228.0 144 1489.6
60 321.7 156 1504.6
72 429.5 168 1526.8
84 566.1 180 1556.7
O foguete cheio de combustível possui uma massa de
2, 03× 106 Kg.
Faça um script para Matlab/Octave para
(a) estimar por integração numérica (R. Simpson) o des-
locamento total do foguete;
(b) usando derivação numérica (centrada) estime os valo-
res das acelerações e
(c) faça um gráfico das acelerações em função do tempo.
2.6- A energia potencial gravitacional de um sistema de dois corpos pode ser expressa por
U(r) = −GMm
r
onde G = 6, 7 × 10−11Nm2/kg2 é a constante universal da gravitação e r é a distância entre os
corpos de massas M = 6, 0 × 1024Kg e m = 80Kg. Sabemos ainda que a força gravitacional
pode ser obtida a partir do gradiente de U(r). Considerando o exposto e que r = 6, 4× 106 + z,
onde z é a altitude em metros:
a) determine analiticamente a expressão da força gravitacional que atua nos corpos;
b) desenvolva um script para Matlab em que os dados de entrada sejam as massas; que crie um
vetor contendo os valores de r igualmente espaçados; que faça um gráfico superposto de
F (r) analítico (pontos) e de F (r) numérico (linha cheia) usando a aproximação centrada da
derivada primeira de U ;
c) refaça os cálculos para derivada primeira de ordem superior (centradas Tab 2.3) e armazene
os resultados.
d) Calcule o erro numérico entre a aproximação da derivada primeira e o resultado analítico de
F e faça um gráfico superposto de F × r e do Erro × r para comparação.
2.7- Desprezando a resistência do ar faça um script que calcule numericamente o trabalho W da
força gravitacional ao deslocar um paraquedista com massa de 80 Kg de uma altitude z igual
a 7500 m até uma altitude final de 300 m considerando que (a) a aceleração da gravidade é
constante e igual a 9, 81m/s2 e (b) considerando que a gravidade varie com a altitude segundo
a expressão
g(z) = G
M
(R+ z)2
,
onde R = 6, 4 × 106m é o raio médio da Terra e M = 6, 0 × 1024Kg é a massa da Terra. (c)
Calculea diferença relativa entre os resultados: |Wb −Wa|/Wb.
28 CAPÍTULO 2. INTEGRAÇÃO E DIFERENCIAÇÃO NUMÉRICAS
2.8- Segundo RESNICK; HALLIDAY (1983, cap. 8) a solução completa de problemas unidimensionais
com forças dependentes da posição pode ser descrita pela equação
t2 − t1 = ±
√
m
2
ˆ x2
x1
dx√
E − U(x)
onde m é a massa da partícula, E a energia mecânica (constante) do sistema e U(x) a energia
potencial do sistema. Os tempos t1 e t2 estão relacionados com as posições x1e x2 da partícula,
respectivamente. O sinal + ou − deve ser escolhido segundo o sinal da velocidade da partícula em
eixo orientado. Sendo a energia do sistema invariante ocorre divergência do integrando quando
a energia cinética tende a zero: E −U(x) −→ 0. As regras do trapézio e de Simpson não podem
ser usadas para encontrar o valor de t2 pois o integrando não é definido nas situações em que
E = U(x). A integração numérica por quadratura pode contornar este problema visto que os
pontos de integração sempre estão localizados no interior do intervalo a ser integrado.
a) Considere o sistema massa mola com m = 0, 2Kg, constante elástica k = 30N/m e que em
t1 = 0s a posição inicial da partícula seja x1 = 0, 10cm e que a mesma parta do repouso.
Construa um script para Matlab/OCTAVE em que o limite superior varie de 0,10 m até -0,10
m com incremento de 0,01 m para obter a resposta t2. Este processo permitirá construir dois
vetores x2 e os respectivos tempos t2 associados. (Observ.: dentro do laço para o cálculo de
t2 a integral deve ser resolvida numericamente por um dos processos de integração discutidos
neste capítulo)
b) Compare os resultados obtidos no item anterior com os valores advindos da equação x(t) =
A cos(
√
k
m t) sendo A a amplitude do movimento. Para isso use o vetor t2 na equação e faça
um gráfico superposto dos resultados, i.e., use a instrução plot(t2, x2, ’ro’, t2, x, ’b’).
Capítulo 3
Raízes de Funções de uma Variável
Frequentemente desejamos encontrar o valor de uma variável x que anula ou satisfaz uma determinada
equação, por exemplo: x + 2 = 0; x2 + 3x − 4 = 0; sen(x) = 2x. Existem vários métodos numéricos
que podem ser aplicados para tal fim, dentre eles os clássicos a saber: o Método da Iteração Linear, o
da Bissecção e o de Newton Raphson, Régula Falsi (ou Mét. da Secante).
Um teorema muito utilizado na busca de raízes de funções é o Teorema de Bolzano que afirma
“Se f(x) é uma função contínua sobre um intervalo fechado [a, b] e f(a) e f(b) têm sinais
contrários, então, existe pelo menos um ponto c ∈ [a, b], tal que, f(c) = 0”.
Devemos ter em mente a possibilidade de mais de uma raiz da função f(x) no intervalo mencionado.
Outro caso ocorre quando os sinais de f(a) e f(b) são iguais: neste caso ou existe um número par de
raízes ou nenhuma raiz de f(x) no intervalo [a, b].
3.1 Método da Iteração Linear
Para encontrar o valor de x que zera a função f(x) isolamos x da função escrevendo uma equação
alternativa na forma x = g(x):
• se f(x) = 2x2 − x+ 1 temos: x = 2x2 + 1 = g(x)
• se f(x) = sin(x) + 2x temos: x = −0.5 sin(x) = g(x)
O processo numérico para determinar o valor
(aproximado) da raiz da função f(x) é simples
basta escrever
xi+1 = g(xi) (3.1)
e inserir um valor inicial para o x1 na fun-
ção g(x). O resultado proveniente da função
é atribuído a x2. O valor de i é incrementado
de uma unidade de modo que podemos subs-
tituir em g(x) o valor anteriormente obtido,
x2. Este processo Itera até a precisão numé-
rica desejada:
Passo 1 dado x1 obtemos x2 = g(x1)
Passo 2 com x2 obtemos x3 = g(x2), ...
Passo n com xn obtemos xn+1 = g(xn)
Figura 3.1: Ilustração gráfica do Método da
Iteração Linear.
29
30 CAPÍTULO 3. RAÍZES DE FUNÇÕES DE UMA VARIÁVEL
Observaçao: Nos métodos iterativos usamos com freqüência um ou dois critérios de parada: se
desejamos que o valor da raiz possua precisão de três algarismos significativos então
caso |xi+1 − xi| ≤ � pare o cálculo
onde � = 10−3 ≡ 1.0e − 3. Outro critério de parada pode ser pelo número de iterações,
ou seja, quando n > N - número máximo de iterações. Dependendo do chute inicial x1 o
método pode não convergir para o valor da raiz da função f(x) é o programa terminará
pelo número máximo de iterações estipulado.
3.1.1 Convergência
Suponha que escrevamos a equação f(x) = 0 na forma x = g(x) sendo g(x) a função de iteração. Para
haver convergência do método existem três condições necessárias mas não suficientes a saber:
1. g(x) e g′(x) devem ser contínuas no intervalo [a, b];
2. |g′(x)| < 1 ∀x ∈ [a, b]; a derivada de g deve ser inferior a 1 para todo x no intervalo
3. x1 ∈ [a, b]; o chute inicial deve estar dentro do intervalo considerado;
3.2 Método de Newton-Raphson
Seja f(x) uma função contínua no intervalo [a, b] e
c o único zero da função no intervalo. Se a primeira
derivada e a segunda derivada também são contí-
nuas em [a, b] escrevemos a expansão em Série de
Taylor da função na forma truncada
f(xi + ∆xi) ≈ f(xi) + 1
1!
f ′(xi)∆xi
ou ainda, igualando a zero
∆xi = − f(xi)
f ′(xi)
−→ xi+1 = xi − f(xi)
f ′(xi)
(3.2)
que fornece um método iterativo para a estimativa
do zero da função f(x) no intervalo definido [a, b].
Um problema deste método ocorre quando f ′(xi)
é zero - caso que deve ser evitado pois causa uma
indeterminação no segundo termo à direita da igual-
dade na equação.
Figura 3.2: Ilustração gráfica do Método de
Newton-Raphson
3.2.1 Convergência
Considerando a função f(x) no intervalo [a, b] as condições suficientes de convergência do método de
Newton-Raphson são:
1. f(a).f(b) < 0
2. f ′(x) 6= 0, ∀x ∈ [a, b]
3. f ′′(x) > 0 ou f ′′(x) < 0 (i.e. f ′′(x) não muda de sinal em [a, b])
4.
∣∣∣ f(a)f ′(a) ∣∣∣ < (b− a) e ∣∣∣ f(b)f ′(b) ∣∣∣ < (b− a)
para qualquer valor inicial de x ∈ [a, b].
3.3. MÉTODO DA BISSECÇÃO 31
3.3 Método da Bissecção
Para se aproximar de uma raiz c, o princípio da bisseção consista em reduzir o intervalo inicial testando
o sinal de f(x) para o ponto médio do intervalo, i.e. x = (a+ b)/2. Considerando um intervalo inicial
[a, b] usamos
• Se f(a).f(a+b2 ) < 0 o intervalo [a, b] é substituído por [a, (a+ b)/2]
• Se f(b).f(a+b2 ) < 0 o intervalo [a, b] é substituído por [(a+ b)/2, b]
com o novo intervalo calcula-se novamente o valor médio de x, realiza-se o teste e redefine-se o intervalo
até que o critério de parada seja satisfeito. A desvantagem deste método é a sua lenta convergência
para obter uma aproximação da raiz da equação (x = c). Outra desvantagem é a necessidade de
fornecer ao método dois pontos, a e b, com valores de f(a) e f(b) possuindo sinas contrários.
3.4 Método da Secante (Régula Falsi)
Método usado para obter zeros de funções sua-
ves. Neste método os pontos de cruzamento da
secante de uma aproximação são usados como
pontos de base inicial para uma nova secante,
mais próxima do zero real. É uma modificaçao
do Método de Newton visto que a derivada da
função f ′(x) é substituida por uma aproxima-
ção numérica.
Regula falsi: caso especial do Método da Se-
cante, em que se condiciona que o zero pro-
curado fique sempre entre os pontos de cruza-
mento da secante utilizada na iteração ante-
rior:. Diferente do Método da Secante neste o
processo de busca do zero da função inicia-se
com fa · fb < 0, obrigatóriamente.
Para aplicar o método iterativo usamos a se-
guinte expressão
xi+1 = xi − f(xi) xi − xi−1
f(xi)− f(xi−1)
Figura 3.3: Ilustração gráfica do Método da
Secante.
e, observando melhor, verificamos que a fração à direita da igualdade é o inverso da aproximação
numérica para a derivada primeira (Série de Taylor).
O inconveniente deste método é a escolha de dois valores iniciais para iniciar o processo: se i = 0
necessitamos de x0 e x−1 um à esquerda da raiz e outro à direta dela. Porém se a relação f(a).f(b) < 0
for satisfeita podemos escolher x0 = ae x−1 = b.
3.5 Exercícios
3.1- Dois automóveis partem simultaneamente de cidades vizinhas A e B afastadas de 100 Km em
linha reta. As equações horárias de suas posições são
xa(t) = 12t+ 4t
2,
xb(t) = 10
5 − 6t2,
respectivamente, com a origem do sistema de coordenadas sobre a cidade A. Para encontrar o
tempo em que os móveis se cruzam na estrada devemos igualar as equações, i.e., fazer xa = xb.
32 CAPÍTULO 3. RAÍZES DE FUNÇÕES DE UMA VARIÁVEL
(a) Encontre analiticamente o tempo decorido entre as suas partidas e o instante em que se
cruzam na estrada. No prompt do Matlab defina t=0:1:120; xa= ..., xb=... e faça um
gráfico com o comando plot(t, xa, t, xb). Compare o resultado encontrado algébricamente
com o tempo no ponto que as curvas se cruzam. No gráfico é possível estimar a que distância
da cidade A os carros se encontrarão.
(b) implemente o Método da Iteração em um script para Matlab/Octave para calcular este
tempo. Confronte o resultado com o obtido algébricamente.
(c) Crie outro script para encontrar o tempo decorrido através do método de Newton-Raphson
e verifique se está de acordo com o resultado analítico.
(d) Sem muitas alterações o script do item (c) pode ser modificado para obter o tempo pelo
Método da Secante. Faça isso e compare o tempo obtido numéricamente com o resultado
analítico.
(e) Finalmente escreva, execute e compare o tempo numérico com o analitico implementando o
Método da Bissecção em um novo script para Matlab/OCTAVE.
3.2- Considere o automóvel A do exercício anterior com uma equação horária mais complexa, a saber
xa(t) = 12 t+ t
2 exp(0, 05 t).
Neste caso não é possível obter analiticamente o tempo decorrido entre a partida dos móveis e o
cruzamento dos mesmos na rodovia.
a) No prompt do Matlab defina um vetor t=0:1:70, xa=.... e xb=.... Faça um gráfico com
o comando plot para estimar em que tempo eles se cruzarão e a posição do cruzamento a
partir da cidade A.
b) Altere a equação do automóvel A nos códigos dos itens (b), (c), (d) e (e) do exercício anterior.
Encontre o tempo numericamente por cada método comparando os resultados com o valor
estimado gráficamente.
3.3- No capítulo anterior (secção (2.2.1)) foram mencionados os Polinômios de Legendre. A quadratura
Gauss-Legendre de N pontos exige o conhecimento de N raízes e pesos para sua implementação
e uma função GaussLegendre.m foi apresentada para este fim. Desenvolva um script em que
permita o usuário entrar com o valor de N entre 1 e 4 retornando os pontos da quadratura (zeros
dos polinômios) e seus respectivos pesos usando:
a) o método da Iteração Linear;
b) o método de Newton-Raphson.
(Observação: as raízes e os pesos dos polinômios de Legendre encontram-se limitados ao
intervalo [-1, 1].)
(Sugestões: use a estrutura de decisão if N==1... elseif N==2 ...else .. ...end; defina o
polinômio usando o comando inline (ex. P=inline(’z’) se N=1) e a expressão do peso W (zi)
em cada caso da estrutura. Encontre os zeros e calcule os respectivos pesos.)
3.4- A força de interação entre duas moléculas ou dois átomos neutros pode ser obtida a partir do
potencial de Lennard-Jones
V (r) =
A
r12
− B
r6
,
onde são constantes que dependem do menor valor do potencial V0 e da distância entre moléculas
em que o potencial se anula r0. Considere as distâncias r em Angstrons e os valores A =
2, 27× 105eV·m12 e B = 104, 5eV·m6.
a) Encontre o valor de r0 usando Newton-Raphson ;
3.5. EXERCÍCIOS 33
b) Obtenha a expressão analítica da força F (r) a partir do potencial acima. Use um dos métodos
apresentados neste capítulo para calcular o valor de r em que a força se anula (mínimo do
potencial).
3.5-
34 CAPÍTULO 3. RAÍZES DE FUNÇÕES DE UMA VARIÁVEL
Capítulo 4
Ajuste de Curvas
4.1 Mínimos Quadrados
O ajuste dados provenientes de medidas experimentais com o intuito de prever valores ou leis Físicas
que representem o fenômeno em questão são comuns na Física Experimental. Ao longo desta seção
iremos assumir um conjunto conhecido de dados (xi, yi), com i = 1, 2, . . . , N .
Uma técnica muito empregada minimiza o quadrado dos resíduos entre as observações (yi) e os valores
preditos por uma função de ajuste (f(xi)), i. e.
min
{
R2
}
= min
{
N∑
i=1
[f(xi)− yi]2
}
. (4.1)
Esta minimização pode ser obtida após a definição do grau do polinômio (p) que deve ajustar os dados.
Para esclarecer o método dos mínimos quadrados vamos admitir um polinômio de grau p na forma
f(x) = c0 + c1x+ c2x
2 + · · ·+ cpxp (4.2)
onde desejamos encontrar os coeficientes {c0, c1, . . . , cp}. Como visto na disciplina de cálculo, o
mínimo de R2 pode ser obtido pela derivada da função R2 com respeito aos coeficientes do polinômio,
i.e.
∂(R2)
∂c0
= 2
∑
[f(xi)− yi] ∂f(xi)
∂c0
= 2
∑
[f(xi)− yi] = 0
∂(R2)
∂c1
= 2
∑
[f(xi)− yi] ∂f(xi)
∂c1
= 2
∑
xi [f(xi)− yi] = 0
∂(R2)
∂c2
= 2
∑
[f(xi)− yi] ∂f(xi)
∂c2
= 2
∑
x2i [f(xi)− yi] = 0
∂(R2)
∂cp
= 2
∑
[f(xi)− yi] ∂f(xi)
∂cp
= 2
∑
xpi [f(xi)− yi] = 0
substituindo f(xi) e expandindo temos
c0
∑
1 + c1
∑
xi + c2
∑
x2i + · · ·+ cp
∑
xpi =
∑
yi
c0
∑
xi + c1
∑
x2i + c2
∑
x3i + · · ·+ cp
∑
xp+1i =
∑
xiyi
c0
∑
x2i + c1
∑
x3i + c2
∑
x4i + · · ·+ cp
∑
xp+2i =
∑
x2i yi
...
...
...
c0
∑
xpi + c1
∑
xp+1i + c2
∑
xp+2i + · · ·+ cp
∑
x2pi =
∑
xpi yi
35
36 CAPÍTULO 4. AJUSTE DE CURVAS
e que pode ser expresso na forma matricial a saber
N
∑
xi
∑
x2i · · · · · ·
∑
xpi∑
xi
∑
x2i
∑
x3i · · · · · ·
∑
xp+1i
...
...
...
...∑
xpi
∑
xp+1i
∑
xp+2i · · · · · ·
∑
x2pi


c0
c1
c2
...
cp
 =

∑
yi∑
xiyi∑
x2i yi
...∑
xpi yi
 . (4.3)
Podemos ainda renomear as matrizes e vetores e reescrever o sistema de equações em uma forma
compacta
M c¯ = d¯ (4.4)
cuja solução desejada é obtida pela multiplicação da inversa de M pelo vetor d¯,
c¯ = M
−1
d¯.
4.1.1 Exemplos
Nos vários ramos da física é comum encontrarmos leis, relativamente simples, que prevêem os compor-
tamentos dos sistemas. Tais leis são funções como retas e parábolas, exponenciais, logaritmos entre
outras. Algumas dessas funções não são polinômios porém com uma simples transformaçao de variáveis
tais funções podem ser linearizadas.
Exemplo 1: Ajuste de dados através de uma função linear do tipo: y = A+Bx.
Neste caso possuímos uma tabela de dados na forma
xi yi
x1 y1
x2 y1
...
...
xN yN
e desejamos ajustar uma reta
pelos dados experimentais (xi, yi). Pelo que foi exposto na seção 4.1 o sistema a ser resolvido é[
N
∑
xi∑
xi
∑
x2i
] [
A
B
]
=
[ ∑
yi∑
xiyi
]
(4.5)
Sendo um sistema 2× 2 é possível obter uma solução analítica simples para os coeficientes dese-
jados do polinômio A e B.
O seguinte algoritmo contém as etapas básicas para o cálculo dos coeficientes da reta
1- limpar memória
2- definir os vetores x e y com os dados provenientes das medidas
3- definir N ← número de medidas
4- definir soma_x, soma_x2, soma_y, soma_xy como zero
5- para i de 1 até N faça
soma_x ←soma_x + x (i)
soma_x2 ←soma_x2 + x (i) * x (i)
soma_y ←soma_y + y(i)
soma_xy ←soma_xy + x(i) * y(i)
fim do laço 5
6- gerar a matriz M e o vetor d (sistema M·c = d)
7- obter o vetor de coeficientes c = M−1· d
8- mostrar os coeficientes.
Exemplo 2: Ajuste de dados através de uma função logarítmica na forma:
y = A+B lnx.
4.1. MÍNIMOS QUADRADOS 37
Para obtermos os coeficientes A e B reescrevemos a equação na forma y = A+Bw onde w = lnx.
O sistema que permite o ajuste dos dados é similar ao do Ex1, i.e.[
N
∑
wi∑
wi
∑
w2i
] [
A
B
]
=
[ ∑
yi∑
xiyi
]
ou seja, basta criar mais uma coluna de dados na tabela contendo ao invés de xi os valores de
wi= lnxi.
Exemplo 3: Ajuste de dados por uma função exponencial na forma:
y = A exp(Bx).
Para obtermos os coeficientes A e B do ajuste reescrevemos a equação na forma ln y = lnA+Bx
ou ainda z = a+bx. Assim a tabela de dados necessita de uma coluna com os valores de ln yi = zi.
O procedimento do Ex1 é então repetido de modo a obter os coeficientes a e b por[
N
∑
xi∑
xi
∑
x2i
] [
a
b
]
=
[ ∑
zi∑
xizi
]
.
Os reais coeficientes da função exponencial serão então A = exp(a) e B = b.
Exemplo 4: Ajuste de dados por uma lei de potência na forma:
y = AxB.
Neste caso aplicamos a função “ln” na equação para obtermos ln y = lnA + B lnx ou ainda
z = a+ bw. O sistema resultante será um misto dos exemplos 2 e 3 ou seja[
N
∑
w∑
wi
∑
w2i
] [
a
b
]
=
[ ∑
zi∑
wizi
]
.
Os coeficientes desejados para o ajuste serão A = exp(a) e B = b.
Exemplo 5:
Em uma experiência como pêndulo
simples conseguimos os seguintes perío-
dos em função do comprimento (tabela
ao lado). Após efetuar os procedimentos
indicados nos exemplos de 1 a 4 obtive-
mos os seguintes coeficientes:
Tabela 4.1 - dados de comprimento e pe-
ríodo em experiência de pêndulo simples.
Li(m) Ti(s)
0.50 1.43
1.00 2.00
2.00 2.85
3.00 3.48
4.00 3.98
- ajuste linear: A= 1.23 e B= 0.722 y = 1, 23 + 0, 722x
- ajuste logaritmo: A= 2.14 e B= 1.23 y = 2, 14 + 1, 23 lnx
- ajuste exponencial: A=exp(a)=1.43 ; b=B=0.281 y = 1, 43 exp(0, 281x)
- ajuste potência de x: A=exp(a)=2.01 ; b=B=0.495 y = 2, 01x0,495
O gráfico abaixo ilustra o resultado dos ajustes. Podemos observar que a Lei de Potência forneceu
o melhor ajuste entre os apresentados. Entretanto a potência deveria ser igual a 0,5. Isto indica
que as medidas experimentais não foram realizadas com a necessária precisão fornecendo um erro
de (0,500-0,495)/0,500= 0,01 = 1%.
O valor da gravidade local pode ser estimado a partir do coeficiente A fazendo g = (2pi/A)2 =
9, 77m/s2. Assumindo o valor real de 9,81m/s2 o erro na estimativa de g é (9,81-9,77)/9,81=0,004
= 0,4%.
38 CAPÍTULO 4. AJUSTE DE CURVAS
Figura 4.1: Diferentes ajustes obtidos com os dados provenientes de uma experiência com o Pêndulo
Simples.
Exercício 2.1: Monte o sistema p/ obter os coeficientes A e B considerando um ajuste linear (reta)
do seguinte conjunto de dados: x = [0.00, 0.20, 0. 50, 0.90, 1.40] e y = [2.32, 1.46, 0.40, -0.88,
-2.32]. Resolva o sistema analiticamente, i. e., forneça a equação y = A + B*x.
Exercício 2.2: Faça um ’script’ usando a linguagem do Matlab para realizar um ajuste linear dos
pontos fornecidos no Exercício 2.1. No final do código faça um gráfico com os dados da tabela
(círculos vermelhos) e, sobreposto, o ajuste conseguido (linha cheia em azul).
4.1.2 Mínimos Quadrados Ponderados
Imagine que tenhamos a disposição um cronômetro cuja menor leitura seja de um segundo. Ao reali-
zarmos medidas com este cronômetro os resultados devem ser expressos com a devida imprecisão do
equipamento: valor medido±0, 5s. Em certos equipamentos a incerteza da medida é um percentual do
valor medido (e.g., 2, 5± 0, 02 · 2, 5 Volts), ou seja, as incertezas nos registros de voltagem são de 2%
dos valores medidos.
Vamos assumir um conjunto de N medidas (xi, yi) e que os valores de yi possuam imprecisões co-
nhecidas σi, respectivamente. Para uma dada função ajuste f(xi, c) polinomial de grau p desejamos
encontrar os coeficientes c = {co, c1, . . . , cp} que minimizam a soma dos quadrados dos resíduos pon-
derados
R2w =
N∑
i=1
[
f(xi, c)− yi
σi
]2
=
N∑
i=1
w2i [f(xi, c)− yi]2 , (4.6)
onde wi = 1/σi sendo que w2i representa o quadrado do peso atribuído ao i-ésimo resíduo. Assim,
quanto menor for a incerteza de uma dada medida maior será o seu peso na estimativa dos coeficientes
da função ajuste.
Sabemos que as medidas experimentais podem sofrer flutuações por imperícia do experimentador ao
manusear o equipamento de medida, má calibração do equipamento, flutuações térmicas do laboratório
ou mesmo na rede elétrica. Estas flutuações (aleatórias) alteram as medidas além do normal aumen-
tando suas incertezas. Na figura abaixo duas medidas se afastam consideravelmente do ajuste linear
estipulado para os dados. Usando Mínimos Quadrados convencional os valores dos coeficientes da reta
4.1. MÍNIMOS QUADRADOS 39
são alterados por estes dois desvios exagerados. Uma forma de contornar este problema seria refazer
as duas medidas ou, na impossibilidade, descartá-las antes do ajuste.
Figura 4.2: Comparação entre ajustes realizados por Mínimos Quadrados com e sem ponderação.
Na tentativa de minimizar os efeitos de medidas suspeitas no ajuste podemos optar pelo Método dos
Mínimos Quadrados Ponderados com uma técnica iterativa. O procedimento inicial é usar Mínimos
Quadrados convencional para obter um primeiro conjunto de coeficientes da função ajuste, c(1) =
{c0, c1, . . . , cp}(1). Na seqüência os resíduos são computados por ri = f(xi, c) − yi e então usados
para definir os pesos wi. Os pesos recém computados são usados com MQP para estimar os novos
coeficientes da função ajuste, c(2) = {c0, c1, . . . , cp}(2). O processo pode iterar mais algumas vezer até
que um critério de convergência (ou parada) seja satisfeito.
Para facillitar as deduções mudaremos a formulação adotada para uma matricial muito utilizada em
artigos e livros que tratam de ajuste, regressão e inversão.
Desejamos minimizar R2w e para isso usando uma matriz de pesos diagonal Wii = wi e a definição
f(xi, c) =
[
1 xi x
2
i ... x
p
i
]

c0
c1
c2
...
cp
 =
m∑
j=1
Xijcj−1, com m = p+ 1
reescrevemos a Eq.(4.6) a saber
R2w =
N∑
i=1
Wii
 m∑
j=1
Xijcj−1 − yi
 ·
Wii
 m∑
j=1
Xijcj−1 − yi

Como cada termo entre colchetes é o i-ésimo elemento de um vetor, temos
R2w = [W (Xc− y)]T [W (Xc− y)] . (4.7)
Usando as identidades (XY)T = YTXT e WTW = W2, obtemos
R2w = (Xc− y)T W2 (Xc− y) ,
= (Xc)TW2 (Xc− y)− yTW2 (Xc− y) ,
= cTXTW2Xc− cTXTW2y− yTW2Xc + yTW2y,
R2w = c
TXTW2Xc− 2yTW2Xc + yTW2y.
40 CAPÍTULO 4. AJUSTE DE CURVAS
Com o auxílio das identidades
dxTx
dx
= 2x ;
daTx
dx
= a ;
dxTAx
dx
= 2Ax
derivamos R2w em relação aos parâmetros e igualamos a zero, ou seja
dR2w
dc
= 2XTW2Xc− 2 (yTW2X)T + 0 = 0,
2XTW2Xc− 2XTW2y = 0,
XTW2Xc = XTW2y,
=⇒ c = (XTW2X)−1 XTW2y. (4.8)
A Eq.(4.8) fornece os coeficientes que minimizam R2w pelo método dos Mínimos Quadrados Ponde-
rados. O caso particular em que a matriz diagonal W2 é igual a matriz identidade os coeficientes
retornados pela equação serão os mesmos obtidos pelo método dos Mínimos Quadrados convencional:
c =
(
XTX
)−1 XTy = M−1d, com M = (XTX) e d = XTy.
4.2 Ajuste por uma função qualquer (Newton e Gauss-Newton)
Em alguns casos os dados experimentais não possuem um bom ajuste com as funções exploradas nas
seções precedentes. Exemplos destes casos são
P (w, A, wo, Γ) =
A
(w − wo)2 + Γ2 , potência em exp. de ressonância,
G(x, σ, xo) =
1√
2piσ2
exp(
−(x− xo)2
2σ2
), função Gaussiana,
f(k, λ) =
λk exp(−λ)
k!
, distribuição de Poisson
e cujos gráficos apresentam-se como nas figuras abaixo.
(a)
−5 −4 −3 −2 −1 0 1 2 3 4 5 
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
x0= 0,0; σ
2
 =0,2
x0= 0,0; σ
2
 =1,0
x0= 0,0; σ
2
 =5,0
x0=−2,0; σ
2
 =0,5
(b)
Figura 4.3: Gráficos ilustrando (a) distribuições de Poisson e (b) distribuições Normal (Gaussiana).
Nesta seção será mostrado um procedimento conhecido pelo Método de Newton que permite ajustar
dados experimentais por uma função qualquer - função esta definida previamente - a exemplo das
equações acima apresentadas e normalmente não lineares.
4.2. AJUSTE POR UMA FUNÇÃO QUALQUER (NEWTON E GAUSS-NEWTON) 41
Neste método desejamos minimizar o quadrado dos resíduos como