Programação para as
Ciências Experimentais
2006/7
Teórica 4
Ludwig Krippahl, 2007
Na aula de hoje...

Revisão:

Encontrar o zero de um polinómio

Encontrar o zero de uma função contínua
Encontrar o mínimo de uma função contínua

• function, if, for, while
• Precisão, representação de valores numéricos
Ludwig Krippahl, 2007
2
Funções



Para usar a função
Se retorna uma variável:
x = funcaoqq(arg1, arg2)
Se retorna mais que uma:
[x,y,z] = outrafn(arg1, arg2)
Ludwig Krippahl, 2007
3
Funções

O que o Octave faz
• funcaoqq – não há nada com este nome em
•
•
memória.
Procura ficheiro funcaoqq.m
Nesse ficheiro executa a função
Ludwig Krippahl, 2007
4
Funções

O que nós fazemos
• Criamos o ficheiro funcaoqq.m
• Nesse ficheiro declaramos a função:
function res=funcaoqq(arg1,arg2)
....
endfunction
Ludwig Krippahl, 2007
5
Funções
function res = funcaoqq( arg1 , arg2 )
Indica que é a declaração de uma
função.
Ludwig Krippahl, 2007
6
Funções
function res = funcaoqq( arg1 , arg2 )
Nome da variável com o valor a
devolver
Ludwig Krippahl, 2007
7
Funções
function [res1, res2] = funcaoqq( arg1 ..
Se devolve vários valores, usamos
um vector de variáveis
Ludwig Krippahl, 2007
8
Funções
function res = funcaoqq( arg1 , arg2 )
Variáveis (locais) para onde são
copiados os valores dados à
função como argumentos
Ludwig Krippahl, 2007
9
Funções
function res = funcaoqq( arg1 , arg2 )
Todas a variáveis declaradas aqui
e no corpo da função são locais.
Só existem dentro da função e não
afectam nem são afectadas por
variáveis externas.
endfunction
Ludwig Krippahl, 2007
10
Controlo de execução: if
if condição
Isto é executado se a condição for
verdadeira (não for 0)
else
Isto é executado se a condição for
false (0). O else é opcional
endif
Ludwig Krippahl, 2007
11
Controlo de execução: while
while condição
Executado enquanto a condição
for verdadeira (não for 0).
É preciso garantir que dentro do
ciclo o valor da condição muda,
senão o ciclo não acaba...
endwhile
Ludwig Krippahl, 2007
12
Controlo de execução: while
for var=vector
Esta parte é repetida uma vez para
cada valor no vector.
A variável var toma cada um dos
valores do vector a cada iteração.
endfor
Ludwig Krippahl, 2007
13
Função: polinomio

Polinómio:
Y= k1 + k2*x + k3*x2 + k4*x3 ...
Ludwig Krippahl, 2007
14
Função: polinomio

Polinómio:
Y= k1 + k2*x + k3*x2 + k4*x3 ...
Coeficientes num vector:
[ k1 , k2 , k3 , k4 ... ]
Ludwig Krippahl, 2007
15
Função: polinomio
function y=polinomio(coefs,x)
xx=x;
y=coefs(1);
for f=2:length(coefs)
y=y+coefs(f)*xx;
indentação, torna
xx=x*xx;
mais legível.
endfor
endfunction
Ludwig Krippahl, 2007
16
Função: polinomio
function y=polinomio(coefs,x)
xx=x;
y=coefs(1);
for f=2:length(coefs)
y=y+coefs(f)*xx;
Erro se x for um
xx=x*xx;
vector.
endfor
endfunction
Ludwig Krippahl, 2007
17
Função: polinomio
function y=polinomio(coefs,x)
xx=x;
y=coefs(1);
for f=2:length(coefs)
y=y+coefs(f)*xx;
já funciona com
xx=x.*xx;
vectores
endfor
endfunction
Ludwig Krippahl, 2007
18
Função: polinomio (exemplos)

Calcular o valor de
• y = 2+3x-x2
para x=3
octave:16> polinomio([2,3,-1],3)
ans = 2

Traçar o gráfico de
• y = 2+3x-x2 entre -10 e 10:
• plot(-10:10, polinomio([2,3,-1], -10:10))
Ludwig Krippahl, 2007
19
Uma raiz de um polinómio

Método da bissecção.
Ludwig Krippahl, 2007
20
Uma raiz de um polinómio


y = 0.3+x – x2 + x3
x = -0.23309
Ludwig Krippahl, 2007
21
Uma raiz de um polinómio

Começamos com um intervalo que inclui
o zero: [-1,1]
Ludwig Krippahl, 2007
22
Uma raiz de um polinómio

Começamos com um intervalo que inclui
o zero: [-1,1]
+

Os extremos têm sinal diferente:
-
Ludwig Krippahl, 2007
23
Uma raiz de um polinómio

Dividir ao meio: (-1 +1)/2 = 0
+
-
Ludwig Krippahl, 2007
24
Uma raiz de um polinómio

Calculamos y(0)= 0.3
+
+
-
Ludwig Krippahl, 2007
25
Uma raiz de um polinómio

Calculamos y(0)= 0.3
-

+
+
O intervalo com sinais opostos nos
extremos contém o zero
Ludwig Krippahl, 2007
26
Uma raiz de um polinómio
-
Ludwig Krippahl, 2007
-
+
+
27
Uma raiz de um polinómio
-
Ludwig Krippahl, 2007
-
-
+
+
28
Uma raiz de um polinómio

Quando paramos?
• Quando o intervalo for pequeno
• Quando no ponto médio y próximo de 0.
• Precisão.
Ludwig Krippahl, 2007
29
Uma raiz de um polinómio

Algoritmo:
• Dado: x1, x2, precisão
• Enquanto abs(x2-x1)>precisão repetir
• ym = valor no ponto médio xm
• Se abs(ym)<precisão, pára
• Caso contrário, escolhe o intervalo x1,xm ou xm,x2
onde os sinais sejam opostos
Ludwig Krippahl, 2007
30
Representação de números

A que precisão podemos ir?
Ludwig Krippahl, 2007
31
Representação de números

Um número no Octave é representado
com 64 bits (double precision):
• Sinal (+, -) : 1 bit
• Expoente: 11 bits
• Mantissa: 52 bits

(sinal) mantissa * 2 expoente
Ludwig Krippahl, 2007
32
Representação de números
• Expoente: 11 bits,
• Mantissa: 52 bits
• (sinal) mantissa * 2 expoente

Máximo valor:

Precisão (épsilon)
• 1.7977x 10308 (realmax, 1.7977e308 )
• 2.2204x 10-16 (eps, 2.2204e-16)
Ludwig Krippahl, 2007
33
Representação de números

Precisão (épsilon)
• 2.2204x 10-16 (eps, 2.2204e-16)
• O menor número que somado a 1 dá um
resultado diferente de 1:
octave:17> (1+eps)==1
ans = 0
octave:18> (1+eps/2)==1
ans = 1
Ludwig Krippahl, 2007
34
Representação de números

Importante:
• Todos os dados no computador são
sequências de bits.
• A memória é limitada (64 bits para os
números), por isso a precisão é limitada.
• Normalmente não há problema, mas atenção
aos arredondamentos:
octave:20> sqrt(2)^2==2
ans = 0
Ludwig Krippahl, 2007
35
O zero de uma função
Suponhamos que uma função y=f(x)
pode ser especificada por um vector de
parâmetros (constantes) e pelo nome.
 e.g:
function y=polinomio(coefs,x)

Ludwig Krippahl, 2007
36
O zero de uma função
A nossa função genérica será.
y = nome(params,x)


Para a função que encontra o zero
temos que enviar o nome, os
parâmetros, o intervalo, e a precisão.
function xm=zerofn(func,params,x1,x2,prec)
Ludwig Krippahl, 2007
37
O zero de uma função
Para avaliar a função func usamos a função do
Octave feval:
feval(nome,arg1,arg2, arg3) é o mesmo que
nome(arg1, arg2, arg3)

octave:22> sin(1)
ans = 0.84147
octave:23> feval("sin",1)
ans = 0.84147
Ludwig Krippahl, 2007
38
O zero de uma função
Em vez de:
y1=polinomio(coefs,x1);
y2=polinomio(coefs,x2);
 Fica
y1=feval(func,params,x1);
y2=feval(func,params,x2);

Ludwig Krippahl, 2007
39
O zero de uma função

Para calcular uma raíz do polinómio:
• z=zerofn("polinomio",coefs,-1,1,0.0001)
Ludwig Krippahl, 2007
40
O mínimo de uma função

Método da razão dourada
Ludwig Krippahl, 2007
41
O mínimo de uma função

Tal como “encurralámos” a raiz num
intervalo, vamos fazer o mesmo com o
mínimo, mas precisamos de 3 pontos:
a
c
b
Ludwig Krippahl, 2007
42
O mínimo de uma função

Se x1<x2<x3 e y2<y1 e y2<y3 então tem
que haver um mínimo local entre x1 e x3
x1
Ludwig Krippahl, 2007
x2
x3
43
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2007
x2
x3
44
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2007
x2
x3
45
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2007
x2
x3
46
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2007
x2
x3
47
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2007
x2
x3
48
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2007
x2
x3
49
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2007
x2
x3
50
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2007
x2
x3
51
O mínimo de uma função

O algoritmo é (novamente) partir os
intervalos, testar, e repetir até que seja
suficientemente pequeno
x1
Ludwig Krippahl, 2007
x2
x3
52
O mínimo de uma função

Guardar sempre os 3 pontos
consecutivos em que o do meio é menor
que os extremos.
x1
Ludwig Krippahl, 2007
x2
x3
53
O mínimo de uma função

Como dividir o intervalo:
• O ideal é manter as proporções. Dividir ao
meio não é bom.
x1
Ludwig Krippahl, 2007
x3
x2
54
O mínimo de uma função

Como dividir o intervalo:
x1
Ludwig Krippahl, 2007
x4
x5
x3
x2
55
O mínimo de uma função

Como dividir o intervalo:
• Escolher o ponto novo no intervalo maior e
• Partir pela razão dourada:
(a+b)/a = a / b
a= 0.618 (a+b)
b= (1-0.618) (a+b)
Ludwig Krippahl, 2007
56
O mínimo de uma função

Detalhes:
• Escolher o lado maior
• if abs(x1-xm)>abs(x2-xm)
• false, lado maior xm x2
x1
Ludwig Krippahl, 2007
x2
xm
57
O mínimo de uma função

Detalhes:
• Calcular novo ponto xn
• xn=c1*xm+c2*x2
x1
Ludwig Krippahl, 2007
xm
(c1=0.618; c2=1-c1)
x2
xn
58
O mínimo de uma função

Detalhes:
• Calcular yn=f(xn)
• yn=feval(func,params,xn);
x1
Ludwig Krippahl, 2007
xm
x2
xn
59
O mínimo de uma função

Detalhes:
• se yn<xn, xn será o novo xm, e xm o novo x1
• ym=yn
• x1=xm
• xm=xn
x1
Ludwig Krippahl, 2007
xm
x2
xn
60
O mínimo de uma função

Detalhes:
• caso contrário xn será o novo x2
x1
Ludwig Krippahl, 2007
xm
x2
xn
61
O mínimo de uma função

Se o lado maior for entre x1 e xm,
mesma coisa, mas trocando o x1 e o
x2...
Ludwig Krippahl, 2007
62
O mínimo de uma função

Como podemos seguir o cálculo:
guardar o xn e yn num vector (pontos).
function [xm,pontos]=minfnpts(func,params,x1,xm,x2,prec)
...
guardar os 3 pontos iniciais:
y1=feval(func,params,x1);
y2=feval(func,params,x2);
ym=feval(func,params,xm);
pontos=[x1,x2,xm;y1,y2,ym];
Ludwig Krippahl, 2007
63
O mínimo de uma função

Como podemos seguir o cálculo:
guardar o xn e yn num vector (pontos).
while abs(x2-x1)>prec
...
...
durante o ciclo guardar cada xn, yn:
pontos=[pontos,[xn;yn]];
Ludwig Krippahl, 2007
64
O mínimo de uma função

Como podemos seguir o cálculo: guardar o xn e yn num
vector (pontos).
Depois de chamada a função, fazer o plot:
coefs=[0.3,-5,1]
v=-10:10;
clearplot
plot(v,polinomio(coefs,v));
traça o polinómio
[x,p]=minfnpts("polinomio",coefs,-10,0,10,0.001)
hold on
plot(p(1,:),p(2,:),"+")
traça os pontos
Ludwig Krippahl, 2007
65
O mínimo de uma função
Ludwig Krippahl, 2007
66
Mais informação:




Numerical Recipes
http://www.nrbook.com/a/bookcpdf.php
Raiz: 9.1
Mínimo: 10.1
Ludwig Krippahl, 2007
67
Próxima aula:

Apresentação do trabalho prático:
• A partir de uma lista de strings com equações
de reacções químicas, constantes de
equilíbrio, e concentrações iniciais, calcular
as concentrações de equilíbrio de todas as
espécies tendo em conta todas as reacções.

Revisões
• function, if, while, for: muito importante ter isto
bem sabido...
Ludwig Krippahl, 2007
68
Dúvidas...
Ludwig Krippahl, 2007
69