Introdu¸cão aos Sistemas Dinâmicos usando Maximamargarete/download/MET200-0/sistemasdinamicos.pdf · Introdu¸cão aos Sistemas Dinâmicos usando Maxima Jaime E. Villate Faculdade

Introducao aos Sistemas Dinamicos usando Maxima

Jaime E. VillateFaculdade de Engenharia da

Universidade do Porto

Janeiro de 2004 (versao 0.2)

Indice

1 Introducao 11.1 Sistemas dinamicos . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 Sistemas Computacionais Algebricos (CAS) . . . . . . . . . . . 21.3 Maxima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41.4 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2 Computacao algebrica usando Maxima 62.1 A interface do Maxima . . . . . . . . . . . . . . . . . . . . . . . 62.2 Entrada e saıda . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.3 Variaveis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82.4 Constantes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92.5 Equacoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.6 Funcoes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.7 Algebra e trigonometria . . . . . . . . . . . . . . . . . . . . . . 162.8 Passar informacao entre sessoes . . . . . . . . . . . . . . . . . . 172.9 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3 Alguns conceitos de calculo 213.1 Derivada de uma funcao . . . . . . . . . . . . . . . . . . . . . . 21

3.1.1 Interpretacao geometrica da derivada . . . . . . . . . . . 223.1.2 Algumas propriedades . . . . . . . . . . . . . . . . . . . 23

3.2 Primitivas e integrais . . . . . . . . . . . . . . . . . . . . . . . . 273.3 Serie de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . 303.4 A transformada de Laplace . . . . . . . . . . . . . . . . . . . . 323.5 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4 Equacoes diferenciais ordinarias 354.1 Equacoes de primeira ordem . . . . . . . . . . . . . . . . . . . . 354.2 Campo de direccoes . . . . . . . . . . . . . . . . . . . . . . . . 37

ii Introducao aos Sistemas Dinamicos usando Maxima

4.3 Problemas com condicao inicial . . . . . . . . . . . . . . . . . . 384.4 Resolucao de equacoes de primeira ordem . . . . . . . . . . . . 39

4.4.1 Um caso simples . . . . . . . . . . . . . . . . . . . . . . 394.4.2 Equacoes autonomas . . . . . . . . . . . . . . . . . . . . 394.4.3 Outras equacoes . . . . . . . . . . . . . . . . . . . . . . 39

4.5 Classificacao das equacoes de primeira ordem . . . . . . . . . . 434.5.1 Equacoes de variaveis separaveis . . . . . . . . . . . . . 434.5.2 Equacoes exactas . . . . . . . . . . . . . . . . . . . . . . 444.5.3 Equacoes lineares . . . . . . . . . . . . . . . . . . . . . . 454.5.4 Equacoes homogeneas . . . . . . . . . . . . . . . . . . . 46

4.6 Equacao de Bernoulli . . . . . . . . . . . . . . . . . . . . . . . . 474.7 Equacoes de segunda ordem . . . . . . . . . . . . . . . . . . . . 48

4.7.1 Equacoes autonomas . . . . . . . . . . . . . . . . . . . . 484.7.2 Equacoes lineares . . . . . . . . . . . . . . . . . . . . . . 51

4.8 Espaco de fase . . . . . . . . . . . . . . . . . . . . . . . . . . . 534.9 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

5 Resolucao numerica de equacoes diferenciais 565.1 Metodo de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . 565.2 Metodo de Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . 625.3 Sistemas de equacoes de primeira ordem . . . . . . . . . . . . . 645.4 Eliminacao de singularidades . . . . . . . . . . . . . . . . . . . 665.5 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

6 Sistemas oscilantes 706.1 O oscilador harmonico simples . . . . . . . . . . . . . . . . . . 706.2 Osciladores amortecidos . . . . . . . . . . . . . . . . . . . . . . 736.3 Osciladores forcados . . . . . . . . . . . . . . . . . . . . . . . . 756.4 O pendulo simples . . . . . . . . . . . . . . . . . . . . . . . . . 766.5 Osciladores acoplados . . . . . . . . . . . . . . . . . . . . . . . 796.6 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

7 Ondas e difusao 857.1 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

8 Estabilidade e bifurcacao 888.1 Teoria linear da estabilidade . . . . . . . . . . . . . . . . . . . . 888.2 Bifurcacao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 918.3 Pontos fixos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

Indice iii

8.4 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

9 Osciladores nao-lineares e caos 989.1 Oscilador cubico, forcado . . . . . . . . . . . . . . . . . . . . . 989.2 Oscilador de van der Pol . . . . . . . . . . . . . . . . . . . . . . 989.3 Equacoes de Lorenz . . . . . . . . . . . . . . . . . . . . . . . . 999.4 Caos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1019.5 Equacoes de Rossler . . . . . . . . . . . . . . . . . . . . . . . . 1019.6 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

10 Sistemas dinamicos discretos 10410.1 Resolucao Numerica de Equacoes . . . . . . . . . . . . . . . . . 10510.2 Representacao grafica . . . . . . . . . . . . . . . . . . . . . . . 10810.3 Pontos fixos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11210.4 Pontos periodicos . . . . . . . . . . . . . . . . . . . . . . . . . . 11310.5 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

11 Sistemas em duas dimensoes e fractais 11811.1 Equacoes de diferencas de segunda ordem . . . . . . . . . . . . 11811.2 Fractais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12011.3 Sistemas iterativos de funcoes . . . . . . . . . . . . . . . . . . . 12211.4 Problemas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

Indice Remissivo 126

iv Introducao aos Sistemas Dinamicos usando Maxima

Lista de figuras

1.1 Espaco de fase de uma esfera ligada a uma mola, mostrandoduas possıveis trajectorias correspondentes a duas oscilacoes. . 3

2.1 Xmaxima, versao 5.9.0. . . . . . . . . . . . . . . . . . . . . . . 72.2 Grafico do polinomio 3x3 + 5x2 − x + 6. . . . . . . . . . . . . . 11

3.1 Grafico de log(x) e a sua derivada, 1/x. . . . . . . . . . . . . . 223.2 Integral de 1/(1 + x2) entre 0 e 1. . . . . . . . . . . . . . . . . . 293.3 Area entre y = 2x2 − 3 e y = 2x. . . . . . . . . . . . . . . . . . 30

4.1 Campo de direccoes de y′ = x2 + 3y, e seis solucoes particulares. 374.2 Aumento da populacao, no modelo logıstico. . . . . . . . . . . . 42

5.1 Aproximacao numerica a solucao de x = x, com diferentes numerosde passos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

5.2 Solucoes particulares da equacao x = t − x2, com x = 3, 1, 0,−0.7 e −0.75, em t = 0. . . . . . . . . . . . . . . . . . . . . . . 60

5.3 Solucao caotica para n = 18000 (esquerda), e desaparicao daregiao caotica usando uma precisao maior n = 26000. . . . . . . 62

5.4 Solucao de x = t− x2 obtida pelo metodo de Runge-Kutta. . . 645.5 Campo de direccoes do sistema x = y, y = −x e curva integral

em (3, 0). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 675.6 Solucao do sistema x = y, y = −x, com x(0) = 3 e y(0) = 0. . . 68

6.1 Esfera suspensa por uma mola vertical. . . . . . . . . . . . . . . 706.2 Evolucao do oscilador harmonico simples, no espaco de fase e no

domınio do tempo. . . . . . . . . . . . . . . . . . . . . . . . . . 726.3 Oscilador harmonico amortecido. . . . . . . . . . . . . . . . . . 736.4 Evolucao do oscilador harmonico com amortecimento fraco. . . 746.5 Oscilacoes forcadas. . . . . . . . . . . . . . . . . . . . . . . . . 76

LISTA DE FIGURAS v

6.6 Pendulo simples e coordenadas polares. . . . . . . . . . . . . . 776.7 Algumas trajectorias do pendulo simples, no espaco de fase e no

domınio do tempo. . . . . . . . . . . . . . . . . . . . . . . . . . 786.8 Dois osciladores acoplados. . . . . . . . . . . . . . . . . . . . . 79

7.1 Perturbacao na pressao do gas num cilindro. . . . . . . . . . . . 857.2 Onda de presao z(t, x) com velocidade no sentido positivo de x. 86

8.1 Pendulo simples ligado a uma mola e diagrama de corpo livre. 898.2 Diagrama de bifurcacao do ponto de equilıbrio do pendulo com

mola e energia potencial. . . . . . . . . . . . . . . . . . . . . . . 928.3 Os tres tipos de pontos fixos, quando os valores proprios sao

reais. Ponto repulsivo, atractivo e de sela. . . . . . . . . . . . . 948.4 Os tres tipos de ponto fixo, quando os valores proprios sao com-

plexos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

9.1 Oscilacoes do sistema de Lorentz para dois valores muito proximosdo valor inicial: x(0) = 5 (vermelho) e x(0) = 5.005 (azul).Parametros: a = 10, b = 8/3, r = 28, y(0) = 5, z(0) = 5. . . . . 99

9.2 Trajectoria no espaco de fase, projectada no plano xz. Osparametros sao os mesmos da figura 9.1, com x(0) = 5. . . . . . 100

9.3 O caso c = 2.3 (com parametros a = b = 0.2) conduz a um ciclode fronteira com perıodo simples. . . . . . . . . . . . . . . . . . 101

9.4 O caso c = 3.3 (com parametros a = b = 0.2) conduz a um ciclode fronteira com perıodo duplo. . . . . . . . . . . . . . . . . . . 102

10.1 Metodo de Newton para aproximacao a uma raız. . . . . . . . . 10610.2 Solucao de xn+1 = cos(xn) com x0 = 2. . . . . . . . . . . . . . 10810.3 Diagrama de degraus para xn+1 = cos(xn) com x0 = 2. . . . . . 10910.4 Solucao do sistema xn+1 = x2

n − 0.2 com valor inicial 1.1 (es-querda) e 1.5 (direita). . . . . . . . . . . . . . . . . . . . . . . . 110

10.5 Orbitas do modelo logıstico com valor inicial 0.1. Para c = 2(esquerda) a sequencia converge, mas para c = 4 (direita) ocomportamento e caotico. . . . . . . . . . . . . . . . . . . . . . 111

11.1 Evolucao do sistema de Fibonacci . . . . . . . . . . . . . . . . . 11911.2 Mapa de Henon. . . . . . . . . . . . . . . . . . . . . . . . . . . 12011.3 Triangulo de Sierpinski e quadrado fractal. . . . . . . . . . . . . 12211.4 Fractal em forma de folha, gerada com 4 transformacoes. . . . 124

vi LISTA DE FIGURAS

1. Introducao

Os sistemas dinamicos tem um papel bastante importante na engenha-ria. Os metodos utilizados na teoria dos sistemas dinamicos podem ser aplica-dos a uma grande variedade de problemas; alguns exemplos sao o crescimentodemografico de uma populacao, a interaccao entre duas especies biologicas di-ferentes, as orbitas planetarias, o sistema de amortecimento de um automovele um circuito electrico receptor de ondas de radio.

1.1 Sistemas dinamicos

Os sistemas dinamicos surgem em analogia com o problema central dadinamica. Nomeadamente, conhecidas as causas do movimento (forcas exter-nas) de um objecto, e conhecida a sua posicao e velocidade num instante inicial,calcular a sua posicao e velocidade (dois vectores ~r e ~v) em qualquer instanteposterior.

Em termos matematicos, as equacoes que permitem resolver esse pro-blema sao:

d~v

dt=

~F

m(1.1)

d~r

dt= ~v (1.2)

onde ~F e a forca externa resultante e m a massa do objecto. A resolucaodessas equacoes diferenciais permite encontrar expressoes para a posicao ~r e avelocidade ~v em funcao do tempo t.

Um sistema dinamico em geral sera caracterizado por um conjunto devariaveis de estado que dependem do tempo, e um conjunto de equacoesde evolucao. Esse conjunto de variaveis definem o estado do sistema num

Introducao 1

instante. Dado um estado inicial, as equacoes de evolucao permitem calcularo estado num instante posterior.

Os pormenores de cada sistema dinamico poderao ser muito diversos,mas podemos estudar algumas propriedades comuns a muitos sistemas dife-rentes. Por exemplo, na sua evolucao temporal os sistemas dinamicos podemapresentar alguns comportamentos tıpicos:

· O estado muda em forma periodica.

· O estado do sistema nao muda.

· O sistema apresenta uma tendencia a regressar a um determinado estado.

· O estado do sistema muda constantemente, passando por todos os esta-dos possıveis.

A figura 1.1 mostra um exemplo de um sistema com estado a mudarperiodicamente. As variaveis de estado sao a posicao vertical x da esfera,medida a partir da posicao de repouso da mola, e a velocidade vertical daesfera. O estado em cada instante e um ponto (x, v) num espaco de fase.O grafico da direita mostra a evolucao do estado para dois oscilacoes comamplitudes A1 e A2. A mola e esticada ate x = A1, e a esfera e largada apartir do repouso; o estado inicial e o ponto (A1, 0) no espaco de fase. Avelocidade comeca a aumentar em modulo, para cima, o qual corresponde auma velocidade negativa que diminui ate chegar a um valor mınimo em x = 0.Uma oscilacao completa corresponde a uma elipse no espaco de fase.

1.2 Sistemas Computacionais Algebricos (CAS)

Existem programas de computador, designados de CAS (Computer Al-gebra System), que permitem resolver equacoes do genero das equacoes 1.1 e1.2 em forma analıtica. Neste texto vamos fazer uso intensivo dum desses sis-temas e nao vamos estudar sistematicamente as tecnicas matematicas usadaspara resolver equacoes diferenciais ou equacoes de diferencas.

O uso dum sistema CAS nao elimina a necessidade de raciocınio ma-tematico por parte dos alunos nem torna o ensino puramente tecnico. Umadas complicacoes inerentes a computacao algebrica e o facto de nao existirem

2 Introducao aos Sistemas Dinamicos usando Maxima

x

0

x

v

A1 A2

Figura 1.1: Espaco de fase de uma esfera ligada a uma mola, mostrando duaspossıveis trajectorias correspondentes a duas oscilacoes.

respostas unicas. Diferentes ferramentas podem produzir respostas a um mes-mo problema que aparentemente sao diferentes, mas que representam a mesmafuncao. Ou as respostas podem ser funcoes de facto diferentes, que so coincidempara uns determinados valores dos parametros do problema. Inclusivamenteas mesma ferramenta CAS pode produzir duas respostas aparentemente dife-rentes se o problema for resolvido com metodos diferentes; em alguns casos osistema CAS pode nao dar nenhuma resposta ou a resposta dada pode estarerrada.

E preciso ganhar alguma experiencia para poder utilizar as ferramentasCAS de forma a obter as respostas na forma mais simples e poder confiarna sua veracidade. No processo de obter essa experiencia, o utilizador acabapor adquirir algum conhecimento sobre os metodos matematicos usados pelosistema CAS.

Hoje em dia a grande maioria dos profissionais de engenharia e cienciasexactas dependemos das maquinas de calcular para obter a raiz quadrada de umnumero, por exemplo 3456. Nao me parece que essa dependencia seja grave,nem que todos os profissionais tenham que aprender a calcular

√3456 com

papel e caneta antes de passar a usar a calculadora, enquanto existam aindalivros que expliquem metodos para calcular a raiz quadrada caso cheguemos aprecisar deles.

Por outro lado, quem possui uma calculadora que permite obter raızes

Introducao 3

quadradas pode avancar mais rapidamente para outros temas como o calculo deequacoes quadraticas e a traves desse estudo pode ate ganhar um bom conhe-cimento da funcao

√x, embora nao saiba como calcular raızes sem calculadora.

Os seus problemas matematicos deslocam-se da area do calculo de raızes parao problema da obtencao de solucoes para uma equacao quadratica. No casodas equacoes diferenciais e equacoes de diferencas, o uso de ferramentas CASpermite avancar mais facilmente para outros temas como o caos e as fractais.

1.3 Maxima

Maxima (http://maxima.sourceforge.net) e um dos sistemas CAS maisantigos. Foi criado pelo grupo MAC no MIT, na decada de 1960, e inicialmentechamava-se Macsyma (project MAC’s SYmbolic MAnipulator). Macsyma foidesenvolvido inicialmente para os computadores de grande escala DEC-PDP-10que eram usados em varias instituicoes academicas

Na decada de 1980 foi portado para varias novas plataformas, e uma dasnovas versoes foi designada de Maxima. Em 1982 o MIT decidiu comercializarMacsyma e, em forma paralela, o professor William Schelter da Universidade deTexas continuou a desenvolver Maxima. Na segunda metade da decada de 1980apareceram outros sistemas CAS proprietarios, como Maple e Mathematica, eem 1999 a versao proprietaria de Macsyma foi vendida e retirada do mercado,incapaz de concorrer contra os outros sistemas CAS propietarios. Em 1998 oprofessor Schelter obteve autorizacao do DOE (Department of Energy), quetinha os direitos de autor sobre a versao original de Macsyma, para distribuiro codigo fonte de Maxima.

Apos a morte do professor Schelter em 2001, formou-se um grupo de vo-luntarios que continuam a desenvolver e distribuir Maxima como software livre.Maxima partilha as muitas vantagens do software livre: acesso ao codigo fonte,criacao de um espırito de comunidade e partilha, criacao de patrimonio cultu-ral da humanidade, maior seguranca, preco mais baixo. No caso dos sistemasCAS algumas dessas vantagens sao ainda mais importantes. Como foi referidona seccao anterior, as respostas dos sistemas de CAS nao sao unicas nem per-feitamente determinadas; a implementacao de um metodo algebrico pode serfeita de muitas formas diferentes. Quando um metodo falha ou da respostasmuito complicadas e bastante util ter acesso aos pormenores da implementacao


http://maxima.sourceforge.net

subjacente ao sistema CAS. Por outra parte, no momento em que comecarmosa depender dos resultados dum sistema CAS, e desejavel que a documentacaodos metodos envolvidos esteja disponıvel e que nao existam impedimentos le-gais que nos proibam tentar descobrir ou modificar esses metodos.

1.4 Problemas

1. Tendo em conta a definicao de sistema dinamico dada neste capıtulo,proponha 3 exemplos de sistemas dinamicos. Em cada caso identifiqueas variaveis de estado.

2. Outros dois sistemas CAS com licenca de software livre sao Axiom eYacas. Procure informacao sobre Maxima, Axiom e Yacas; observe asfuncionalidades de cada um deles e o estado actual de cada projecto.Conhece algum outro sistema CAS?

Introducao 5

2. Computacao algebrica usando Maxima

2.1 A interface do Maxima

Existem varias interfaces diferentes para trabalhar com Maxima. Algu-mas delas podem estar incorporadas dentro de um editor de texto e a repre-sentacao das equacoes no ecran pode ter um aspecto grafico simplista ou maissofisticado em diferentes interfaces. Neste documento vamos admitir que se estaa usar a interface xmaxima, versao 5.9.0 (figura 2.1), que e a que normalmenteaparece por omissao e funciona igual em diferentes plataformas.

2.2 Entrada e saıda

Para lancar a interface, existira provavelmente no seu menu uma opcao“xmaxima” ou simplesmente “maxima”; desde um terminal tipo Unix podeusar o comando xmaxima. A interface de xmaxima reconhece-se por ter umabarra de menu com sub-menus “File”, “Edit”, “Options”, “Maxima” e “Help”(ver figura 2.1). Sugiro que antes de comecar a trabalhar entre no sub-menu“Options” e seleccione a opcao “Toggle Browser Visibility”; enquanto nao es-tiver a ler o tutorial ou o manual, pode aproveitar a area completa da janelapara introduzir comandos do Maxima. Se algum momento seleccionar “Ma-xima Help”, devera lembrar-se de fazer aparecer novamente o “browser” parapoder consultar a ajuda. Uma outra opcao que me parece muito util e “PlotWindows”; aconselho seleccionar “Separate” para que quando seja criado umgrafico apareca numa janela separada.

Quando se inicia uma sessao do maxima, aparece um sımbolo (C1) se-guido do prompt. Os sımbolos C1, C2, C3, . . . representam os comandos in-


Figura 2.1: Xmaxima, versao 5.9.0.

seridos pelo utilizador; e as respostas obtidas a cada comando representam-sepor D1, D2, D3, . . .. Comecemos por exemplo por fazer umas contas simples,2, 5× 3, 1 + 5, 2 log(2), guardando o resultado numa variavel res:

(C1) 2.5*3.1;

(D1)

7.75

(C2) 5.2*log(2);

(D2)

5.2 log 2

(C3) D2, numer;

Computacao algebrica usando Maxima 7

(D3)

3.604365338911716

(C4) res: d1 + D3;

(D4)

11.35436533891172

Podiamos ter feito tudo numa linha so sem ter que o fazer por partes.Repare que cada comando deve ser finalizado com “;”. Se nao for assim, naoaparecera o resultado (D1) pois Maxima pensara que estamos a escrever umcomando de varias linhas. Os numeros irracionais, como log(2) nao sao calcu-lados como um numero decimal, pois Maxima trabalha com a maior precisaonumerica possıvel. Para forcar log(2) a dar um numero decimal, usamos a“D2, numer” na alınea C3, que significa “o valor numerico do resultado D2”.Na alınea C4 repare que para dar um valor a uma variavel usa-se “:” e nao“=”, que sera utilizado para outros fins.

2.3 Variaveis

O nome das variaveis podera ser qualquer combinacao de letras, numerose os sımbolos % e . O primeiro caracter no nome da variavel nao pode serum numero. Maxima nao faz distincao entre maiusculas e minusculas. Algunsexemplos:

(C5) x1 : 2;

(D5)

2

(C6) area : 5$

(C7) %d_23 : 8;

(D7)

8


(C8) a%2 : (x1 : x1 + 2, x1*x1);

(D8)

16

Na alınea C6 usamos $ em vez de ponto e vırgula para terminar o co-mando. O $ no fim faz com que o comando seja executado, mas o resultadonao seja apresentado na interface. Na alınea C8 foi escrita uma sequencia decomandos, entre parentesis e separados por coma; os comandos sao executadosem sequencia, e o resultado e armazenado na variavel a%2; o primeiro comandona sequencia incrementa o valor de x1 em 2, e logo calcula o quadrado de x1.

Alguns nomes de variaveis nao podem ser usados por estarem reservados.Ja vimos que os nomes Cn e Dn, onde n e um inteiro positivo, estao reservadospara referir os comandos inseridos numa sessao, e os resultados obtidos. Umavariavel tambem nao pode ter o mesmo nome de algum comando do Maxima;por exemplo FOR, WHILE e SUM.

2.4 Constantes

Existem algumas constantes importantes ja predefinidas em Maxima.Os seus nomes comecam sempre por %. Tres constantes importantes saoo numero π, representado por %pi, o numero e, base dos logaritmos natu-rais, representado por %e, e o numero imaginario i =

√−1, representado por

%i.

(C9) %pi, numer;

(D9)

3.141592653589793

(C10) %e, numer;

(D10)

2.718281828459045


(C11) (3 + %i*4) * (2 + %i*5), expand;

(D11)

23i− 14

No ultimo comando calculamos o produto entre dois numeros complexos,utilizando o comando expand que obriga a que as expressoes entre parentesissejam multiplicadas.

2.5 Equacoes

O sinal de igualdade e usado para equacoes matematicas. Uma equacaomatematica pode ser memorizada numa variavel.

Exemplo 2.1Encontre as solucoes da equacao cubica 3x3 + 5x2 − x + 6 = 0.

(C12) eq1: 3*x^3 + 5*x^2 -x + 6 = 0;

(D12)

3x3 + 5x2 − x + 6 = 0

(C13) solve(eq1, x);

(D13)x =34(√

3i2 − 1

2

)81(√

1333154√

3− 1843

1458

) 13

+

(√13331

54√

3− 1843

1458

) 13(−√

3i

2− 1

2

)− 5

9,

x =

(√13331

54√

3− 1843

1458

) 13(√

3i

2− 1

2

)+

34(−√

3i2 − 1

2

)81(√

1333154√

3− 1843

1458

) 13− 5

9,

x =

(√13331

54√

3− 1843

1458

) 13

+34

81(√

1333154√

3− 1843

1458

) 13− 5

9


-3 -2 -1 1

-25

-20

-15

-10

-5

5

10

x

3x3 + 5x

2 - x + 6

Figura 2.2: Grafico do polinomio 3x3 + 5x2 − x + 6.

O comando solve foi usado para resolver a equacao para a variavel x.O resultado e uma lista com tres elementos; cada um desses 3 elementos e umaequacao que define uma das solucoes da equacao. Ha duas raızes complexas euma raiz real. Se quisermos isolar a raiz real, podemos usar o comando part;a raiz real e a terceira parte da lista anterior, e dentro dessa terceira partequeremos isolar o lado direito da equacao (segunda parte):

(C14) part(%, 3, 2), numer;

(D14)

−2.221834293486761

A figura 2.2 mostra o grafico do polinomio, obtido com os comandos:

(C15) pol: 3*x^3 - 5*x^2 - x + 6;

(D15)

3x3 − 5x2 − x + 6


(C16) "plot2d(pol, [x, -3, 1])"$

Exemplo 2.2Calcule a corrente em cada uma das resistencias do circuito.

4 kΩ 2 kΩ

1 kΩ

5 kΩ

4 V

8 kΩ6 V

Admitindo que I1 e I2 sao as correntes de malha, medidas em mili am-peres, em sentido dos ponteiros do relogio, as equacoes das duas malhas nocircuito sao

(C17) malha1: (4 + 8)*I1 - 8* I2 = 6 + 4;

(D17)

12I1− 8I2 = 10

(C18) malha2: (2+ 8 + 5 + 1)*I2 - 8*I1 = -4;

(D18)

16I2− 8I1 = −4

podemos usar mais uma vez o comando solve para resolver o sistemade duas equacoes; as duas equacoes e as duas variaveis a calcular deverao serincluıdas dentro de uma lista:

(C19) solve([malha1, malha2], [I1, I2]);

(D19) [[I1 = 1, I2 =

14

]]


A corrente na resistencia de 4 kΩ e I1 = 1 mA, a corrente nas resistenciasde 2 kΩ, 5 kΩ e 1 kΩ e I2 = 0, 15 mA, e a corrente na resistencia de 8 kΩ eigual a I1 − I2 = 0, 75 mA.

Observe que a resolucao do sistema de equacoes mostra os valores paraI1 e I2, mas nao guarda nenhum valor nas variaveis I1 e I2. O sistema anteriortambem poderia ter sido resolvido mais rapidamente com linsolve, em vez desolve, por tratar-se de um sistema de equacoes lineares.

2.6 Funcoes

Para definir funcoes usa-se o sımbolo :=

(C20) f(x) := 3 + x^2;

(D20)

f (x) := 3 + x2

(C21) f(5);

(D21)

28

(C22) g(x,y,z) := x*y^2 + z;

(D22)

g (x, y, z) := xy2 + z

(C23) g(2,3,4);

(D23)

22

Tambem existem funcoes que definem listas. Por exemplo

(C24) cubo[n] := n^3;

(D24)

cubon := n3


(C25) makelist(cubo[i], i, 1, 8);

(D25)

[1, 8, 27, 64, 125, 216, 343, 512]

o comando makelist foi usado para criar uma lista com os oito primeiroscubos. Ha que ter algum cuidado com as funcoes que definem listas, poisquando um determinado elemento da lista ja tem sido calculado, a funcaonao voltara a ser usada para o calcular. Por exemplo, se agora redefinirmoscubo

(C26) cubo[n] := 4*n^3;

(D26)

cubon := 4n3

e pedissemos o valor de cubo[3] continuara a ser 27, e nao 108, ja quequando criamos a lista dos oito primeiros cubos ja foi calculado cubo, ficandocom o valor 27. Se quisermos redefinir uma lista, sera melhor primeiro apaga-lausando kill

(C27) cubo[n] := 4*n^3;

(D27)

cubon := 4n3

(C28) cubo[3];

(D28)

27

(C29) kill(cubo);

(D29)

DONE

(C30) cubo[n] := 4*n^3;

(D30)

cubon := 4n3


(C31) cubo[3];

(D31)

108

Exemplo 2.3O banco empresta-nos ¤500, a uma taxa de juros de 5% anual, com prazo de20 meses. Qual sera a prestacao mensal que temos que pagar?

A taxa de juro mensal sera

(C32) J : 0.05/12$

Se M[n] for o montante em dıvida no mes numero n, este sera igual a omontante que ja tınhamos em dıvida no mes anterior, mais os juros devidos nes-te perıodo, menos o a prestacao (P) paga neste mes:

(C33) M[n] := expand(M[n-1] + J*M[n-1] - P);

(D33)

Mn := EXPAND(Mn−1 + JMn−1 − P )

A funcao expand foi usada porque como P e uma variavel desconheci-da, Maxima deixa varios produtos indicados, tornando a expressao para M[n]

complicada, a menos que os produtos em M[n-1] ja tenham sido calculados.O nosso problema consiste em encontrar o valor de P que faz com que

M[20] seja nula, comecando com M[0]= 500, que e o valor inicialmente emdıvida.

(C34) M[0] : 500$

(C35) sol: solve(M[20] = 0, P);

RAT replaced 543.3579448601711 by 95631//176 = 543.3579545454545

RAT replaced -20.8118135328821 by -15151//728 = -20.8118131868132

(D35) [P =

8702421333322

]


(C36) part(sol, 1, 2), numer;

(D36)

26.10815067712302

A prestacao mensal sera de ¤26,11.

2.7 Algebra e trigonometria

Maxima facilita a manipulacao de expressoes algebricas. Por exemplo,vamos expandir um polinomio:

(C37) (x + 4*x^2*y + 2*y^2)^3;

(D37) (2y2 + 4x2y + x

)3(C38) expand(%);

(D38)

8y6 +48x2y5 +96x4y4 +12xy4 +64x6y3 +48x3y3 +48x5y2 +6x2y2 +12x4y+x3

Se quisermos agora substituir x por 1/(z − 3), escrevemos

(C39) %, x=1/(z-3);

(D39)

12y4

z − 3+

48y5

(z − 3)2+

6y2

(z − 3)2+

48y3

(z − 3)3+

1(z − 3)3

+96y4

(z − 3)4+

12y

(z − 3)4+

48y2

(z − 3)5+

64y3

(z − 3)6+ 8y6


e para reduzir tudo a um denominador comum usamos a funcao ratsimp

(o resultado ocupa varias linhas e nao vamos apresenta-lo)

(C40) ratsimp(%)$

Podemos factorizar a expressao anterior, para obter um resultado maissimples

(C41) factor(%);

(D41) (2y2z2 − 12y2z + z + 18y2 + 4y − 3

)3(z − 3)6

A funcao trigexpand serve para expandir senos ou co-senos de somas oudiferencas de angulos; trigreduc tenta expandir de forma que cada termo sotenha uma funcao seno ou co-seno.

(C42) sin(u+v)*cos(u)^3;

(D42)

cos3 u sin (v + u)

(C43) trigexpand(%);

(D43)

cos3 u (cos u sin v + sinu cos v)

(C44) trigreduce(%);

(D44)

sin (v + 4u) + sin (v − 2u)8

+3 sin (v + 2u) + 3 sin v

8

2.8 Passar informacao entre sessoes

Para guardar o conteudo de uma sessao em Xmaxima, existe a opcao “Sa-ve Console to File” no menu “Edit”. Essa opcao guarda toda a informacao que


apareceu no ecran, incluindo os sımbolos (Cn) e (Dn). Para gravar os comandosexecutados, numa forma que possa ser aproveitada em sessoes posteriores, usa-se o comando stringout. Vejamos um exemplo1

(C45) stringout("/home/villate/trig.mac", C41, D42)$

(C46) stringout("/home/villate/prestacao.mac", [31, 35])$

(C47) stringout("/home/villate/capitulo2.mac", input)$

No ficheiro /home/villate/trig.mac fica armazenado:

SIN(u+v)*COS(u)^3;

COS(u)^3*(COS(u)*SIN(v)+SIN(u)*COS(v));

No ficheiro /home/villate/prestacao.mac ficam guardados os passosnecessarios (C31, C32, . . ., C35) para resolver o Exemplo 2.3 acima:

J:0.05/12;

M[n]:=EXPAND(M[n-1]+J*M[n-1]-P);

M[0]:500;

sol:SOLVE(M[20] = 0,P);

EV(PART(sol,1,2),NUMER);

esses passos podem ser executados posteriormente por meio do comando

(C48) batch("/home/villate/prestacao.mac")$

Finalmente, o ficheiro /home/villate/capitulo2.mac tera uma copiade todos os comandos usados neste capıtulo. O conteudo desse ficheiro e textosimples, que pode ser modificado com um editor de texto e executado emMaxima com o comando batch.

1Em Windows sera preciso usar algo como C:\\MeusDocumentos\\trig.mac para os nomesdos ficheiros (com barras a dobrar).


2.9 Problemas

1. A equacao quadratica x2 = x + 1 aparece frequentemente em problemasde dinamica de sistemas. Essa equacao tem duas raızes reais, positiva enegativa. Usando Maxima, demonstre que a raiz positiva e o numero

φ =1 +

√5

2

(em Maxima esse numero e a constante %phi) e a raiz negativa e −1/φ.Use Maxima para demonstrar as propriedades:

φn+2 = φn+1 + φn 1φ

= φ− 1

O numero φ, designado de relacao aurea, e conhecido desde a epocados gregos, e era considerada a relacao geometrica perfeita para obteros melhores resultados esteticos.

2. Encontre a equacao da circunferencia que passa pelos pontos (−2, 7),(−4, 1) e (4,−5).

Sugestao: Procure uma solucao da forma:

(x− a)2 + (y − b)2 = r2

3. Desenhe o grafico de cada uma das seguintes funcoes, usando intervalosque mostrem bem a forma das funcoes.

y = x3 − 5x2 + 2x + 3, y =sin(x)

x. y =

√20− x2, y =

3x2 + 2x2 − 4

4. Defina um procedimento em maxima para calcular os numeros de Fi-bonacci, Fn = 1, 1, 2, 3, 5, 8, . . ., definidos, para (n = 0, 1, 2, 3, . . .), por:

F0 = 1 F1 = 1 Fn = Fn−1 + Fn−2

Calcule a relacao Fn+1/Fn para alguns valores crescentes de n, e mostreque a relacao aproxima-se do limite φ, onde φ e a relacao aurea (problema1). (A sucessao de Fibonacci foi proposta por Leonardo de Pisa, tambemconhecido por Filius Bonacci, no seculo XIII, para estudar o crescimentoanual de uma populacao de coelhos).


5. O grafico da funcao y = x3−6x2+7x+2 apresenta dois pontos extremos(designados de mınimo local e maximo local). Desenhe o grafico dessafuncao e faca uma estimativa das coordenadas x e y desses dois pontos.Melhore a sua estimativa usando o comando subst (consulte o sistemade ajuda em linha) para encontrar os valores maior ou menor de y nessasvizinhancas.

Sugestao: E util primeiro definir uma expressao f: x^3 - 6*x^2 + 7*x

+ 2. Faca varias tentativas ate conseguir melhorar o seu resultado.

Pode sugerir algum outro metodo grafico de obter uma melhor estimativainicial para os valores do maximo e o mınimo local?

6. A nota final nesta disciplina e calculada udando a formula

0.3 Freq + Trab (1− 0.015 Freq)

onde Freq e Trab sao as notas do trabalho feito nas aulas e o trabalhofeito em casa, ambas entre 0 e 20 valores. Demonstre que a nota finalsera sempre maior o igual do que a nota do trabalho de casa.


3. Alguns conceitos de calculo

3.1 Derivada de uma funcao

Para calcular a derivada de uma funcao, usa-se o comando diff. Oprimeiro argumento devera ser uma funcao de uma ou mais variaveis, o segun-do argumento e a variavel em ordem a que vai ser derivada a funcao, e umterceiro argumento optativo, que indica a ordem da derivacao (se nao apare-cer entender-se-a que se trata de uma derivada de primeria ordem). Algunsexemplos:

(C1) diff(x^n, x);

(D1)

nxn−1

Podemos usar o apostrofo para que maxima nos mostre, em notacao ma-tematica, o calculo que esta a ser feito:

(C2) ’diff(exp(a*x), x) = diff(exp(a*x), x);

(D2)deax

dx= aeax

(C3) ’diff(sin(x), x, 2) = diff(sin(x), x, 2);

(D3)d2 sinx

dx2= − sinx

Se a funcao a ser derivada depende de varias variaveis, pode escrever-seuma lista de variaveis de derivacao, seguida cada uma da sua ordem (nesse caso

Alguns conceitos de calculo 21

nao pode ser omitida a ordem de derivacao), para calcular derivadas parciais.Por exemplo:

(C4) ’diff(x^3/y, y, 1, x, 2) = diff(x^3/y, y, 1, x, 2);

(D4)∂3(x3/y

)∂2x∂y

= −6x

y2

3.1.1 Interpretacao geometrica da derivada

Em termos geometricos, a derivada de uma funcao e uma outra funcaoque em cada ponto e igual ao declıve da funcao original nesse mesmo ponto.Por exemplo, a figura 3.1 mostra o grafico da funcao log(x) e a sua derivada,1/x, obtido com o comando:

(C5) plot2d([log(x), 1/x], [x, 0.1, 10], [y, -1, 4])$

Perto da origem, o declıve de log(x) tem um valor elevado, mas paravalores maiores de x o declıve aproxima-se de zero.

2 4 6 8

-1

1

2

3

x

f(x)

log(x)1/x

Figura 3.1: Grafico de log(x) e a sua derivada, 1/x.


3.1.2 Algumas propriedades

Se u e v forem duas funcoes que dependem de x, podemos usar o comandodepends para definir essa dependencia

(C6) depends([u, v], x);

(D6)

[u (x) , v (x)]

e a seguir poderemos calcular propriedades gerais para quaisquer funcoesu e v. Por exemplo, a derivada do produto e do quociente das funcoes sao1:

(C7) ’diff(u*v, x) = diff(u*v, x);

(D7)d (uv)

dx= u

dv

dx+

du

dxv

(C8) ’diff(u/v, x) = ratsimp(diff(u/v, x));

(D8)d (u/v)

dx= −

u dvdx −

dudxv

v2

A regra da cadeia usa-se quando queremos calcular a derivada de umafuncao composta, isto e, uma funcao y que depende de uma outra funcao x,que pela sua vez depende da variavel t:

(C9) depends(y, x, x, t);

(D9)

[y (x) , x (t)]

Assim, y depende implicitamente da variavel t e a sua derivada em funcaode t sera:

(C10) ’diff(y, t) = diff(y, t);

(D10)dy

dt=

dx

dt

dy

dx1A funcao ratsimp usa-se para agrupar as fraccoes sobre um denominador comum.


Exemplo 3.1Constroi-se uma caixa de altura x, a partir de uma lamina quadrada, comaresta d = 30 cm, cortando quatro quadrados de aresta x nos cantos da lamina.Encontre o valor que devera ter x para que o volume da caixa seja o maximopossıvel.

d

d

x

x

O volume da caixa, em funcao de x, sera

(C11) vol(x) := x*(30 - 2*x)^2;

(D11)

vol (x) := x (30− 2x)2

Podemos fazer uma tabela de valores para ter alguma ideia do compor-tamento da funcao vol(x). Primeiro criamos uma lista de possıveis valores dex, que pode estar compreendido entre 0 e 15

(C12) xi : makelist(i, i, 0, 15);

(D12)

[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]

A seguir, vamos usar o comando map para aplicar a funcao vol a ca-da elemento da lista xi, e com a lista de valores de x e do volume cria-mos uma matriz onde cada linha representa um par de valores (x, vol(x))


(C13) yi : map(vol, xi);

(D13)

[0, 784, 1352, 1728, 1936, 2000, 1944, 1792, 1568, 1296, 1000, 704, 432, 208, 56, 0]

(C14) transpose(matrix(xi, yi));

(D14)

0 01 7842 13523 17284 19365 20006 19447 17928 15689 129610 100011 70412 43213 20814 5615 0

o comando matrix aceita uma ou varias listas com o mesmo numerode elementos e cria uma matriz onde cada lista e uma linha. O comandotranspose foi usado para que as linhas da matriz passem a ser colunas.

Aparentemente o valor maximo obtem-se para x = 5. Podemos demons-trar em forma mais precisa que essa e a resposta do problema, se calcularmosos pontos onde a derivada da funcao vol e nula

(C15) solve(diff(vol(x), x) = 0, x);

(D15)

[x = 5, x = 15]

Para determinar quais desses pontos sao maximos ou mınimos locais,calcula-se a segunda derivada


(C16) diff2 : diff(vol(x), x, 2);

(D16)

8x− 8 (30− 2x)

no ponto x = 5 o valor da segunda derivada e

(C17) diff2, x=5;

(D17)

−120

o seu valor negativo indica que a curvatura da funcao nesse ponto apontapara baixo e, portanto, a funcao tem um maximo local em x = 5.

Exemplo 3.2A posicao de uma partıcula, em funcao do tempo t, e definida pelo vector

~r = 5 cos(t)~ux + 4 sin(t)~uy + 2e−t~uz

calcule a velocidade e a aceleracao da partıcula em t = 3.

Vamos representar as tres coordenadas do vector posicao por meio deuma lista

(C18) pos: [5*cos(t), 4*sin(t), 2*exp(-t)];

(D18) [5 cos t, 4 sin t, 2e−t

]

A velocidade e a derivada do vector posicao, em funcao do tempo, e a ace-leracao e a derivada da velocidade.

(C19) vel: map(lambda([u], diff(u, t)), pos);

(D19) [−5 sin t, 4 cos t,−2e−t

](C20) acel: map(lambda([u], diff(u, t)), vel);

(D20) [−5 cos t,−4 sin t, 2e−t

]


nos comandos anteriores definimos uma funcao anonima, por meio docomando lambda, que, dada uma funcao u, calcula a sua derivada em ordema t. A seguir aplicamos essa funcao a cada uma das coordenadas do vectorposicao e do vector velocidade, usando o comando map do maxima. Em t = 3a velocidade e a aceleracao sao

(C21) vf : ev(vel, t=3, numer);

(D21)

[−0.70560004029934,−3.959969986401782,−0.09957413673573]

(C22) af : ev(acel, t=3, numer);

(D22)

[4.949962483002227,−0.56448003223947, 0.09957413673573]

Os modulos da velocidade e da aceleracao podem ser calculados, usando oproduto interno entre vectores, representado por um ponto:

(C23) modvel : sqrt(vf.vf);

(D23)

4.023574122441392

(C24) modacel: sqrt(af.af);

(D24)

4.983039363544433

3.2 Primitivas e integrais

A primitiva de uma funcao e qualquer outra funcao que quando derivadada como resultado a funcao original. Por exemplo, vamos encontrar as primiti-vas da funcao xn, onde n e um parametro constante


(C25) ’integrate(x^n, x) = integrate(x^n, x);

Is n+1 zero or nonzero?

nonzero;

(D25) ∫xn dx =

xn+1

n + 1

Maxima perguntou-nos se n + 1 e nula, isto e, se n e igual a −1. Paran diferente de −1, obtemos o resultado acima, que representa apenas uma dasprimitivas da funcao xn. Todas as outras primitivas obtem-se somando umaconstante arbitraria ao resultado anterior.

Um integral definido, por exemplo

(C26) ’integrate(1/(1 + x^ 2), x, 0, 1);

(D26) ∫ 1

0

1x2 + 1

dx

representa a area entre a funcao e o eixo dos x, e entre os limites x = 0 ex = 1 (ver a figura 3.2). Neste caso o valor do integral e:

(C27) ’integrate(1/(1 + x^ 2), x, 0, 1) =

integrate(1/(1 + x^ 2), x, 0, 1);

(D27) ∫ 1

0

1x2 + 1

dx =π

4

O valor infinito, representa-se em maxima pela variavel inf. Por exem-plo a funcao acima tambem pode ser integrada desde menos infinito ate infini-to:

(C28) ’integrate(1/(1 + x^ 2), x, -inf, inf) =

integrate(1/(1 + x^ 2), x, -inf, inf);

(D28) ∫ ∞

−∞

1x2 + 1

dx = π


0 x

1/(1 + x2)

0-1 1

0.5

1.0

Figura 3.2: Integral de 1/(1 + x2) entre 0 e 1.

Exemplo 3.3Calcule a area da regiao delimitada pela parabola y = 2x2 − 3 e a recta y = 2x

Primeiro temos que encontrar os pontos de interseccao entre a parabolae a recta

(C29) parab(x) := 2*x^2 - 3$

(C30) recta(x) := 2*x$

(C31) solve(parab(x) = recta(x), x);

(D31) [x = −

√7− 12

, x =√

7 + 12

]

A figura 3.3 mostra o grafico das funcoes, obtido com:

(C32) plot2d([parab(x), recta(x)], [x, -2, 3])$

Entre os dois pontos de interseccao, a recta esta por cima da parabola.A area sera


x

f(x)

-1 1 2

4

8

12

Figura 3.3: Area entre y = 2x2 − 3 e y = 2x.

(C33) area : integrate(recta(x) - parab(x), x,

part(d31,1,2), part(d31,2,2));

(D33)7√

7 + 106

+7√

7− 106

(C34) area, numer;

(D34)

6.173419725817378

3.3 Serie de Taylor

Uma funcao f(x) diz-se analıtica num ponto a, se a funcao e todas assuas derivadas, f(x), f ′(x), f ′′(x), f ′′′(x), . . . existirem no ponto a. Nesse caso


a funcao pode ser escrita como uma serie de potencias

f(x) = f(a) + f ′(a)(x− a) +f ′′(a)

2!(x− a)2 +

f ′′′(a)3!

(a− x)3 + . . . (3.1)

esse tipo de serie e designada de serie de Taylor. Em maxima, o comandopowerseries da a forma geral da serie de Taylor de uma funcao analıtica e ocomando taylor pode ser usado para obter os primeiros termos da expansao deuma funcao em serie de Taylor. 2

(C35) ’sin(x) = niceindices(powerseries(sin(x), x, 0));

(D35)

sinx =∞∑

i=0

(−1)ix2i+1

(2i + 1)!

(C36) ’cos(x) = niceindices(powerseries(cos(x), x, 0));

(D36)

cos x =∞∑

i=0

(−1)ix2i

(2i)!

Os quatro primeiros termos da serie binomial sao

(C37) taylor((1+x)^n, x, 0, 3);

(D37)

1 + nx +

(n2 − n

)x2

2+

(n3 − 3n2 + 2n

)x3

6+ · · ·

Factorizando cada termo e reagrupando como uma serie infinita (pormeio do comando trunc) obtemos

(C38) map(factor,%);

(D38)(n− 2) (n− 1) nx3

6+

(n− 1) nx2

2+ nx + 1

(C39) (1 + x)^n = trunc(%);

(D39)

(x + 1)n = 1 + nx +(n− 1) nx2

2+

(n− 2) (n− 1) nx3

6+ · · ·

2O comando niceindices substitui os ındices I1, I2, . . ., usados por omissao, por ındicesi, j, . . .


Uma serie de potencias pode ser derivada ou integrada termo a termo.Por exemplo, o integral que calculamos em D21 pode ser calculado por separadopara cada termo da serie de potencias da funcao e obtemos, assim, uma serie quepermite calcular π/4

(C40) serie : niceindices(integrate(

powerseries(1/(1+x^2),x,0),x));

(D40)∞∑

i=0

(−1)ix2i+1

2i + 1

(C41) integrate(1/(1 + x^2), x, 0, 1) = ev(serie, x=1) -

ev(serie, x=0);

(D41)π

4=

∞∑i=0

(−1)i

2i + 1

3.4 A transformada de Laplace

A transformada de Laplace de uma funcao f(x) e uma outra funcao f(s)definida por

f(s) =∫ ∞

0

f(x)e−sx dx (3.2)

por exemplo, a transformada da funcao seno e:

(C42) laplace(sin(x), x, s);

(D42)1

s2 + 1

Um outro exemplo:

(C43) laplace(t^2*exp(3*t)*sin(4*t), t, s);


(D43)8 (2s− 6)2

(s2 − 6s + 25)3− 8

(s2 − 6s + 25)2

A transformada inversa calcula-se usando o comando ilt. Por exem-plo:

(C44) ilt(1/(s^2 + 1), s, t);

(D44)

sin t

3.5 Problemas

1. Calcule as derivadas das funcoes seguintes:

(a) y = e2x cos(3x) sin(x)

(b) y =√

3x sinh(2/x)

(c) y =x3 + 4x2

x2 − 1

2. Calcule os pontos estacionarios (onde a derivada e nula) da funcao

3x4 − 4x3 − 24x2 + 48x + 5

e diga quais desses pontos sao maximos ou mınimos locais.

3. Calcule a primitiva de (log x)20 e comprove o resultado por meio dederivacao.

4. Calcule os seguintes integrais definidos:

(a)∫ 2

0

(4x4 − 5x2 + 8

)dx

(b)∫ π/3

0

sin(3x) dx


(c)∫ ∞

−∞e−x2

dx

5. Calcule a derivada da funcao f(x), onde:

(a) f(x) =∫ x5

1

√t2 + 1 dt

(b) f(x) =∫ sin(x)

1

(t2 + 1)3 dt

6. Em cada caso calcule a serie de Taylor da funcao, no ponto a, e ate otermo de ordem n.

(a) f(x) =1

(x− 4)2a = 5 n = 5

(b) f(x) = tan(x) a = π n = 6

(c) f(x) = x3/2 a = 1 n = 8

7. Calcule as transformadas de Laplace das seguintes funcoes:

(a) (1 + et)2

(b) sin(2t) sinh(2t)

(c) t2 + e2t − 2

(d) t2et sin(t)


4. Equacoes diferenciais ordinarias

Neste capıtulo, em algumas ocasioes usaremos a notacao simplificaday′, y′′,. . . para representar as derivadas de uma funcao y. Devera ficar claro,dentro do contexto do problema, qual a variavel em relacao a que se deriva afuncao. Por exemplo, se y for uma funcao da variavel x, y′ e y′′ representam:

y′ =dy

dxy′′ =

d2y

dx2(4.1)

Uma equacao diferencial ordinaria —ou em forma curta, uma EDO—e uma equacao que relaciona uma funcao y com as suas derivadas y′, y′′, . . .

e com a variavel x. Por exemplo, 3y′′ − (x − 2)y = y′. As funcoes y(x) queverifiquem a EDO serao as solucoes da equacao.

Em muitos problemas, quando andamos a procura de uma funcao quedefina o valor de uma determinada grandeza fısica, as leis fısicas do sistemaconduzem a equacoes diferenciais que deverao ser resolvidas para encontrar afuncao que procuramos.

4.1 Equacoes de primeira ordem

A ordem de uma EDO e a ordem da derivada de ordem superior queapareca na equacao. Se a unica derivada da funcao y na equacao e a suaprimeira derivada y′, a equacao e de primeira ordem.

A forma geral das equacoes de primeira ordem e1

dy

dx= f(x, y) (4.2)

1O problema de escrever uma equacao como (3y′)2 = x sin(y′) na forma geral 4.2 e umproblema algebrico que nao vamos abordar aqui. Em alguns casos uma EDO pode conduzira um sistema de varias equacoes com a forma geral 4.2.

Equacoes diferenciais ordinarias 35

Esta equacao pode ser escrita numa outra forma equivalente, designadade equacao inversa

dx

dy=

1f(x, y)

(4.3)

Na equacao inversa y passa a ser a variavel independente e x a variaveldependente; no entanto, as solucoes implıcitas das duas equacoes sao as mes-mas. Por exemplo, a solucao de uma equacao pode ser y = x2; esta solucaoe uma funcao de x, na forma explıcita y = f(x), mas se quisermos escrever x

em funcao de y, a mesma solucao conduz a duas solucoes explıcitas x =√

y ex = −√y.

Para escrever uma EDO em Maxima pode usar-se o apostrofo para que asderivadas fiquem indicadas, sem serem calculadas; por exemplo:

(C1) (2*x^2)/(y^3)*’diff(y,x) = (3+y)*sin(x);

(D1)

2x2

y3

dy

dx= sinx (y + 3)

(C2) 2*x + 2*y*’diff(y,x) = 0;

(D2)

2ydy

dx+ 2x = 0

Uma outra forma de escrever as equacoes diferenciais em Maxima consis-te em usar o comando depends para declarar que y e uma variavel dependentede x:

(C3) depends(y, x);

(D3)

[y (x)]

(C4) (2*y + exp(x)*cos(y))*diff(y, x) = -exp(x)*sin(y);

(D4)

(ex cos y + 2y)dy

dx= −ex sin y


4.2 Campo de direccoes

Na equacao 4.2, a funcao f(x, y) define, em cada ponto do plano (x, y),o declıve que devera ter uma curva y(x) que seja solucao da EDO. O campode direccoes e um desenho do plano (x, y), onde em cada ponto aparece umvector com declıve igual a f(x, y).

A figura 4.1 mostra o campo de direccoes da equacao diferencial y′ =x2 + 3y e algumas solucoes da equacao. As solucoes sao as curvas que em cadaponto seguem o sentido do campo de direccoes.

1

1

-2.5

-2.5

0

0

2.5

2.5

-2.5 -2.5

0 0

2.5 2.5

Figura 4.1: Campo de direccoes de y′ = x2 + 3y, e seis solucoes particulares.

Para desenhar um campo de direccoes em Maxima, e preciso primei-ro definir uma funcao que chamaremos plotdf. Entra-se numa sessao doLisp, por exemplo carregando simultaneamente em Ctrl e g, ou seleccionan-do “Interrupt”, no menu “File” (devera aparecer o prompt MAXIMA>>). A


seguir, define-se a funcao assim2:

(defun $plotdf() (show-open-plot "plotdf"))

Regressa-se a consola do Maxima (com :q na versao GCL do Maxima, ou com:rt seguido de 0, na versao Clisp) e executa-se a funcao que se acabou de definir:

(C5) plotdf();

Aparecera uma janela com um campo de direccoes predefinido. Passandoo cursor sobre os numeros no canto superior esquerdo, aparece o menu; a opcao“Config” no menu permite escrever a EDO que se quer representar. Selecciona-se “dy/dx” e escreve-se a funcao f(x, y) no campo dy/dx. Depois de seleccionarOK, sera preciso usar a opcao “Replot” no menu, para que seja desenhado ocampo da funcao que acaba de ser escrita. Carregando com o rato sobre umponto, aparece desenhada a solucao y(x) que passa por esse ponto.

4.3 Problemas com condicao inicial

Uma EDO tem infinitas solucaoes. Por exemplo, a figura 4.1 mostra seispossıveis solucoes da equacao y′ = x2 +3y. Para determinar uma solucao unicade uma equacao de primeira ordem

dy

dx= f(x, y) (4.4)

e preciso impor uma condicao inicial

y = y0 em x = x0 (4.5)

No desenho do campo de direccoes, a condicao inicial corresponde a um ponto(x0, y0), e a solucao do problema sera a curva y(x) que passa por esse ponto esegue as direccoes do campo de direccoes.

2Esta definicao pode tambem ser copiada num ficheiro, em algum dos directorios onde o Ma-xima procura funcoes externas. Por exemplo, $home/.maxima/plotdf.lisp em GNU/Linux,e podera ser recuperada em sessoes posteriores com o comando load("plotdf")


4.4 Resolucao de equacoes de primeira ordem

4.4.1 Um caso simples

Se a funcao no lado direito da equacao 4.2 depende unicamente de x,isto e,

dy

dx= f(x) (4.6)

a equacao representa simplesmente a definicao da derivada de uma funcao y(x).A solucao obtem-se facilmente por primitivacao

y =∫

f(x) dx + C (4.7)

onde C e podera ser qualquer constante arbitraria. Cada possıvel valor deC define uma solucao particular; a famılia de funcoes definida por 4.7 edesignada de solucao geral.

4.4.2 Equacoes autonomas

Se a funcao no lado direito da EDO depende unicamente de y,

dy

dx= f(y) (4.8)

designamos a EDO de equacao autonoma. Neste caso a equacao inversa e dotipo da equacao 4.6 na seccao anterior e resolve-se facilmente por primitivacao.

x =∫

1f(y)

dy + C (4.9)

4.4.3 Outras equacoes

No caso geral, pode usar-se o comando ode2 de Maxima para encontrara solucao geral duma equacao diferencial, e para aplicar a condicao inicial asolucao geral obtida, usa-se o comando ic1

Exemplo 4.1A partir dos dados demograficos para Portugal em Julho de 1993:

populacao: 10 486 140taxa de natalidade anual: 11,59 por mil hab.taxa de mortalidade anual: 9,77 por mil hab.taxa de migracao anual: 1,8 por mil hab.

Faca uma estimativa da populacao em 2010, admitindo taxas constantes (mo-


delo de Malthus) e taxa de mortalidade directamente proporcional a populacao(modelo logıstico).

Se t for o tempo, em anos, e p(t) o numero de habitantes num determina-do instante, p′ sera a taxa de aumento anual da populacao. A taxa de aumentoda populacao sera igual a p vezes a taxa de aumento por cada habitante. Comos dados do problema, a taxa de aumento por um habitante e

(C6) r: (11.59 - 9.77 + 1.8)/1000$

A taxa de aumento da populacao e, por um lado, rp, e por outro lado,a derivada p′. Assim, obtemos a equacao diferencial

(C7) eq1: ’diff(p,t) = r*p;

(D7)dp

dt= 0.00362p

a solucao geral obtem-se com o comando ode2

(C8) sol1: ode2(eq1, p, t);

(D8)

p = %Ce181t50000

Se o tempo t for medido a partir de 1993, a condicao inicial escreve-se

(C9) malthus: ic1(sol1, t=0, p=10486140);

(D9)

p = 10486140e181t50000

A populacao em 2100 sera,

(C10) malthus, t=17, numer;

(D10)

p = 1.1151727127034325 · 107

Segundo este modelo, a populacao cresce indefinidamente


(C11) limit(malthus, t, inf);

(D11)

p = ∞

Um modelo mais realista, designado de modelo logıstico, admite quea taxa de mortalidade, por cada habitante, aumenta em forma directamenteproporcional ao numero de habitantes. Isto e,

taxa mortalidade = k × p

onde k e uma constante; usando os dados do problema,

(C12) k: (9.77/1000)/10486140$

As taxas de nascimento e migracao, por habitante, admitem-se constantes:

(C13) a: (11.59 + 1.8)/1000$

Assim, a equacao diferencial para o aumento da populacao e

(C14) eq2: ’diff(p, t) = (a - k*p)*p;

(D14)

dp

dt=(0.01339− 9.3170604245222741 · 10−10p

)p

e a resolucao do problema e:

(C15) sol2: ode2(eq2, p, t)$

(C16) logistic: ic1(sol2, t=0, p=10486140);

(D16)

−39059 log (39059p− 561335846131)− 39059 log p

523=

523t + 39059 log 10486140− 39059 log (−151757703871)523


10 50

5

0

0

25

25

50

50

75

75

0 0

10 10

20 20

Tempo, em anos a partir de 1993

Pop

ulaç

ão, e

m m

ilhõe

s

Figura 4.2: Aumento da populacao, no modelo logıstico.

Aqui vamos dar alguma ajuda ao Maxima para resolver a equacao para p:

(C17) plogistic: solve(p/(39059*p-561335846131)=

exp(rhs(logistic*523/39059)),p);

(D17) [p =

5886246269548124340e523t39059

409578142260e523t39059 + 151757703871

]

(C18) plogistic[1], t=17, numer;

(D18)

p = 1.1096894862054735 · 107


Neste caso a populacao nao cresce indefinidamente, mas alcanca um limi-te de aproximadamente 14 milhoes de habitantes:

(C19) limit(plogistic[1], t, inf), numer;

(D19)

p = 1.4371485346040605 · 107

A figura 4.2 mostra o mapa de direccoes e a solucao particular.

4.5 Classificacao das equacoes de primeira ordem

O comando ode2 do Maxima classifica a EDO dada e se estiver dentrode alguma classe conhecida, sera resolvida usando o metodo conhecido pararesolver esse tipo de equacao. Para saber o tipo de metodo utilizado (e assim otipo de equacao), pede-se o valor da variavel method, logo a seguir ao resultadode ode2. A continuacao vamos ver alguns tipos de equacoes de primeira ordem.

4.5.1 Equacoes de variaveis separaveis

Sao equacoes que podem ser escritas na forma

dy

dx=

f(x)g(y)

(4.10)

a solucao geral e ∫g(y) dy =

∫f(x) dx + C (4.11)

Exemplo 4.2Encontre a solucao geral da equacao

y′ sin(x) = 4y3

(C20) ode2(’diff(y,x)*sin(x) = 4*y^3, y, x);

(D20)

− 18y2

= %C− log (cos x + 1)− log (cos x− 1)2


(C21) solve(%, y);

(D21) [y = − 1

2√

log (cos x + 1)− log (cos x− 1)− 2%C,

y =1

2√

log (cos x + 1)− log (cos x− 1)− 2%C

]

(C22) method;

(D22)

SEPARABLE

4.5.2 Equacoes exactas

Uma EDO de primeira ordem pode sempre ser escrita na forma

M(x, y) + N(x, y)dy

dx= 0 (4.12)

que e semelhante a expressao obtida para a derivada de uma funcao F (x, y),de duas variaveis (x, y), com y dependente de x:

dF

dx=

∂F

∂x+

∂F

∂y

dy

dx(4.13)

se conseguirmos encontrar uma funcao F (x, y) com derivadas parciais

∂F

∂x= M(x, y)

∂F

∂y= N(x, y) (4.14)

entao a equacao diferencial 4.12 diz simplesmente que F (x, y) e uma funcaoconstante, nomeadamente, a solucao geral de 4.12 sera:

F (x, y) = C (4.15)

A condicao necessaria e suficiente para que exista a funcao F que verificaas equacoes 4.14 e

∂M

∂y=

∂N

∂x(4.16)


se esta condicao for valida, diz-se que a EDO 4.12 e uma equacao exacta. Aequacao pode ser resolvida encontrando a funcao F e substituindo em 4.15.

Obviamente que existem varias formas de escrever uma equacao na forma4.12. Se multiplicarmos a equacao toda por uma expressao qualquer, as funcoesM e N serao modificadas. Se uma equacao nao for exacta, por vezes e possıveldescobrir uma funcao (factor integrante) que quando multiplicar a equacaotorna-la-a numa equacao exacta.

Exemplo 4.3Resolva a seguinte equacao

dy

dx=

9x2 + y − 14y − x

(4.17)

(C23) ode2(’diff(y, x) = (9*x^2 + y - 1)/(4*y - x), y, x);

(D23)

2y2 − xy − 3x3 + x = %C

(C24) method;

(D24)

EXACT

4.5.3 Equacoes lineares

Sao equacoes que podem ser escritas na forma

dy

dx+ p(x)y = f(x) (4.18)

o importante e que o unico termo com a variavel dependente y seja uma funcaoso de x a multiplicar a y. Por exemplo, nao poderao aparecer termos da formayy′, y2 ou sin(y).

Pode ser demonstrado que exp(P (x)), onde P (x) e uma primitiva3 dep(x), e um factor integrante. Nomeadamente, e facil ver que

eP dy

dx+ eP (yp− f) = 0 (4.19)

e uma equacao exacta e pode ser resolvida pelo metodo para equacoes exactas.3Isto e, P ′ = p.



dy

dx=

x3 − 2y

x

(C25) ode2(’diff(y,x) = (x^3 - 2*y)/x, y, x);

(D25)

y =

x5

5+ %C

x2

(C26) method;

(D26)

LINEAR

4.5.4 Equacoes homogeneas

Uma equacao de primeira ordem diz-se homogenea se tiver a seguinteforma geral

dy

dx= f

(y

x

)(4.20)

para resolver este tipo de equacao usa-se a substituicao

v =y

x−→ dy

dx= v + x

dv

dx(4.21)

que transforma a equacao numa equacao de variaveis separaveis.


dy

dx=

y2 − xy

x2

(C27) ode2(’diff(y, x)=(y^2 - x*y)/x^2, y, x);

(D27)

%Cx = elog( y−2x

x )−log( yx )

2

(C28) %, radcan;


(D28)

%Cx =√

y − 2x√

y

(C29) solve(%^2, y);

(D29) [y = − 2x

%C2x2 − 1

]

(C30) method;

(D30)

HOMOGENEOUS

O comando radcan foi usado para simplificar a expressao com logaritmose exponenciais.

4.6 Equacao de Bernoulli

Um tipo de equacao diferencial que pode ser reduzida a equacao linear,e a chamada equacao de Bernoulli:

dy

dx+ p(x)yn = f(x)y (4.22)

onde o ındice n e um numero racional, diferente de 0 e de 1. A substituicao

v = y1−n −→ v′ = (1− n)y−ny′ (4.23)

transforma a equacao de Bernoulli numa equacao linear.


dy

dx=

y2 − xy

x2

(C31) ode2(’diff(y, x)=(y^2 - x*y)/x^3, y, x);


(D31)

y =e

1x

%C−∫

e1x

x3dx

(C32) method;

(D32)

BERNOULLI

(C33) odeindex;

(D33)

2

Neste caso, a resposta foi dada em funcao de um integral que nao podeser calculado analıticamente e que define alguma funcao especial, que nao podeser escrita em termos das funcoes elementares. O comando odeindex da o valordo ındice n na equacao 4.22.

4.7 Equacoes de segunda ordem

Nesta seccao vamos usar t para a variavel independente e x para a va-riavel dependente. As derivadas de x(t) serao escritas em forma abreviadacomo x, x, . . .. Esta notacao sera util para fazer analogias com os sistemasmecanicos, onde t sera o tempo e x a posicao.

Uma EDO de segunda ordem e uma equacao com a forma geral

x = f(t, x, x) (4.24)

onde f e uma funcao dada que depende da variavel independente t, da funcaox e da sua primeira derivada x.

4.7.1 Equacoes autonomas

Se a funcao f , na forma geral 4.24, nao depende da variavel independentet, diz-se que a equacao e autonoma:

x = f(x, x) (4.25)


se designarmos de v(x) a funcao x(t), a equacao assume uma forma mais simples

vdv

dx= f(x, v) (4.26)

que e uma EDO de primeira ordem. Cada solucao dessa EDO sera uma (oumais) funcao g(x) que representa a v em funcao da variavel x. Para calcularx(t) resolve-se a seguir

dx

dt= g(x) (4.27)

que e tambem uma equacao autonoma, mas de primeira ordem. O problema eque em alguns casos nao e facil escrever v na forma explıcita g(x).

Na analogia mecanica, v e a velocidade e a funcao f na equacao 4.25 ea forca resultante, por unidade de massa.

Exemplo 4.7A aceleracao de uma partıcula, em funcao do tempo e x = −3x − 5x. Noinstante t = 0, a partıcula parte do repouso no ponto x = 1. Calcule a posicaoe a velocidade da partıcula em funcao do tempo, para t > 0.

Em funcao da velocidade v = x, a equacao da aceleracao e

vv = −3x− 5v

A condicao inicial e v = 0 em x = 1. A solucao obtida com Maxima e:

(C34) eq3: v*’diff(v,x) = -3*x - 5*v;

(D34)

vdv

dx= −3x− 5v

(C35) sol3: ode2(eq3, v, x)$

(C36) ic1(sol3, x=1, v=0), radcan;

(D36)

25

2√

13 x(5√

13 + 19) 5

2√

13 3√

13−52√

13

=x((

5−√

13)x + 2v

) 52√

13((√13 + 5

)x + 2v

) 52√

13√

3x2 + 5vx + v2


Nessa expressao nao e facil obter v(x) em forma explıcita. Para cal-cular x(t) neste caso sera mais facil resolver directamente a equacao inici-al

(C37) eq4: ’diff(x, t, 2) = -3*x -5*’diff(x, t);

(D37)d2x

dt2= −5

dx

dt− 3x

(C38) sol4: ode2(eq4, x, t)$

(C38) ic2(sol4, t=0, x=1, diff(x,t)=0);

(D39)

x =

(5√

13 + 13)e

(√

13−5)t

2

26−(5√

13− 13)e

(−√

13−5)t

2

26

Repare no uso da funcao ic2 para especificar as condicoes iniciais deuma equacao de segunda ordem. Sao precisos quatro, argumentos: a solucaogeral, o valor da variavel onde sao conhecidas as condicoes iniciais, e os valoresda funcao e da sua derivada nesse ponto. Para dar o valor da derivada nao fazfalta usar apostrofo.

No exemplo anterior, a substituicao de x por uma funcao v(x) complicouo problema. A equacao 4.26 sera muito util nos casos em que a aceleracao x

nao depende da velocidade x, isto e, forcas que nao dependem da velocidade.Veremos um exemplo.

Exemplo 4.8A aceleracao de uma partıcula, em funcao do tempo e x = −3x. No instantet = 0, a partıcula parte do repouso no ponto x = 1. Calcule a posicao e avelocidade da partıcula em funcao do tempo, para t > 0.

(C40) eq5: v*’diff(v,x) = -3*x;

(D40)

vdv

dx= −3x

(C41) sol5: ode2(eq5, v, x)$


(C42) ic1(sol5, x=1, v=0);

(D42)

−v2

6=

x2 − 12

(C43) vel: solve(%,v);

(D43) [v = −

√3√

1− x2, v =√

3√

1− x2]

obtivemos duas solucoes para a velocidade. A equacao D47 e a lei daconservacao da energia mecanica. A equacao diferencial para a posicao emfuncao da velocidade e:

(C44) eq6: ’diff(x,t) = v;

(D44)dx

dt= v

podemos substituir uma das expressoes obtidas para a velocidade, eresolver a equacao com as condicao incial dada

(C45) sol6: ode2(ev(eq6, vel[1]), x, t)$

(C46) solve(ic1(sol6, t=0, x=1), x);

(D46) [x = cos

(√3t)]

se usarmos a segunda expressao obtida para a velocidade, obtem-se amesma resposta para x em funcao de t.

4.7.2 Equacoes lineares

A forma geral de uma EDO linear de segunda ordem e

x + p(t)x + q(t)x = f(t) (4.28)

onde p(t), q(t) e f(t) sao tres funcoes que dependem de t. Se a funcao f fornula, diz-se que a equacao e linear homogenea.


A solucao geral de uma EDO linear de segunda ordem depende de duasconstantes independentes. As solucoes particulares definem um espaco vectorialde funcoes, com dimensao igual a 2. Para fixar os valores das duas constantesarbitrarias podemos impor condicoes iniciais, isto e, fixar o valor da funcaoe da sua derivada num mesmo ponto t0

x(t0) = t0 x(t0) = v0 (4.29)

se as tres funcoes p, q e f estiverem bem definidas, existira sempre uma solucaounica com essas condicoes iniciais. Para substituir as condicoes iniciais emMaxima usa-se a funcao ic2, como ja vimos no exemplo 4.7.

Uma outra forma de obter uma solucao particular e impor condicoesfronteira, que consiste em fixar o valor da funcao x(t) em dois pontos diferen-tes.

x(t1) = x1 x(t2) = x2 (4.30)

nesse caso, usa-se em Maxima a funcao bc2, como veremos no exemplo que sesegue.

Exemplo 4.9A corrente I, no circuito LCR representado no diagrama, verifica a equacaodiferencial (t em milisegundos e I em miliamperes)

I + 4I + 4I = cos(2t)

Os valores medidos da corrente, em t = 0 e t = π, sao 0 e 5 respectivamente.Encontre uma expressao para a corrente em funcao do tempo

4 kΩ

1 H250 nF

0.5 sin(2t) V

(C47) eq7: ’diff(I,t,2) + 4*’diff(I,t) + 4*I = cos(2*t);

(D47)d2I

dt2+ 4

dI

dt+ 4I = cos (2t)


(C48) sol7: ode2(eq7, I, t);

(D48)

I =sin (2t)

8+ (%K2t + %K1) e−2t

(C49) bc2(sol7, t=0, I=0, t=%pi, I=5);

(D49)

I =sin (2t)

8+

5te2π−2t

π

4.8 Espaco de fase

As equacoes lineares de ordem superior a dois podem ser resolvidas pormetodos analıticos semelhantes aos metodos usados para o caso de ordem 2. Asequacoes nao-lineares de ordem 2 ou superior geralmente tem que ser resolvidasem forma numerica. A resolucao numerica de equacoes diferenciais e o temado proximo capıtulo; mas antes de prosseguir, vamos definir alguns conceitosimportantes.

Todas as equacoes de ordem 1 ou superior que temos estudado nestecapıtulo podem ser escritas na forma de um sistema de EDOs autonomas deprimeira ordem:

x1 = f1(x1, . . . , xn) (4.31)

x2 = f2(x1, . . . , xn)... (4.32)

xn = fn(x1, . . . , xn)

As variaveis do problema, (x1, x2, . . . , xn), designam-se de variaveis deestado, e definem um espaco de n dimensoes, chamado espaco de fase.

Vejamos alguns exemplos. Para escrever a equacao y′ = x2 + 3y comoum sistema autonomo, basta substituir x = t na EDO; a definicao x = t podeser escrita como equacao diferencial, derivando os dois lados em ordem a t, e


obtemos o sistema autonomo:

y = x2 + 3y (4.33)

x = 1

As variaveis de estado sao x e y; o espaco de fase ja foi representado na figura4.1.

A equacao de segunda ordem x = −3x − 5x, do exemplo 4.7, escreve-se facilmente como um sistema autonomo se substituirmos a aceleracao peladerivada da velocidade:

v = −3x− 5v (4.34)

x = v

O espaco de fase e um espaco a duas dimensoes e as variaveis de estado sao aposicao x e a velocidade v.

Finalmente, para escrever a equacao do circuito do exemplo 4.9 como umsistema autonomo de primeira ordem, para alem de definir I como uma outravariavel de estado J , sera tambem preciso considerar o tempo t como variavelde estado, por exemplo, s = t. Assim, a equacao e equivalente ao sistema:

J = cos(2s)− 4J − 4I (4.35)

I = J

s = 1

o espaco de fase, neste caso, e um espaco com 3 dimensoes.

4.9 Problemas

1. Demonstre em cada caso que a funcao y(x), na coluna direita, e solucaoda equacao diferencial, na coluna esquerda (a e c sao constantes).

(a)dy

dx=

x√x2 + a2

(a 6= 0) y =√

x2 + a2

(b)14

(d2y

dx2

)2

− xdy

dx+ y = 1− x2 y = x2

(c) yy′ = e2x y2 = e2x


(d) y′ =y2

xy − x2y = c ey/x

Sugestao: quando for dada a funcao y = f(x) explicitamente, podesubstituir directamente na equacao e mostrar que se chega a uma identi-dade. Se nao tiver y = f(x), derive a equacao que define y e simplifiqueate obter a equacao diferencial dada.

2. Encontre a solucao geral das seguintes equacoes diferenciais ordinarias eem cada caso diga a que classe pertence a equacao.

(a)dy

dt=

y2 − 2ty

y2

(b) (2y + ex cos y)y′ = −ex sin y

(c)dy

dx=

x2 − 1y2 + 1

(d)dy

dt+ 2ty = 2t3

√y

(e)dy

dx=

x3 − 2y

x

(f)dy

dx=

x

x2y + y3

(g)dy

dx=

x(2y + 1)y − x2

(h)dy

dx=

y − x2

y2 − x

3. Faca o desenho do campo de direccoes de cada uma das equacoes que seseguem, e desenhe a solucao que verifica a condicao inicial dada (encontreum domınio que mostre bem o comportamento do campo)

(a)dy

dt+ y = 1 + t2 y(1) = 2

(b)dy

dx= − x + y

x + 2yy(2) = 3

(c)dy

dtcos y = − t sin y

1 + t2y(1) =

π

2

4. Encontre a solucao das seguintes equacoes de segunda ordem, com ascondicoes iniciais, ou de fronteira, dadas

(a) y′′ + y′ − 2y = 3− 6x y(1) = 2 y′(1) = 1

(b) y′′ − y = 0 y(0) = 2 y(1) = 3

(c) x2y′′ + xy′ − 4y = x2 + x4 y(−1) = 0 y′(−1) = −2


5. Resolucao numerica de equacoes diferenciais

Os metodos de resolucao numerica consistem em substituir o problemade encontrar a solucao x(t), num domınio real, por um outro problema dis-creto onde e calculada uma sequencia xi(ti), i = 0, 1, . . . que seja uma boaaproximacao a funcao de variavel real. A equacao diferencial

x = f(t, x) (5.1)

conduz ao sistema discreto:

ti = a, a + h, a + 2h, a + 3h, . . . , b

xi = x(a), x(a + h), . . . , x(ti)

xi = g(h, x1, x2, . . . , xi−1)

a ultima dessas equacoes obtem-se a partir de x = f , com alguma aproximacaopara estimar a derivada em funcao dos valores discretos de t e de x. Vamos veralguns metodos de aproximacao.

5.1 Metodo de Euler

A definicao da derivada x e

x(t) = f(t, x(t)) = limh→0

x(t + h)− x(t)h

(5.2)

se h for suficientemente pequeno, uma boa aproximacao sera

x(t + h) = x(t) + hf(t, x(t)) (5.3)

xi+1 = xi + hf(ti, xi) (5.4)


Para ilustrar o metodo, consideremos um problema simples com solucaoja conhecida:

x = x x(0) = 1 (5.5)

vemos facilmente que a solucao e x = et, mas vamos ver como seria resolvidopelo metodo de Euler. A funcao f e x neste caso e a equacao 5.4 fica xi+1 =xi + hxi.

Em Maxima, um possıvel programa para calcular a solucao pelo metodode Euler e

(C1) prog1(h) := (x:1, for t:h step h thru 1 do x:x+h*x)$

O programa aceita como argumento o valor de h. O comando for in-crementa o valor de t desde h ate o valor final 1, com intervalos iguais a h, eem cada passo calcula o valor de x nesse intervalo, usando a aproximacao 5.4.

Repare na sintaxe para incluir mais do que um comando na definicao dafuncao; os comandos escrevem-se entre parentesis, separados por vırgulas, massem vırgula apos o ultimo comando. Nao devem ser usados outros caracterescomo “;” ou “$” para terminar um comando entre os parentesis.

Vejamos os resultados obtidos, e o erro em cada caso, para varios valoresde h

(C2) for h in [0.1, 0.01, 0.001, 0.0001, 0.00001] do

(prog1(h), er: float(%e-x),

print("h=", h, " x= ", x, " erro= ", er))$

h= 0.1 x= 2.5937424601 erro= 0.12453936835905

h= 0.01 x= 2.678033494476758 erro= 0.04024833398229

h= 0.001 x= 2.714209722513383 erro= 0.00407210594566

h= 0.0001 x= 2.718145926825227 erro= 0.00013590163382

h= 0.00001 x= 2.718268237174493 erro= 0.00001359128455

Convem generalizar o programa admitindo tambem como parametrosde entrada a funcao f que define a equacao diferencial e os valores iniciais.Tambem e aconselhavel usar variaveis locais dentro do programa (usando ocomando block), e guardar os resultados intermedios de t e x para poderestudar o comportamento da solucao. Finalmente, em vez de dar o valor doincremento h, daremos o numero de iteracoes que o programa devera realizar

(C3) eulermet(f, t0, x0, tf, n) :=

block

Resolucao numerica de equacoes diferenciais 57

([t:t0, x:x0, l:[t0, x0], h:(tf - t0)/n,

display2d: false, numer: true],

for i thru n do

(x: x + h*ev(f), t: t0 + i*h, l: append(l,[t,x])),

l

)$

O programa usa o comando block que permite definir variaveis locais(uma lista imediatamente no inıcio do bloco) com possıveis valores iniciais. Nes-te caso as variaveis locais sao t e x, com valores iniciais iguais aos parametros deentrada t0 e x0, uma lista l, inicializada com a condicao inicial, que guardaraos valores de t e x calculados em cada passo, e duas variaveis que controlama forma como os resultados serao apresentados; display2d devera ter um va-lor falso, para que nao sejam inseridos espacos em branco entre numeros esımbolos que produz uma ma interpretacao de valores negativos na saıda, e ovalor verdadeiro de numer faz com que sejam apresentados resultados de vırgulaflutuante, quando possıvel. O ındice i tambem sera definido implicitamentecomo variavel local, por fazer parte do comando for.

O incremento h e calculado dividindo o intervalo de t0 ate tf em n

partes iguais. Repare que primeiro devera ser actualizado o valor de x, antesde se actualizar t, pois a funcao f(t, x) deve ser calculada em (ti, xi) e nao em(ti+1, xi).

Para acrescentar os valores de t e x a lista l, usou-se o comando append.No fim do programa, escreve-se a lista completa de valores (t, x). Podemosusar o comando openplot curves para mostrar um grafico com os resultados(ver figura 5.1). Na funcao f(t, x) x e t deverao ser variaveis e nao constantes;como em C2 ja demos valores a x sera preciso comecar por fazer uma limpeza.

(C4) kill(x)$

(C5) for n IN [10, 20, 1000] do

concat(l, n) :: eulermet(x, 0, 1, 5, n)$

(C6) openplot_curves([["plotpoints 1 nolines"], l10,

l20, l1000]);

No ultimo caso, n = 1000, os resultados ja sao indistinguıveis da solucaoexacta no grafico. E preciso ter algum cuidado com a sintaxe do comandoopenplot curves. Quem usar a interfaz Xmaxima, nao podera terminar o co-mando com $ pois a interface fica pendurada (a comunicacao do Xmaxima comcomandos externos e um pouco problematica na versao 5.9.0). O argumento


1 2 3 4

50

100

t

x

Figura 5.1: Aproximacao numerica a solucao de x = x, com diferentes numerosde passos.

do comando devera ser sempre uma lista, com uma ou mais listas dentro. Umadas listas pode ter opcoes de configuracao para o openplot; as opcoes devemvir dentro de aspas e chavetas, deixando espaco entre chavetamkdir s vizinhas.

Usemos agora o programa anterior para resolver uma equacao mais com-plexa, a equacao nao linear x = t − x2. Vamos resolver a equacao para cincocondicoes iniciais diferentes:

(C7) i:0$

(C8) for x0 in [3, 1, 0, -0.7] do

(i:i+1,

concat(l, i) :: eulermet(t-x^2, 0, x0, 9, 180))$

(C9) l5: eulermet(t-x^2, 0, -0.75, 1.8, 36)$

(C10) openplot_curves([["plotpoints 1 nolines"], l1, l2,

l3, l4, l5]);

O resultado e apresentado na figura 5.2. Foi usado o comando concat

que permite ligar duas constantes para criar uma constante de texto, e o ope-


rador :: que permite atribuir um valor ao nome de variavel obtido com ocomando no lado esquerdo. O caso x0 = −0.75 teve que ser tratado a parte,porque nesse caso a solucao decresce rapidamente para −∞ dando origem aerros numericos se o intervalo de resolucao for muito grande.

1 2 3 4 5 6 7 8

-2

-1

1

2

t

x

Figura 5.2: Solucoes particulares da equacao x = t− x2, com x = 3, 1, 0, −0.7e −0.75, em t = 0.

Neste exemplo podemos ver um comportamento tıpico das equacoes naolineares, nomeadamente a sua sensibilidade aos valores iniciais. Com dois pon-tos iniciais muito poximos, (0,−0.7) e (0,−0.75), obtemos dois comportamentostotalmente opostos.

Seria de esperar que o tempo que o programa eulermet demorasse paracalcular 2n interaccoes fosse o dobro do tempo necessario para n interaccoes.Mas nao e assim, o que acontece e que o tempo cresce rapidamente e por voltadas 1000 interacoes pode ser ja muito elevado, mesmo num computador rapido.Isto e devido ao uso do comando append que faz com que seja preciso carregarem memoria uma lista que cresce cada vez mais.

Para acelerar o programa, no caso de milhares de iteracoes, em vez deguardar os dados intermedios numa lista, podemos arquiva-los num ficheiro ex-


terno, e no fim criar a lista com todos os valores. Um programa que implementaesse metodo e assim:

(C11) ode_euler(f, t0, x0, tf, n) :=

block

([t:t0, x:x0, h:(tf - t0)/n, euler_dat,


with_stdout

(

"euler.dat", print("euler_dat: [ ", t, ", ", x),

for i thru n do

(x:x + h*ev(f), t:t0 + i*h, print(", ", t, ", ", x)),

print("]$")

),

batchload("euler.dat"),

euler_dat

)$

Podemos agora resolver novamente a equacao x = t − x2 mas para umnumero elevado de passos.

(C12) caos: ode_euler(t-x^2, 0, 0, 900, 18000)$

(C13) openplot_curves([["plotpoints 1 nolines

pointsize 0.7"], caos]);

O resultado obtido mostra um comportamento surpreendente (figura 5.3,lado esquerdo); em t ≈ 400 a solucao comeca a oscilar entre dois valores; emt ≈ 700, a solucao oscila entre varios valores em forma caotica. E de salientarque esse comportamento caotico nao e devido a erros numericos aleatorios, poiscada ponto e calculado em forma precisa, a partir do ponto anterior, com asequacoes de recorrencia xi+1 = xi + h(ti − x2

i ), ti+1 = ti + h. Por exemplo,com t0 = 750, x0 = 20 e h = 0.05, obtem-se para xi a sequencia de valores: 20,37.5, 4.69, 41.10, . . .

No entanto, o que acontece neste problema e que a equacao discretaxi+1 = xi + h(ti − x2

i ) deixa de ser uma boa aproximacao para a equacaocontınua x = t − x2. A solucao exacta da equacao contınua nao apresentaoscilacoes, como podemos conferir usando outros metodos melhores do que ometodo de Euler, ou aumentando o numero de passos. O lado direito da figura


0 900-60

60

t

x

0 900-60

60

-30

0

30

300 600 0 900-60

60

tx

0 900-60

60

-30

0

30

300 600

Figura 5.3: Solucao caotica para n = 18000 (esquerda), e desaparicao da regiaocaotica usando uma precisao maior n = 26000.

5.3 mostra o resultado com n = 26000. Ate t ≈ 800 obtem-se uma curvacontınua.

5.2 Metodo de Runge-Kutta

A aproximacao 5.4 usada na seccao anterior calcula o valor da funcaono ponto ti + h a partir da derivada no ponto ti. Para melhorar o metodo serapreciso usar uma melhor aproximacao ao valor medio da derivada, no intervalo(ti, ti + h). No metodo de Runge-Kutta de quarta ordem, o valor da derivada,d, obtem-se a partir da media das derivadas em 4 pontos diferentes, com pesosdiferentes:

d1 = f(ti, xi) (5.6)

d2 = f(ti + h/2, xi + (h/2)d1) (5.7)

d3 = f(ti + h/2, xi + (h/2)d2) (5.8)

d4 = f(ti + h, xi + hd3) (5.9)

d =16(d1 + 2d2 + 2d3 + d4) (5.10)


com o valor da derivada, calcula-se o valor de xi+1

xi+1 = xi + hd (5.11)

O programa que se segue implementa esse metodo. Vamos aproveitarpara mostrar outra forma de escrever programas em Maxima, usando o coman-do read para pedir os argumentos interactivamente, em vez de serem aceitescomo argumentos da funcao. Para simplificar as contas, em vez de se calcu-larem as derivadas, calculam-se as derivadas multiplicadas por h: k1, k2, k3 e k4

(C14) rungekutta():=

block

(

[f, t, x, t0, tf, n, h, k1, k2, k3, k4, rk4_dat,


f: read("Escreva a func~ao f(t,x); termine com ;"),

print("f(t,x) = ", f),

t0: read("Valor inicial de t?"),

x: read("Valor inicial de x?"),

tf: read("Valor final de t?"),

n: read("Numero de passos?"),

t: t0, h: (tf - t0)/n,

with_stdout

(

"rk4.dat", print("rk4_dat: [ ", t, ", ", x),

for i thru n do

(

k1: h*ev(f),

k2: h*ev(f, t = t + h/2, x = x + k1/2),

k3: h*ev(f, t = t + h/2, x = x + k2/2),

k4: h*ev(f, t = t + h, x = x + k3),

x: x + (k1 + 2*k2 + 2*k3 + k4)/6,

t: t0 + i*h,

print(", ", t, ", ", x)

),

print("]$")

),

batchload("rk4.dat"),


rk4_dat

)$

Podemos usar este programa para encontrar a solucao exacta do proble-ma representado na figura 5.3: x = t − x2, t0 = 0, x0 = 0, tf = 900. Comn = 20000 ja e possıvel obter uma boa solucao (figura 5.4). Se o numero depassos nao for suficientemente elevado, Maxima queixar-se-a de valores naonumericos, quando o valor de x crescer demasiado.

0 9000

30

t

x

0 9000

30

5

10

15

20

25

300 600

Figura 5.4: Solucao de x = t− x2 obtida pelo metodo de Runge-Kutta.

5.3 Sistemas de equacoes de primeira ordem

No capıtulo anterior vimos como escrever qualquer equacao diferencialordinaria, ou sistemas delas, como um sistema de equacoes diferencias de pri-meira ordem. Para resolver sistemas de equacoes, basta generalizar o metodode Runge-Kutta, admitindo que nas equacoes 5.10 a funcao f e as derivadas d1,


d2, d3 e d4 sao vectores com m componentes, onde m e o numero de equacoesno sistema.

(C15) ode_rk4(odes, st, init, domain) :=

block

(

[t, t0, xv:[], f:[], var:[], m, n, h,

k1:[], k2:[], k3:[], k4, k:[], rungekutta_dat,


if listp(odes) then f: copylist(odes) else f: [odes],

if listp(st) then var: copylist(st) else var: [st],

if listp(init) then xv: copylist(init) else xv: [init],

k1: copylist(xv), k2: copylist(xv), k3: copylist(xv),

k4: copylist(xv), k: copylist(xv),

m: length(f),

n: domain[4],

h: (domain[3] - domain[2])/n,

t0: domain[2],

t: t0,

with_stdout

(

"rungekutta.dat", print("rungekutta_dat: [ ", t),

for j thru m do print(", ", xv[j]),

for i thru n do

(

for j thru m do

(

k1[j]: h*ev(f[j], domain[1] = t),

for l thru m do k1[j]: ev(k1[j], var[l] = xv[l])

),

for j thru m do

(

k2[j]: h*ev(f[j], domain[1] = t + h/2),

for l thru m do k2[j]:

ev(k2[j], var[l] = xv[l] + k1[l]/2)

),

for j thru m do

(


k3[j]: h*ev(f[j], domain[1] = t + h/2),


ev(k3[j], var[l] = xv[l] + k2[l]/2)

),

for j thru m do

(

k4[j]: h*ev(f[j], domain[1] = t + h),


ev(k4[j], var[l] = xv[l] + k3[l]),

k[j]: (k1[j] + 2*k2[j] + 2*k3[j] + k4[j])/6

),

for j thru m do xv[j]: xv[j] + k[j],

t: t0 + i*h,

print(", ", t),

for j thru m do print(", ", xv[j])

),

print("]$")

),

batchload("rungekutta.dat"),

rungekutta_dat

)$

Este programa pode ser usado para resolver uma unica EDO de primei-ra ordem, ou um sistema dessas equacoes. Na proxima seccao veremos umexemplo.

5.4 Eliminacao de singularidades

Os metodos para resolver equacoes diferenciais estudados nas seccoesanteriores calculam o valor da solucao a partir do valor da derivada num pontoinicial. Se a derivada no ponto inicial for infinita, o metodo falha. Por exemploconsideremos o problema de valor inicial

dy

dx= −x

yy(3) = 0 (5.12)


no ponto inicial, como y = 0, a derivada e infinita e nao vai ser possıvel calcularo valor de y para um valor de x perto de 3.

Em geral, no diagrama do campo de direccoes de uma EDO, nos pontosonde o declıve for vertical, os metodos numericos nao dao bom resultados.Mas o problema pode ser resolvido facilmente, se definirmos a equacao comoum sistema livre de singularidades. Por exemplo, a equacao 5.12 e equivalenteao sistema de equacoes

dx

dt= y (5.13)

dy

dt= −x (5.14)

-10 2

2

10

-10

-10

-5

-5

0

0

5

5

-5 -5

0 0

5 5

10

Figura 5.5: Campo de direccoes do sistema x = y, y = −x e curva integral em(3, 0).

A figura 5.5 mostra o campo de direccoes desse sistema. O campo seriaidentico se tivesse sido definido pela equacao 5.12; no entanto, com essa defi-nicao, se tentassemos obter a curva integral que passa por (3, 0) nao teriamos


sucesso (experimente). O mesmo campo de direccoes definido com o sistemaanterior, permite obter facilmente as curvas integrais

Com o nosso programa ode rk4 podemos obter a solucao do sistema deequacoes 5.13 e 5.14

(C16) a: ode_rk4([y,-x],[x,y],[3,0],[t,0,7,700])$

A lista a contem uma sequencia de pontos (t, x, y). Para desenhar umgrafico em duas dimensoes sera preciso extrair os valores que queremos; porexemplo, para desenhar um grafico de y vs x e um outro grafico de x vs t, usare-mos:

(C17) b:makelist(([a[3*i+2],a[3*i+3]]),i,0,700)$

(C18) c:makelist(([a[3*i+1],a[3*i+2]]),i,0,700)$

e a seguir podemos obter os graficos das listas b e c usando o comandoopenplot curves.

-3 -2 -1 1 2 3

-3

-2

-1

1

2

3

x

y

1 2 3 4 5 6

-3

-2

-1

1

2

3

t

x

Figura 5.6: Solucao do sistema x = y, y = −x, com x(0) = 3 e y(0) = 0.

5.5 Problemas

1. Escreva o programa eulermet da seccao 5.1, grave-o num ficheiro, euse-o para reproduzir a figura 5.2.


2. Escreva o programa ode euler da seccao 5.1, grave-o num ficheiro, euse-o para obter o grafico da solucao (no intervalo pedido) dos seguintesproblemas com condicao inicial:

(a) y =t2 − y2

yy(0) = 1 0 ≤ t ≤ 3

(b) x = cos t + 2x x(−2) = 0 − 2 ≤ t ≤ 1

(c) y + y = 1 + t2 y(0) = 2 0 ≤ t ≤ 4

3. Carregue o programa ode rk4, usando o ficheiro ode rk4.mac disponıvelcom estes apontamentos, e use-o para resolver o sistema

x′1 = x1 + x2

x′2 = 4x1 + x2

com condicoes iniciais x1(0) = 1, x2(0) = 2, no intervalo 0 ≤ t ≤ 6.Desenhe os graficos de x1 e x2 em funcao do tempo.

4. Usando o programa ode rk4 (ver problema anterior), resolva os seguintesproblemas, de segunda ordem:

(a) (x2 +1)y′′− (x+3)y′+ y = 0 y(−3) = 0, y′(−3) = 1, −3 ≤x ≤ 1

(b) y′′+xy′+x2y = sinx, y(−1) = 0, y′(−1) = −2 −1 ≤ x ≤ 5

5. Em cada caso, escreva a EDO como um sistema livre de singularidades,faca o desenho do campo de direccoes e desenhe a solucao que verifica acondicao inicial dada (encontre um domınio que mostre bem o compor-tamento do campo)

(a)dy

dx= − x + y

x + 2yy(2) = 3

(b)dy

dtcos y = − t sin y

1 + t2y(1) =

π

2


6. Sistemas oscilantes

6.1 O oscilador harmonico simples

A figura 6.1 mostra uma mola numa posicao vertical com um extremofixo. Quando uma esfera de massa m e pendurada no outro extremo, a molaestica uma distancia ∆x. O lado direito da figura mostra a forca Fk produzidapela mola, em funcao da posicao x da esfera, segundo o eixo de coordenadasdefinido no diagrama da esquerda. A forca elastica da mola e directamente pro-porcional ao alongamento da mola, com declıve constante −k (o sinal negativoe devido a que a forca e oposta ao alongamento).

0

x

x∆ x

F

−mg

−∆x

Fk

Fr

Figura 6.1: Esfera suspensa por uma mola vertical.

A forca resultante, Fr, obtem-se somando o peso, mg, mais a forcaelastica Fk. Como se mostra no grafico 6.1, Fr e outra recta deslocada umadistancia mg no eixo das ordenadas. Assim, em funcao da coordenada x, aforca resultante e Fr = −kx. A aceleracao da esfera e x. A segunda lei de


Newton implica

mx = −kx (6.1)

Para que fique explıcito o facto de que as duas constantes k e m saopositivas, podemos escrever a equacao na forma

x = −ω2x (6.2)

onde

ω =

√k

m(6.3)

A solucao para x em funcao de t esperamos que seja uma funcao osci-lante. Podemos mostrar que de facto a solucao tem a forma x = Acos(Bt+C):Derivando duas vezes em ordem a t obtemos

(C1) sol1: A*cos(B*t + C);

(D1)

A cos (Bt + C)

(C2) diff(sol1, t, 2);

(D2)

−AB2 cos (Bt + C)

Que e igual a funcao inicial x, multiplicada por −B2. Independentemen-te dos valores das constantes A e C, se B for igual a ω, a funcao sera solucaoda equacao 6.2. Logo a solucao geral da equacao 6.2 e:

x(t) = A cos(wt + C) (6.4)

onde A e uma constante, designada de amplitude, que representa o valorabsoluto maximo do deslocamento x. A constante C, e designada de constantede fase. Tanto A como C dependem das condicoes iniciais do problema (porexemplo, deslocamento e velocidade inicial).

O espaco de fase esta formado pelas variaveis x e v = x. O sistemaautonomo de primeira ordem, equivalente a equacao 6.2 e

x = v (6.5)

v = −ω2x (6.6)

Sistemas oscilantes 71

No caso particular ω = 1, este sistema ja foi estudado em pormenor naseccao 5.4, ja que as equacoes 6.5 e 6.6 tem a mesma forma das equacoes 5.13e 5.14. As figuras 5.5 e 5.6 mostram o espaco de fase e a evolucao de x emfuncao do tempo. Para uma constante ω qualquer, e facil ver que a funcao x(t)(equacao 6.4) completa uma oscilacao completa quando ωt aumentar em 2π, eo valor maximo de v (derivada de x) e Aω. A figura 6.2 mostra a evolucao dex e v e a trajectoria do sistema no espaco de fase.

x

v

A

Aω

t

x, v

P

x(t) v(t)

Figura 6.2: Evolucao do oscilador harmonico simples, no espaco de fase e nodomınio do tempo.

O perıodo P e o tempo que oscilador demora a completar uma oscilacaoe verifica a equacao

P =2π

ω(6.7)

A frequencia angular ω tem unidades de inverso do tempo, e e definida pelaequacao 6.3.

A trajectoria no espaco de fase e uma elipse1 com semi-eixos de compri-mento A e Aω. A equacao da elipse e

x2

A2+

v2

A2ω2= 1 (6.8)

Subsituindo a definicao de ω (equacao 6.3) vemos que esta equacao expressa aconservacao da energia mecanica:

12kx2 +

12mv2 =

12kA2 (6.9)

1Isto pode ser demonstrado matematicamente, resolvendo a equacao 6.2 pelo metodo descritona seccao 4.7.1.


6.2 Osciladores amortecidos

A figura 6.3 mostra o sistema de suspensao de um automovel. Dentro damola ha um amortecedor, isto e, um cilindro cheio de oleo, que produz umaforca oposta ao movimento, directamente proporcional a velocidade: Fv = −av,onde a e uma constante.

x

0

t

x

b<2ω

b=2ω

b>2ω

Figura 6.3: Oscilador harmonico amortecido.

A equacao de movimento, obtem-se somando Fv no lado direito daequacao 6.1. Apos dividir os dois lados da equacao por m, obtemos a equacaodiferencial do oscilador harmonico amortecido:

x + bx + ω2x = 0 (6.10)

onde b = a/m. Trata-se de uma equacao linear, de segunda ordem, que podeser resolvida analıticamente usando o comando ode2.

(C3) eq2: ’diff(x,t,2) + b*’diff(x,t) + w^2*x = 0$

(C4) ode2(eq2, x, t);

Is (2*w - b)*(2*w + b) positive, negative, or zero?

negative;

Sao obtidas 3 tipos de solucoes, segundo (2ω − b) seja nula, positiva ounegativa. O lado direito da figura 6.3 mostra os 3 tipos de solucoes. Quandob < 2ω, dizemos que o amortecimento e fraco, e a solucao e uma funcao que


oscila mas com amplitude que decresce rapidamente. O caso b = 2ω e desig-nado de amortecimento crıtico e conduz a uma solucao que se aproximapara x = 0 (com algumas condicoes iniciais, x pode chegar a mudar de sinalantes de se aproximar para zero). Finalmente, b > 2ω corresponde ao caso desobreamortecimento, em que x decresce rapidamente para zero.

Um oscilador real, por exemplo o sistema da figura 6.1, tem sempre umtermo de amortecimento devido ao atrito com o ar. O que acontece e queo amortecimento pode ser muito fraco, o que faz com que a amplitude dasoscilacoes diminuia lentamente. Assim, a equacao 6.10 e mais realista do quea equacao 6.2. O respectivo sistema autonomo de primeira ordem e

x = v (6.11)

v = −ω2x− bv (6.12)

A figura 6.4 mostra a evolucao do sistema, no espaco de fase (x, v) e nodomınio do tempo.

-10 2

2

10

dx/dt = y , dy/dt = -x-0.3*y

-10

-10

-5

-5

0

0

5

5

-5 -5

0 0

5 5

10 10

For villate Tue Nov 11 13:04:17 WET 2003

5

2

dy/dt=-x-0.3*y, dx/dt=y

10

10

20

20

-5 -5

0 0

5 5


Figura 6.4: Evolucao do oscilador harmonico com amortecimento fraco.

No exemplo do automovel, para que este nao oscile cada vez que passesobre um obstaculo, os amortecedores sao desenhados de forma a produziremsobreamortecimento. Devido ao desgaste, o oleo do amortecedor comeca aperder pressao e quando o sistema entrar no regime de amortecimento fraco, ocarro oscila quando o mecanico empurra a carroceria para baixo, e estara nahora de mudar os amortecedores.


6.3 Osciladores forcados

Uma forca externa F introduz um termo adicional no lado direito daequacao linear 6.10, tornando-la nao homogenea. Por exemplo, consideremosuma forca externa que tambem oscila em forma harmonica, com amplitude C

e frequencia angular d

x + bx + ω2x = c cos(dt) (6.13)

onde c e igual a C/m. Esta equacao ja nao e autonoma. Isso implica queo tempo t tambem sera uma variavel de estado e o espaco de fase tem tresdimensoes (t, x, v). O comando plotdf nao permite representar espacos defase em 3 dimensoes; vamos usar o comando ode2 para resolver a equacaodiferencial num exemplo em particular.

(C5) eq3: ’diff(x,t,2) + ’diff(x,t)/4 + x = cos(t/2);

(D5)

d2x

dt2+

14

dx

dt+ x = cos

(t

2

)

(C6) sol3: ic2(ode2(eq3,x,t), t=0, x=5, diff(x,t)=0)$

(C7) plot2d(rhs(sol3),[t,0,30]);

O resultado aparece na figura 6.5. O movimento obtido e a sobreposicaoda oscilacao forcada produzida pela forca externa e a a oscilacao propria dosistema. Passado um tempo suficientemente grande, a oscilacao propria dosistema acaba por desaparecer, devido ao amortecimento, e o sistema oscilacom a mesma frequencia da forca externa (regime estacionario). No entanto,a amplitude de oscilacao nesse regime estacionario depende da relacao entre afrequencia externa d e a frequencia propria ω. O valor maximo da amplitude eobtido quando d = ω. Nesse caso diz-se que existe ressonancia entre a forcaexterna e o sistema.


5 25

1

5

10

10

20

20

-2.5 -2.5

0 0

2.5 2.5

5 5

Figura 6.5: Oscilacoes forcadas.

6.4 O pendulo simples

Outro exemplo de sistema oscilante e o pendulo simples (figura 6.6).Consiste num objecto pequeno, pendurado por um fio ou por uma barra rıgidamuito leve. As forcas que actuam sobre o objecto sao a tensao no fio, T , e opeso do objecto, mg.

Convem usarmos um sistema de coordenadas polares para analisarmoso movimento do objecto. O lado direito da figura 6.6 mostra as coordenadaspolares (r, θ) do objecto, e os versores ortogonais ur e uθ que apontam nadireccao em que as duas coordenadas aumentam. Os versores polares naopermanecem constantes, mas dependem do angulo θ. Em funcao dos versoresi e j, na direccao dos eixos x e y, temos

ur = cos(θ) i + sin(θ) j (6.14)

uθ = − sin(θ) i + cos(θ) j (6.15)

derivando essas expressoes em funcao de t, e tendo em conta que θ depende det e a sua derivada e a velocidade angular θ, obtemos

ur = θuθ (6.16)

uθ = −θur (6.17)


mg

Tθ

x

y

rθ

u

u

θ

r

Figura 6.6: Pendulo simples e coordenadas polares.

O vector de posicao pode ser escrito como2

r = rur (6.18)

A primeira derivada em funcao do tempo, da a velocidade v; usando as relacoes6.16 e 6.17, obtemos

v = r ur + r θ uθ (6.19)

e derivando novamente em funcao do tempo, obtemos a aceleracao em coorde-nadas polares.

a =(r − rθ2

)ur +

(rθ + 2rθ

)uθ (6.20)

A unica forca que actua na direccao do versor uθ e a componente tan-gencial do peso, mg sin θ, com sinal negativo. Assim, aplicando a segunda leide Newton para as componentes da forca resultante e da aceleracao na direccaode uθ, obtemos

rθ + 2rθ = −g sin θ (6.21)

Para angulos comprendidos entre −π/2 e π/2, o fio permanecera esticadoe, por tanto, r = l (comprimento do fio) e r = 0. O mesmo resultado sera validopara qualquer valor de θ, se em vez de fio o objecto estiver pendurado de umabarra rıgida que pode rodar livremente. A equacao obtida e

θ = −ω2 sin θ (6.22)

2Estamos a usar letra negrita para representar vectores e distinguı-los de escalares, que saorepresentados em letra italica. Assim, r e um vector, mas r e a sua norma.


onde a frequencia angular esta definida por

ω =√

g

l(6.23)

Se o angulo for pequeno, sin θ pode ser aproximado por θ (em radianos)e obtem-se a equacao de um oscilador harmonico simples. Em geral, a equacaonao e linear e as oscilacoes do angulo θ nao sao funcoes harmonicas simples.No espaco de fase (θ, Ω), as equacoes de evolucao sao:

θ = Ω (6.24)

Ω = −ω2 sin θ (6.25)

2

2

dx/dt = y , dy/dt = -sin(x)

0

0

5

5

-2.5 -2.5

0 0

2.5 2.5


2 10

2

dy/dt=-sin(x), dx/dt=y

5

5

10

10

15

15

-2.5 -2.5

0 0

2.5 2.5


θΩ

= 150°θmax

Figura 6.7: Algumas trajectorias do pendulo simples, no espaco de fase e nodomınio do tempo.

A figura 6.7 mostra algumas trajectorias no espaco de fase, para umpendulo que possa rodar livremente mantendo o comprimento l constante. Os-cilacoes com amplitudes nao muito grandes dao resultados parecidos ao oscila-dor harmonico simples: elipses no espaco de fase e funcao sinusoidal para θ(t).Inclusivamente para amplitudes maiores, como o caso θmax = 150, θ(t) pareceser uma funcao sinusoidal; mas a velocidade angular, Ω(t), ja e claramente umafuncao nao harmonica, como se pode apreciar no lado direito da figura.

Se a velocidade angular inicial for muito elevada, o pendulo pode con-tinuar sempre as voltas, num sentido ou no outro (as duas curvas abertas nocampo de direccoes 6.7). Os dois pontos singulares (onde ha curvas que se cru-zam) que aparecem no diagrama do espaco de fase, correspondem a θ0 = ±180,Ω0 = 0.


6.5 Osciladores acoplados

Vamos estudar nesta seccao as oscilacoes de dois sistemas oscilantes aco-plados. Um exemplo e um sistema de duas molas, com a mesma constanteelastica k, ligadas a duas esferas de massa m, como se mostra na figura 6.8.

0

0

x1

x2

Figura 6.8: Dois osciladores acoplados.

As coordenadas x1 e x2 medem-se a partir da posicao de equilıbrio decada esfera e vao representar a posicao de cada esfera, em funcao do tempo.Assim, o alongamento da mola de cima, a partir da posicao de equilıbrio, seraigual a x1, enquanto que o alongamento da mola de baixo sera x2 − x1. Naesfera de cima actuam as forcas das duas molas, enquanto que na esfera debaixo so actua a forca da segunda mola. As equacoes de movimento das duasesferas sao

mx1 = −kx1 + k(x2 − x1) (6.26)

mx2 = −k(x2 − x1) (6.27)

ou, em forma matrizial, x = Mx, onde x e M sao as matrizes

x =(

x1

x2

)M =

(−2ω2 ω2

ω2 −ω2

)(6.28)


com a definicao habitual ω2 = k/m. E possıvel combinar as duas equacoes eredefinir duas novas variaveis y1 e y2 de forma a que na primeira equacao naoapareca y2 e na segunda nao apareca y1.

Dito de outra forma, uma rotacao do sistema de coordenadas (x1, x2)corresponde a multiplicar x por uma matriz ortogonal A, que resulta na ma-triz y = Ax com coordenadas y1 e y2. Multiplicando os dois lados do sis-tema de equacoes diferenciais por A, se transforma em y = AMx, ou ainday = AMA−1y. Existe sempre uma matriz A que faz com que AMA−1 sejadiagonal:

AMA−1 =(

λ1 00 λ2

)(6.29)

e o sistema de equacoes y = AMA−1y sera simplesmente

y1 = λ1y1 (6.30)

y2 = λ2y2 (6.31)

que sao duas equacoes independentes que podem ser resolvidas em forma in-dependente.

As constantes λ1 e λ2 designam-se de valores proprios da matriz M.Existem vectores proprios v1 e v2 que verificam as equacoes

Mv1 = λ1v1 (6.32)

Mv2 = λ2v2 (6.33)

Voltando ao nosso caso das duas molas, os valores e vectores proprios damatriz M (equacao 6.28) podem ser calculados usando o comando eigenvectorsdo Maxima:

(C8) m: matrix([-2*w^2, w^2], [w^2, -w^2]);

(D8) (−2w2 w2

w2 −w2

)

(C9) eigenvectors(m);

(D9)[[[−(√

5 + 3)w2

2,

(√5− 3

)w2

2

], [1, 1]

],

[1,−

√5− 12

],

[1,

√5 + 12

]]


A primeira lista mostra os dois valores proprios, −(√

5 + 3)ω/2 e (√

5−3)ω/2. A segunda lista indica que cada um desses valores proprios aparece umavez na diagonal da matriz diagonal. Finalmente ha outras duas listas que saoas coordenadas dos dois vectores proprios. Aparece aqui o numero φ (relacaoaurea) que ja tinha sido introduzido no problema 1 do capıtulo 2. Em funcaodesse numero, os dois vectores proprios sao

(C10) v1: [1, 1-%phi];

(D10)

[1, 1− φ]

(C11) v2: [1, %phi];

(D11)

[1, φ]

Cada linha da matriz de transformacao, A, sera um dos vectores propriosapos terem sido divididos pela sua norma, para ficarem vectores unitarios. EmMaxima o produto escalar e representado por um ponto; os vectores normali-zados serao:

(C12) v1: v1/sqrt(v1.v1)$

(C13) v2: v2/sqrt(v2.v2)$

E a matriz de transformacao e

(C14) A: matrix(v1,v2)$

(C15) A, numer;

(D15) (0.85065080835204 −0.525731112119130.52573111211913 0.85065080835204

)

Podemos conferir que AMA−1 e de facto uma matriz diagonal. A matrizinversa pode ser obtida com o comando invert


(C16) Diag: A.m.invert(A)$

(C17) map(ratsimp, Diag);

(D17) ((−φ− 1) w2 0

0 (φ− 2) w2

)

Em funcao das variaveis y1 o sistema de equacoes e

y1 = −(1 + φ)ω2 y1 (6.34)

y2 = −(2− φ)ω2 y2 (6.35)

que sao dois osciladores harmonicos simples, com frequencias angulares φω e(φ − 1)ω. As condicoes iniciais para (y1, y2) obtem-se a partir das condicoesiniciais para (x1, x2), multiplicando pela matriz de transformacao

(C18) A.[x1, x2], numer;

(D18) (0.85065080835204x1 − 0.52573111211913x2

0.85065080835204x2 + 0.52573111211913x1

)

Se as duas esferas sao largadas, a partir do repouso desde posicoes coma relacao x1/x2 = 0.526/0.850, y1 sera nula, e so existira o segundo modo deoscilacao, onde y2 oscila com frequencia angular ω2 = 0.38 ω. Se a relacaoinicial entre os deslocamentos for x1/x2 = −0.850/0.526, o segundo modoy2 anula-se, ficando unicamente o primeiro modo normal y1, com frequenciaangular ω1 = 2, 62 ω. Para condicoes iniciais mais gerais, o movimento dosistema ser’a a sobreposicao dos dois modos normais.

6.6 Problemas

1. Uma esfera de 3 gramas esta sujeita a uma forca resultante, F , quedepende da posicao da esfera, x, de acordo com a equacao F = x3 −6x2 + 3x + 10, onde F e medida em newtons e x em metros.


(a) Encontre a posicao x dos pontos de equilıbrio (onde a forca resul-tante e nula).

(b) Desenhe um grafico da forca resultante, mostrando os 3 pontos deequilıbrio. Identifique o ponto de equilıbrio estavel; isto e, o pontoonde a forca passa de positiva para negativa e, portanto, apontapara o ponto de equilıbrio.

(c) Obtenha os primeiros termos da serie de Taylor para F , a vol-ta do ponto de equilıbrio estavel. O primeiro termo determina ocomportamento da forca perto do ponto de equilıbrio e, por sersemelhante a forca num oscilador harmonico simples, produzira os-cilacoes harmonicas. Calcule o perıodo de oscilacao da esfera, pertodo ponto de equilıbrio estavel.

(d) Qual acha que sera o valor maximo da amplitude de oscilacao daesfera a volta do ponto de equilıbrio estavel?

2. Considere o oscilador harmonico amortecido com equacao de movimento:

2x + ax + 3x = 0

onde a e a constante de amortecimento.

(a) Desenhe os graficos do campo de direccoes e de x(t) no caso a = 3,com condicoes iniciais x(0) = 4, x(0) = −1.

(b) Repita a alınea anterior para a = 6

3. A equacao de movimento de um oscilador forcado em particular e:

11x + 40x + 891x = 1555 cos(at− π

4

)− 900 sin(bt)

onde a = 11 e b = 9. Encontre a solucao da equacao com condicoesiniciais x(0) = x(0) = 0 (as constantes a e b so devem ser substituidasdepois de resolver a equacao, para evitar um erro numerico do programaode2 de Maxima). Desenhe o grafico da solucao entre t = 0 e t = 10 (aflutuacao da amplitude e chamada de batimento ou pulsacao).

4. A equacao diferencial do pendulo simples

θ +g

lsin θ = 0

e uma equacao autonoma. Usando o metodo descrito na seccao 4.7.1,substituia θ por

θ = ΩdΩdθ


onde Ω = θ e a velocidade angular.

(a) Encontre a solucao geral da equacao diferencial obtida para Ω emfuncao de θ

(b) Calcule a solucao particular para as condicoes iniciais θ = θ0 eΩ = 0, para t = 0.

(c) A solucao particular obtida representa a trajectoria no espaco de fa-se (θ, Ω). Desenhe a trajectoria usando os seguintes valores numericos:g = 9.81 m/s2, l = 30 cm, θ0 = 70

5. As equacoes de um sistema de osciladores acoplados particular sao:

x1 = −1000x1 + 500x2

x2 = 500x1 − 1000x2

(a) Encontre os valores e vectores proprios da matriz do sistema.

(b) Escreva as equacoes diferenciais dos dois modos normais de osci-lacao, y1 e y2, e encontre a solucao geral de cada equacao.

(c) Defina a matriz de transformacao onde cada linha e um dos vectoresproprios, normalizados.

(d) Com as condicoes iniciais x1 = 2, x2 = 1, x1 = x2 = 0, em t = 0,encontre os valores iniciais de y1, y2, y1 e y2, e encontre as solucoesparticulares para os modos normais y1, y2 em funcao do tempo.

(e) As variaveis iniciais x1 e x2 obtem-se a partir dos modos normais y1,y2 multiplicando pela matriz de transformacao inversa. Encontre asolucao particular x1(t), x2(t).

6. Uma mola esta suspensa na vertical. Pendura-se um objecto de massa m,que estica a mola uma distancia x0 ate uma nova posicao de equilıbrio.A partir da posicao de equilıbrio, um deslocamento adicional do objectoproduz um movimento harmonico simples. Demonstre que o perıodo deoscilacao e

P = 2π√

x0

g

onde g e a aceleracao da gravidade. Pode concluir que o perıodo deoscilacao e independente da massa m do objecto?


7. Ondas e difusao

Uma onda e uma perturbacao de um meio, que se propaga nesse meio.Por exemplo, o som que ouvimos e uma perturbacao na pressao do ar, que sedesloca a uma velocidade propria do estado do ar.

Vamos limitar o nosso estudo ao caso mais simples de ondas em umadimensao. Por exemplo, uma coluna cilındrica cheia de gas (7.1) onde umaexplosao do gas no meio da coluna faz aparecer um aumento da pressao. Ografico de z(x) representa o valor da pressao, em funcao de x.

x

x

z

Figura 7.1: Perturbacao na pressao do gas num cilindro.

A perturbacao na pressao deslocar-se-a em direccao das duas extremi-dades. Se admitirmos que o perfil de pressao permanece constante, a onda quese propaga no sentido positivo de x, com velocidade v, sera a mesma funcaoinicial z(x), com x − vt a substituir a x, nomeadamente, z(x − vt). Vejamosisto com uma funcao concreta. A funcao do grafico 7.1 e

Ondas e difusao 85

(C1) z: 5*exp(-x^2/4) + 0.2;

(D1)

5e−x24 + 0.2

Vamos substituir z por x − 2t e desenhar um grafico em 3 dimensoes,representando z em funcao de t e de x.

(C2) u: ev(z, x=x-2*t);

(D2)

5e−(x−2t)2

4 + 0.2

(C3) plot3d(u, [t,0,5], [x,-9,9]);

x

z

t

x

t

0

Figura 7.2: Onda de presao z(t, x) com velocidade no sentido positivo de x.

7.1 Problemas

1. Ondas estacionarias. Considere uma corda esticada entre dois pontosa uma distancia L (por exemplo uma corda num violino). Se z representao desvio da corda, desde a sua posicao de equilıbrio, num ponto a uma


distancia x de uma das extremidades e no instante t, z verifica a equacaode onda:

∂2z

∂t2= v2 ∂2z

∂x2

onde v e a velocidade de propagacao, que depende da densidade da cordae da tensao.

(a) Demonstre que f(x) cos(ωt) e solucao da equacao de onda, se afuncao f verifica a equacao do oscilador harmonico simples:

d2f

dx2= −ω2

v2f

Encontre a solucao geral para a funcao f(x).

(b) A solucao particular de f(x) devera verificar as condicoes frontei-ra f(0) = 0 e f(L) = 0 (as extremidades da corda estao fixas).Demonstre que a condicao f(0) = 0 implica

f(x) = A sin(ωx

v

)e que a condicao f(L) = 0 implica que a frequencia angular ω

devera ter algum dos valores:

ω =nπv

Ln = 0, 1, 2, 3, . . .

o numero n− 1 representa o numero de nos na corda.

(c) Escreva a solucao para z(t, x) do segundo harmonico (n = 2).

(d) Faca um grafico em 3 dimensoes da funcao z(t, x) na alınea anterior,usando os seguintes parametros: v = 3, L = 4, A = 5, 0 ≤ t ≤ 1,0 ≤ x ≤ 4. Rode o grafico ate a posicao onde o eixo dos x esta nahorizontal e o eixo dos z na vertical (use o rato, ou modifique osvalores de azimute e elevacao no menu de configuracao). Grave ografico num ficheiro PostScript e transforme-lo num ficheiro PDF.

2. Equacao de difusao. Demonstre que a funcao

u(t, x) =1√

1 + 4kte−

x21+4kt

e solucao da equacao de difusao:

∂u

∂t= k

∂2u

∂x2

Desenhe o grafico de u(t, x), no domınio 0 ≤ t ≤ 10, −10 ≤ x ≤ 10, epara k = 2. Rode o grafico para observar o comportamento da funcaodesde diferentes direccoes.

Ondas e difusao 87

8. Estabilidade e bifurcacao

Os sistemas dinamicos podem apresentar pontos fixos, isto e, pontosno espaco de fase onde o sistema permanece sempre no mesmo estado. Paraidentificar os pontos fixos e estudar o comportamento dos sistemas dinamicosna vizinhanca dos pontos fixos, vamos comecar por estudar sistemas mecanicossimples em equilıbrio estatico.

Um primeiro sistema simples, ja estudado por Euler no seculo XVIII,e uma barra flexıvel, por exemplo uma regua plastica, que suporta um pesoP . Se P nao ultrapassar um valor crıtico Pc, a barra permanece directa e emequilıbrio. Se o peso P ultrapassar o valor crıtico Pc, a regua encurva-se, ateficar numa nova posicao de equilıbrio em que o centro da regua esta afastadouma distancia ∆x da vertical. Acontece que o desvıo da barra pode ser paraa direita ou para a esquerda da vertical. Nomeadamente, existem dois pontosde equilıbrio com ∆x positiva ou negativa.

Em funcao de P , o ponto de equilıbrio ∆x = 0, para P < Pc, separa-seem dois pontos de equilıbrio, ∆x > 0 e ∆x < 0, para P > Pc. Este fenomenoe designado de bifurcacao.

8.1 Teoria linear da estabilidade

Outro sistema mecanico que apresenta o fenomeno de bifurcacao do pon-to de equilıbrio, que vamos estudar com maior pormenor, consiste num pendulosimples ligado a uma mola elastica (figura 8.1).

Um objecto de massa m esta pendurado de uma barra rıgida, de com-primento l, com massa desprezavel. O pendulo pode descrever um cırculocompleto de raio l. Uma mola elastica, de comprimento normal l, e ligada aoobjecto desde o ponto mais alto do cırculo.


m

l

d

l

mg

T

Feθ

α

α

Figura 8.1: Pendulo simples ligado a uma mola e diagrama de corpo livre.

Usaremos coordenadas polares, com origem no centro do cırculo e angulomedido a partir do ponto mais baixo do cırculo. A coordenada r do objecto demassa m e igual a l e permanece constante (r = r = 0). Assim, a segunda leide Newton para as componentes polares da forca resultante e:

Fθ = mld2θ

dt2(8.1)

Fr = −ml

(dθ

dt

)2

(8.2)

Na direccao do versor transversal uθ actuam as componentes da forcaelastica, Fe, e do peso1 A mola faz um angulo α com a barra do pendulo e efacil ver na figura 8.1 que α = θ/2 porque o triangulo formado pela barra e amola e isosceles; assim, a forca resultante na direccao transversal e

Fθ = Fe sinθ

2−mg sin θ (8.3)

Para calcular a forca elastica, precisamos do alongamento da mola. Emfuncao do angulo θ, o comprimento d da mola verifica a relacao

d

2= l cos

θ

2(8.4)

1A equacao da forca Fr servira para calcular a tensao no barra, depois de θ(t) ter sidocalculada, mas nao vamos abordar essa parte do problema.

Estabilidade e bifurcacao 89

O alongamento da mola e igual a d− l e a forca elastica obtem-se multiplicandoo alongamento pela constante elastica k

Fe = kl

(2 cos

θ

2− 1)

(8.5)

Subsituindo na equacao 8.1 obtemos uma equacao diferencial ordinariade segunda ordem

d2θ

dt2=

k

msin

θ

2

(2 cos

θ

2− 1)− g

lsin θ (8.6)

As constantes k/m e g/l tem unidades de tempo ao quadrado. Convemescrever a equacao numa forma independente de unidades fısicas; em funcaodo tempo “adimensional”,

u =√

g

lt (8.7)

e substituindo sin θ = 2 sin(θ/2) cos(θ/2), a equacao fica

θ =[2(S − 1) cos

θ

2− S

]sin

θ

2(8.8)

onde θ e a segunda derivada do angulo em funcao de u, e o parametro numericoS e relacao entre a forca elastica e o peso, na posicao mais baixa:

S =kl

mg(8.9)

O sistema estara em equilıbrio nos pontos onde a forca Fθ for nula.Nomeadamente, nos valores de θ que fazem com que o lado direito da equacao8.8 seja nulo. Existem duas possibilidades para que isso aconteca: sin(θ/2) = 0ou a expressao dentro dos parentesis quadrados em 8.8 igual a zero.

O segundo caso sera estudado na proxima seccao. No primeiro caso,sin(θ/2) = 0, o ponto de equilıbrio encontra-se em θ = 0, isto e, na posicaomais baixa do pendulo2. A serie de Taylor da funcao no lado direito de 8.8, apartir de θ = 0 e (

S

2− 1)

θ −(

78S − 1

6

)θ3 + . . . (8.10)

Mantendo unicamente o termo dominante, a equacao diferencial 8.8aproxima-se a

θ =(

S

2− 1)

θ (8.11)

2Neste sistema mecanico o valor de θ nao pode aumentar ate π, pois isso implicaria comprimira mola ate um comprimento nulo, o qual e impossıvel.


As forma das solucoes dessa equacao serao diferentes para valores de S maioresou menores que 2. No caso de uma constante elastica fraca, S < 2, a equacao8.11 e a equacao de um oscilador harmonico simples:

θ = −p2θ p =

√1− S

2(8.12)

As solucoes serao movimentos oscilatorios, a volta de ponto de equilıbrio θ = 0

θ = C1 cos(pu) + C2 sin(pu) (8.13)

dizemos que o ponto θ = 0 e um ponto de equilıbrio estavel. Quando osistema estiver perto do ponto de equilıbrio, a sua tendencia sera regressarpara o ponto de equilıbrio, devido a accao do peso que neste caso e maior doque a forca elastica.

No caso de uma constante elastica forte, S > 2, a equacao 8.11 ja nao ea equacao de um oscilador harmonico simples

θ = p2θ p =

√S

2− 1 (8.14)

As solucoes dessa equacao sao funcoes exponenciais

θ =(

θ0

2+

ω0

2p

)epu +

(θ0

2− ω0

2p

)e−pu (8.15)

o primeiro termo cresce rapidamente e implica que o pendulo se afasta daposicao de equilıbrio3. Trata-se de um ponto de equilıbrio instavel; quandoo sistema e afastado ligeiramente da posicao de equilıbrio θ = 0, a forca elasticapuxa o pendulo para cima; repare que a forca elastica tera que ser pelo menoso dobro do peso, devido a diferenca dos angulos na equacao 8.3.

8.2 Bifurcacao

Consideremos agora o segundo ponto de equilıbrio no sistema estudadona seccao anterior, definido pela equacao

2(S − 1) cosθ

2− S = 0 (8.16)

3quando θ ja nao estiver perto de zero, a aproximacao que fizemos para obter a solucaoanterior ja nao e valida.


a solucao dessa equacao da dois valores de θ, com o mesmo valor absoluto:

θ = ±2 cos−1

[S

2(S − 1)

]se S > 2 (8.17)

Se a constante elastica nao for suficientemente forte, S < 2, so existe o primeiroponto de equilıbrio θ = 0, e esse ponto de equilıbrio e estavel. Se S for maior que2, o ponto de equilıbrio estavel bifurca-se em dois angulos, positivo e negativo,com o mesmo valor absoluto (a demonstracao de que esses pontos de equilıbriosao estaveis, deixa-se ao leitor como um dos problemas propostos no fim docapıtulo).

S

θ

2

2π/3

−2π/3

instável

estável

Figura 8.2: Diagrama de bifurcacao do ponto de equilıbrio do pendulo commola e energia potencial.

Em S > 2, o ponto de equilıbrio θ = 0 continua a existir, mas e agorainstavel. No limite S → ∞, os pontos de equilıbrio situam-se em ±2π/3,nomeadamente, os pontos onde a mola tem o seu comprimento normal l. Afigura 8.2 mostra a posicao do ponto de equilıbrio estavel (curva contınua) edo ponto de equilıbrio instavel (curva tracejada) em funcao do parametro S.Os pontos de equilıbrio estavel correspondem a os “vales” na funcao de energiapotencial, definida como menos a primitiva do lado direito de 8.8.


8.3 Pontos fixos

Consideremos um sistema autonomo com duas equacoes diferenciais deprimeira ordem

x1 = f1(x1, x2) (8.18)

x2 = f2(x1, x2) (8.19)

um ponto fixo e um ponto (x1, x2) no espaco de fase, que verifica as condicoes

f1(x1, x2) = 0 (8.20)

f2(x1, x2) = 0 (8.21)

para um ponto qualquer (c, d), podemos escrever cada uma das funcoes f1 e f2

na forma de uma serie de Taylor; por exemplo, para f1

f1 = k0 +k1(x1−c)+k2(x2−d)+k3(x1−c)(x2−d)+k4(x1−c)2 + . . . (8.22)

se (c, d) for um ponto fixo, a constante k0 sera nula. Os termos dominantes naserie de Taylor sao os termos com as constantes k1 e k2, e essas constantes saoiguais as derivadas parciais da funcao, no ponto (c, d); assim, perto do pontofixo (c, d) uma boa aproximacao para a funcao f1 e:

f1(x1, x2) =∂f1

∂x1

∣∣∣∣(c,d)

(x1 − c) +∂f1

∂x2

∣∣∣∣(c,d)

(x2 − d) (8.23)

o mesmo tipo de aproximacao pode ser feito para f2. Assim, perto do pontofixo (c, d), o sistema de equacoes diferenciais pode ser aproximado pelo sistemalinear [

x1

x2

]=

∂f1

∂x1

∂f1

∂x2

∂f2

∂x1

∂f2

∂x2

(c,d)

[x1 − c

x2 − d

](8.24)

A matriz das derivadas parciais das duas funcoes f1 e f2 designa-se de matrizJacobiana

J(x1, x2) =

∂f1

∂x1

∂f1

∂x2

∂f2

∂x1

∂f2

∂x2

(8.25)


Uma mudanca de coordenadas, y1 = x1 − c, y2 = x2 − d, que equivale adeslocar a origem para o ponto fixo, torna o sistema linear homogeneo

y = J(c, d)y (8.26)

onde y e o vector (matriz coluna) com dois coordenadas y1 e y2. Se numinstante dado o estado do sistema for o vector y0, o produto matrizial Jy serao deslocamento do estado no intervalo infinitesimal dt. Se y0 for um vectorproprio da matriz J

Jy0 = λy0 (8.27)

onde λ e o respectivo valor proprio. Isto implica que o deslocamento do estadoe na mesma direccao do vector proprio y0; o estado do sistema desloca-sena direccao do vector proprio, afastando-se da origem, se o valor proprio forpositivo, ou aproximando-se da origem se o valor proprio for negativo.

Se os dois valores proprios da matriz J forem reais e diferentes, existemdois vectores proprios diferentes. Se os dois valores proprios sao positivos, oponto fixo e um ponto onde saem linhas de campo em todas as direccoes; oponto fixo e repulsivo (figura 8.3). Se os dois valores proprios forem negativos,no ponto fixo entram linhas de campo em todas as direccoes; o ponto fixo eatractivo. Se um dos valores proprios for positivo e o outro negativo, ha duaslinhas de campo que entram no ponto (na direccao do vector proprio do valorproprio negativo) e duas linhas de campo que saem do ponto (na direccao dovector proprio do valor proprio positivo); o ponto fixo e um ponto de sela.

Figura 8.3: Os tres tipos de pontos fixos, quando os valores proprios sao reais.Ponto repulsivo, atractivo e de sela.

As coordenadas (y1, y2) podem ser transformadas em outras duas coor-


denadas normais (z1, z2) onde o sistema de equacoes e mais simples[z1

z2

]=

[λ1 00 λ2

][z1

z2

](8.28)

a solucao geral desse sistema e

z1 = C1eλ1t (8.29)

z2 = C2eλ2t (8.30)

Os dois valores proprios podem tambem ser dois numeros complexos,cada um deles complexo conjugado do outro:

λ = u± iv (8.31)

nesse caso as duas solucoes 8.29 e 8.30 sao

z1 = C1eut [cos(vt) + i sin(vt)] (8.32)

z2 = C2eut [cos(vt)− i sin(vt)] (8.33)

as constantes complexas C1 e C2 deverao ser escolhidas de forma a dar valoresreais para as funcoes y1 e y2; assim, y1 e y2 serao uma combinacao linear dasfuncoes eut cos(vt), eut sin(vt).

Figura 8.4: Os tres tipos de ponto fixo, quando os valores proprios sao com-plexos.

Em funcao dos valores de u e v, podemos classificar 3 casos para as linhasde campo (figura 8.4): se u for nula (valores proprios imaginarios) as linhas decampo serao elipses a volta do ponto fixo; diz-se que o ponto fixo e um centro.Se u for positiva, as linhas de campo serao espirais que se afastam do ponto


fixo; o ponto fixo e um centro repulsivo. Finalmente, se u for negativa, aslinhas de campo serao espirais que entram no ponto fixo; o ponto e um centroatractivo, ou simplesmente um atractor.

Pode tambem existir um unico valor proprio com multiplicidade igual a2. Nesse caso, se existirem dois vectores proprios linearmente independentes,entao qualquer outro vector e vector proprio com o mesmo valor proprio; aslinhas de campo serao rectas a entrar ou sair do ponto fixo, conforme o valorproprio for negativo ou positivo. Se so existe um vector proprio linearmenteindependente, as linhas de campo agrupam-se na direccao do vector proprio.

8.4 Problemas

1. Para demonstrar que os pontos de equilıbrio estudados na seccao 8.2 saopontos de equilıbrio estavel, e preciso mostrar que a derivada da funcaono lado direito da equacao 8.8 e negativa nos pontos de equilıbrio. Deriveo lado direito da equacao 8.8, em funcao de θ, e demonstre que a condicao8.16 implica um valor negativo para a derivada.

2. O sistema de equacoes:

x = x(2− y)

y =y

2(x− 3)

define um modelo de Lotka-Volterra.

(a) Encontre os pontos fixos.

(b) Calcule a matriz Jacobiana do sistema.

(c) Para cada ponto fixo, substituia as coordenadas do ponto na matrizJacobiana e encontre os valores proprios e vectores proprios. Apartir do resultado, classifique cada um dos pontos fixos e diga comoserao as linhas do campo de direccoes na vizinhanca do ponto.

(d) Use plotdf para visualizar os resultados da alınea anterior. Subs-tituia as coordenadas do ponto fixo para os parametros xcenter eycenter, use valores pequenos de xradius e yradius, por exem-plo 0.5, aumente o valor de nsteps (por exemplo ate 1000) e use“forward” no campo direction.


3. Encontre os pontos fixos do sistema de Verhulst:

x = x− x2 − 2xy

y = y − y2 + 5xy

e diga como sao as linhas do campo de direccoes perto de cada ponto.


9. Osciladores nao-lineares e caos

9.1 Oscilador cubico, forcado

A equacao de movimento do oscilador cubico e

x + kx + x3 = A cos(Ωt) (9.1)

E um caso particular da equacao de Duffing.Espaco de fase

t = 1 (9.2)

x = y (9.3)

y = −ky − x3 + A cos(Ωt) (9.4)

9.2 Oscilador de van der Pol

Um oscilador nao linear, que aparece em teoria de circuitos electricos, eo oscilador de van der Pol, com equacao diferencial

x + ε(x2 − 1)x + x = 0 (9.5)

onde ε e um parametro.Este sistema apresenta oscilacoes auto-excitadas, isto e, a partir de

um estado inicial em que x e a sua derivada x sao nulas, o sistema comecaa oscilar de forma espontanea, com amplitude crescente ate atingir um valorcaracterıstico do sistema.

Outra caracterıstica deste sistema e que possui ciclos de fronteira: pa-ra diferentes valores do estado inicial, o comportamento assimptotico e sempreo mesmo ciclo, com amplitude e frequencia caracterısticas do sistema.


9.3 Equacoes de Lorenz

O sistema de equacoes, estudado por Lorenz em 1963,

x = a(y − x) (9.6)

y = rx− y − xz (9.7)

z = xy − bz (9.8)

apresenta solucoes caoticas.

Algumas propriedades deste sistema sao:

· Simetria (x, y, z) −→ (−x,−y, z)

· Eixo z invariante

· Se 0 < r < 1, a origem e ponto fixo (atractor).

· O ponto fixo bifurca-se em r = 1

· Para r > 24.06 aparece caos (atractor estranho).

t

x

10 20

20

-20

Figura 9.1: Oscilacoes do sistema de Lorentz para dois valores muito proximosdo valor inicial: x(0) = 5 (vermelho) e x(0) = 5.005 (azul). Parametros:a = 10, b = 8/3, r = 28, y(0) = 5, z(0) = 5.

Osciladores nao-lineares e caos 99

Figura 9.2: Trajectoria no espaco de fase, projectada no plano xz. Osparametros sao os mesmos da figura 9.1, com x(0) = 5.


9.4 Caos

As manifestacoes do comportamento caotico de um sistema sao:

1. Comportamento nao-periodico.

2. Sensibilidade as condicoes iniciais.

3. Estrutura fractal.

9.5 Equacoes de Rossler

Em 1976 Rossler estudou o sistema:

x = −(y + z) (9.9)

y = x + ay (9.10)

z = b + xz − cz (9.11)

Para valores pequenos do parametro c as trajectorias do sistema, atingemum estado estacionario que e um ciclo com perıodo simples.

line0

1

1

-2.5

-2.5

0

0

2.5

2.5

-2.5 -2.5

0 0

2.5 2.5

line0

5 25

2

0

0

25

25

-2.5 -2.5

0 0

2.5 2.5

Figura 9.3: O caso c = 2.3 (com parametros a = b = 0.2) conduz a um ciclo defronteira com perıodo simples.


Para valores maiores do parametro c o ciclo fronteira pode ter perıododuplo, quadruplo, etc. Finalmente, a partir dum certo valor do parametro c osistema e caotico.

line0

2

2

-5

-5

0

0

5

5

-5 -5

0 0

line0

5 25

2

0

0

25

25

-5 -5

0 0

5 5

Figura 9.4: O caso c = 3.3 (com parametros a = b = 0.2) conduz a um ciclo defronteira com perıodo duplo.

9.6 Problemas

1. As solucoes do sistema de equacoes de Rossler

x = −y − z

y = x + 0, 2 y

z = 0, 2 + (x− c)z

no caso c = 3, sao ciclos de fronteira; isto e, depois de passado umtempo suficientemente grande, as variaveis x, y e z descrevem ciclos quese repetem periodicamente.

(a) Use o programa ode rk4 para encontrar a solucao do sistema, comcondicoes iniciais x(0) = z(0) = 0, y(0) = 4, no intervalo 0 ≤ t ≤200; use 5000 passos (h = 0, 04).

(b) Usando unicamente o intervalo 160 ≤ t ≤ 200 da solucao encontra-da na alınea anterior, desenhe os graficos de y em funcao de x e x


em funcao de t.

(c) Calcule o perıodo dos ciclos representados nos graficos da alıneaanterior (observe o conteudo do ficheiro rungekutta.dat e calculeo intervalo de tempo entre dois pontos sucessivos onde x atinge oseu valor maximo).

2. O sistema de Lorenz

x = 10(y − x)

y = 28x− y − z x

z = −83z + x y

tem solucoes caoticas.

(a) Encontre a solucao do sistema usando o programa ode rk4, comos parametros: 0 ≤ t ≤ 10, x(0) = y(0) = z(0) = 5. Comecepor usar 500 passos (h = 0, 02) e escreva os valores finais obtidospara (x, y, z) em t = 10. Repita o procedimento com 1000 passose compare os valores de (x, y, z) em t = 10; continue a duplicar onumero de passos ate que os valores obtidos coincidam no primeiroalgarismo decimal.

(b) Desenhe as projeccoes da solucao obtida, nos planos xy e xz doespaco de fase.

(c) Demonstre que o sistema tem 3 pontos fixos, um deles na origem,e que a matriz Jacobiana tem valores proprios complexos nos doispontos fixos diferentes da origem

3. A equacao de um pendulo forcado (em unidades adimensionais) e

x + cx + sinx = F cos(ωt)

onde x e o angulo em radianos, c e o coeficiente de amortecimento, e F

e ω sao os parametros da forca externa. Escreva a equacao na forma deum sistema autonomo. Estude as solucoes para os parametros c = 0, 05,F = 0, 7 e ω = 0, 1; 0, 2; 0, 3; . . . ; 1.5. Diga para que valores de ω asolucao e caotica.


10. Sistemas dinamicos discretos

Um sistema dinamico discreto, e um sistema para o qual queremos sa-ber o seu estado durante uma sequencia discreta de tempos t0, t1, t2, . . ..Provavelmente o estado do sistema muda em forma discreta no fim de cadaintervalo, mas tambem e possıvel que o estado mude continuamente mas soestamos interessados em saber o estado no inicio de cada intervalo.

Neste capıtulo estudaremos unicamente sistemas discretos em uma di-mensao e no proximo capıtulo abordaremos os sistemas em duas dimensoes. Oestado de um sistema discreto em uma dimensao e determinado completamentepor uma variavel, x. Os valores da variavel de estado nos instantes t0, t1, t2,. . . sera uma sequencia x0, x1, x2, . . . (esta sequencia costuma-se deisgnarde orbita do sistema). O intervalo de tempo entre dois instantes sucessivos tn

e tn+1 nao tem que ser constante.A equacao de evolucao permite calcular o estado num instante a partir

da evolucao do estado nos instantes anteriores. No caso mais simples, a equacaode evolucao sera uma equacao de diferencas de primeira ordem que permitecalcular o estado xn+1 a partir do estado xn no instante imediatamente anterior:

xn+1 = F (xn) (10.1)

onde F (x) e alguma funcao conhecida. Dado um estado inicial x0, aplicacoessucessivas de funcao F permitem obter facilmente a sequencia de estados xn.Em alguns casos pode ser possıvel obter uma expressao geral para xn em funcaode n.

Exemplo 10.1Durante um perıodo de tempo ∆t, um atomo de uma substancia radioactivatem uma probabilidade k de se desintegrar. Calcule a evolucao do numero deatomos.

A variavel de estado neste caso e o numero de atomos num instante.Este problema pode ser estudado como sistema contınuo, por meio de uma


equacao diferencial simples, mas se nos limitarmos a calcular o numero deatomos unicamente no fim de cada perıodo ∆t, obtemos um sistema discreto.

Seja Ni o numero de atomos num instante ti. Durante um intervalo∆t, cada atomo tem uma probabilidade k de se desintegrar (k devera ser umnumero entre 0 e 1). Assim, vao se desintegrar em media kNi atomos duranteesse intervalo. No instante ti+1 = ti + ∆t o numero de atomos sera

Ni+1 = Ni − kNi (10.2)

Essa relacao de recorrencia permite-nos calcular o numero de atomosem instantes posteriores:

N1 = (1− k)N0 (10.3)

N2 = (1− k)N1 = (1− k)2N0 (10.4)... (10.5)

Ni = (1− k)iN0 (10.6)

10.1 Resolucao Numerica de Equacoes

Outros exemplos de sistemas dinamicos e relacoes de recorrencia sao osmetodos numericos usados para resolver equacoes. No metodo de Newton,para encontrar uma raiz1 de uma funcao contınua e derivavel f(x), comecamospor admitir que a raiz se encontra num ponto x0 e melhoramos a nossa apro-ximacao encontrando o ponto x1 onde a tangente em f(x0) corta o eixo dos x

(ver figura 10.1)

x1 = x0 −f(x0)f ′(x0)

(10.7)

Podemos usar a mesma equacao para calcular uma outra aproximacaox2 a partir de x1. Em geral

xn+1 = xn −f(xn)f ′(xn)

(10.8)

1Por raızes de uma funcao entendemos os pontos onde a funcao seja nula

Sistemas dinamicos discretos 105

x

f

10 20

20

-20

x0x1

f(x0)

Figura 10.1: Metodo de Newton para aproximacao a uma raız.

Exemplo 10.2Calculo da raiz quadrada de 5 usando unicamente produtos divisoes e adicoes.

A raiz quadrada de 5 e uma das solucoes da equacao x2 = 5. Assim,para encontrar a raiz de 5 podemos procurar a raiz positiva da funcao x2 − 5.A derivada da funcao x2

n − 5 e

f ′(xn) = 2xn (10.9)

e substituindo na relacao de recorrencia para o metodo de Newton obtemos

xn+1 = xn −x2

n − 52xn

=12

(xn +

5xn

)(10.10)

Usando uma calculadora ou um programa de computador, podemos ite-rar a relacao de recorrencia ate atingir um valor limite. Em Maxima, teremosque ter cuidado de obter cada valor da sequencia em forma numerica, para naoterminar com expressoes recurssivas complicadas:

(C1) x : 1$

(C2) for i thru 7 do (x: float((x + 5/x)/2), print(x))$

3.0

2.333333333333334

2.238095238095238

2.236068895643363

2.236067977499978

2.23606797749979

2.23606797749979


Exemplo 10.3Encontre a solucao da equacao x = cos x

Neste caso a equacao ja esta numa forma que permite escrever imedia-tamente uma relacao de recorrencia

xn+1 = cos(xn) (10.11)

Comecando com um valor inicial qualquer, calculamos a sequencia atechegar a um limite

(C3) x : 1$

(C4) for i thru 15 do (x: float(cos(x)), print(x))$

0.54030230586814

0.85755321584639

0.65428979049778

0.79348035874257

0.70136877362276

0.76395968290065

0.72210242502671

0.75041776176376

0.73140404242251

0.74423735490056

0.73560474043635

0.74142508661011

0.73750689051324

0.74014733556788

Este metodo designa-se de metodo de iteracao. Tem a vantagem denao ser preciso calcular derivadas, e o calculo e muito simples; numa calculadorabasta com escrever um valor inicial de x e carregar repetidamente na tecla doco-seno ate atingir um valor que nao muda. Mas a convergencia usando estemetodo e lenta. O metodo de Newton pode ser usado para obter uma sequenciaque converge mais rapidamente.


10.2 Representacao grafica

Uma forma grafica de representar a evolucao do sistema consiste emdesenhar um ponto para cada instante, com abcissa igual ao ındice n e ordenadaigual a xn. Em Maxima, o programa que se segue permite desenhar esse tipode diagrama

(C5) evolution(f, x0, n) :=

block

([x:x0, xlist:[0, x0], numer: true],

for i thru n do

(x: ev(f),

xlist: append(xlist, [i, x])),

openplot_curves

([["plotpoints 1 nolines"], xlist])

)$

Deverao ser dados 3 argumentos a o programa; uma expressao com umaunica variavel, que devera ser x (a funcao F (x) no lado direito da equacao dediferencas 10.1), o valor inicial x0 e o numero de elementos da sequencia quedeverao ser representados.

Por exemplo, no caso F (x) = cos x (exemplo 10.3), com valor inicialx0 = 2 e vinte iteracoes, o comando usado e evolution(cos(x),2,20) e oresultado aparece na figura 10.2.

n

xn

5 10 15

1

2

-0.5

Figura 10.2: Solucao de xn+1 = cos(xn) com x0 = 2.


Outro tipo de diagrama que sera muito util para a analisar os siste-mas dinamicos discretos em uma dimensao e o diagrama de degraus, 2 queconsiste em representar as funcoes y = f(x), y = x e uma serie alternada desegmentos verticais e horizontais que unem os pontos (x0,x0), (x0,x1), (x1, x1),(x1,x2), etc. Por exemplo, a figura 10.3 mostra o diagrama de degraus no casoda sequencia representada na figura 10.2

x

y

1-1

1

2

-1

Figura 10.3: Diagrama de degraus para xn+1 = cos(xn) com x0 = 2.

Em Maxima, o grafico de degraus pode ser produzido com um programasimples:

(C6) staircase(f, x0, n) :=

block

([xf, x:x0, xmax:x0, xmin:x0, xm, h, xfun:[],

stair:[x0,x0], numer: true],

for i thru n do

(xf: ev(f),

stair: append(stair, [x, xf, xf, xf]),

x: xf,

if x < xmin then xmin: x,

if x > xmax then xmax: x),

h: (xmax - xmin)/200,

x: xmin,

for i thru 200 do

2Em ingles staircase diagram ou cobweb diagram.


(xfun: append(xfun, [x, ev(f)]), x: x + h),

xm: (xmax + xmin)/2,

openplot_curves([stair, [xmin,xmin,xm,xm,xmax,xmax], xfun])

)$

O grafico 10.3 foi obtido com staircase(cos(x),2,8). O diagrama dedegraus permite-nos compreender quando uma sequencia diverge ou converge equal o valor para onde converge. Por exemplo, consideremos o sistema xn+1 =x2

n − 0.2. Se comecarmos com um valor x0 = 1.1 obtem-se o grafico no ladoesquerdo da figura 10.4; vemos que a sequencia converge para um valor x

negativo que e o ponto de interseccao entre as funcoes y = x2 − 0.2 e y = x,nomeadamente, x = (5− 3

√5)/10.

As duas funcoes interceptam-se num outro ponto positivo x = (5 −3√

5)/10. No grafico podemos observar que a pesar de que o valor inicialestava muito perto do segundo ponto de interseccao, a sequencia afasta-se parao primeiro ponto, devido a que entre os dois pontos de interseccao a funcaox2 − 0.2 encontra-se por baixo de y = x. Se usarmos um valor inicial a direitado segundo ponto de interseccao, x0 = 1.5, a sequencia cresce rapidamenteafastando-se para infinito (lado direito da figura 10.4. Para que as sequenciasconvergesem para o segundo ponto de interseccao, teria sido necessario queentre os dois pontos f(x) estivesse por cima de y = x; isto e, que o declıve def(x) fosse menor que 1, em vez de maior que 1, no segundo ponto de interseccao.

x

y

1

1

x

y

2 4

2

4

Figura 10.4: Solucao do sistema xn+1 = x2n−0.2 com valor inicial 1.1 (esquerda)

e 1.5 (direita).


Exemplo 10.4Analise as solucoes do modelo logıstico, que consiste em considerar uma popu-lacao P com uma taxa de natalidade constante, a, e uma taxa de mortalidadedirectamente proporcional a populacao, bPn, onde a e b sao constantes.

A populacao em questao pode ser por exemplo um grupo de algumaespecie animal, com ciclo de aparelhamento anual, de maneira que P0, P1,P2, . . . representa o numero de especimes durante varios anos sucessivos. SePn for o numero de especimes no inicio do perıodo n. Durante esse perıodonascem, em media, aPn especimes e morrem bP 2

n . No inıcio do perıodo n + 1a populacao sera

Pn+1 = (a + 1)Pn(1− b

a + 1Pn) (10.12)

Convem definir uma variavel x por meio de xn = b Pn/(a+1). Obtemosassim uma equacao com um unico parametro c = a + 1

xn+1 = c xn(1− xn) (10.13)

A figura 10.5 mostra as orbitas obtidas com um valor inicial x0 = 0.1,nos casos em que c = 2 e c = 4. Para c = 2, a orbita converge rapidamentepara o ponto fixo x = 0.5. Para c = 4 a solucao e caotica.

x

y

1

1

0.5

0.5

x

y

1

1

0.5

0.5

Figura 10.5: Orbitas do modelo logıstico com valor inicial 0.1. Para c = 2(esquerda) a sequencia converge, mas para c = 4 (direita) o comportamento ecaotico.


10.3 Pontos fixos

Um ponto fixo do sistema 10.1 e um ponto c onde o estado do sistemapermanece constante. Para isso acontecer sera necessario e suficiente que

f(c) = c (10.14)

Do ponto de vista grafico, os pontos fixos serao todos os pontos ondea curva y = f(x) intersecta a recta y = x. Por exemplo, no caso do modelologıstico, a figura 10.5 mostra que nos casos c = 2 e c = 4 existem dois pontosfixos, um deles em x = 0. Podemos usar o comando solve do Maxima paraencontrar os pontos fixos; no caso c = 4

(C7) flogistic: 4*x*(1-x)$

(C8) fixos: solve(flogistic - x);

(D8) [x =

34, x = 0

]

Os pontos fixos sao 0 e 0.75.Consideremos um ponto fixo, onde a funcao F (x) intersecta a recta com

declıve igual a 1, tal que o valor absoluto da derivada da funcao, F ′(x), formaior que 1 nesse ponto. Para um ponto que esteja a uma distancia ∆x doponto, F (x) produz um ponto que esta mais afastado do ponto fixo; e se aderivada for positiva o ponto muda de esquerda para direita do ponto fixo, ouviceversa. Por tanto, as orbitas perto do ponto fixo afastar-se-ao do ponto.

Se a derivada for maior que 1 e positiva, a orbita afasta-se do pontoformando uma escada no diagrama de degraus, tal como no lado direito dafigura 10.4; designamos o ponto de no repulsivo. Se a derivada for negativa emenor que -1, as orbitas tambem se afastam, mas neste caso alternando de ladopara lado, formando uma “teia de aranha” no diagrama de degraus, dizemosque o ponto fixo e um centro repulsivo.

Se o valor absoluto da derivada for menor que 1, F (x) faz com que ospontos perto do ponto fixo cheguem mais perto. O ponto fixo sera um noatractivo se a derivada for positiva (um exemplo e o lado esquerdo na figura10.5) ou um centro atractivo se a derivada for negativa (por exemplo o pontofixo na figura 10.3). Resumindo, temos os seguintes tipos de pontos fixos xf :


1. No atractivo, se 0 ≤ F ′(xf ) < 1

2. Ponto fixo neutro, se F ′(xf ) = 1

3. No repulsivo, se F ′(xf ) > 1

4. Centro atractivo, se −1 < F ′(xf ) < 0

5. Centro repulsivo, se F ′(xf ) < −1

6. Centro indiferente, se F ′(xf ) = −1

Voltando ao nosso exemplo do modelo logıstico (ver C7 ate D8acima), ovalor da derivada de F nos pontos fixos e:

(C9) dflogistic: diff(flogistic, x);

(D9)

4 (1− x)− 4x

(C10) dflogistic, fixos[1];

(D10)

−2

(C11) dflogistic, fixos[2];

(D11)

4

Assim, os dois pontos sao, no caso c = 4, repulsivos. Em x = 0 ha umno repulsivo, e em x = 0.75 um centro repulsivo.

10.4 Pontos periodicos

Se a sequencia x0, x1, x2, . . . for uma orbita do sistema dinamico

xn+1 = F (xn) (10.15)


um elemento qualquer na sequencia pode ser obtido directamente a partir dex0, por meio da funcao composta Fn

xn = Fn(x0) = F (F (. . . (F︸︷︷︸n vezes

(x)))) (10.16)

Uma orbita ou ciclo de perıodo 2 e uma sequencia com dois valoresrepetidos: x0, x1, x0, x1, . . .. Os pontos x0, x1 sao pontos periodicos comperıodo igual a 2. Para que isto aconteca, e suficiente com que x2 = x0, istoe, F 2(x0) = x0; como F (x0) = x1 temos tambem que F 2(x1) = x1. Assim, x0

e x1 sao pontos fixos de F 2(x). O ciclo sera atractivo ou repulsivo (estavel ouinstavel) segundo o valor que a derivada de F 2 tiver em cada ponto do ciclo.

Para calcular a derivada de F 2 em x0 usa-se a regra da cadeia

(F 2(x0))′ = (F (F (x0)))′ = F ′(F (x0))F ′(x0) = F ′(x0)F ′(x1) (10.17)

a derivada de F 2 e igual nos dois pontos x0, x1 que fazem parte do ciclo, e eigual ao produto da derivada de F nos dois pontos.

Em forma geral, se Fm(x0) = x0, existe um ciclo de perıodo m, formadopelos conjunto x0, x1, . . ., xm. Todos os pontos no conjunto serao pontos fixosde fm. Se o valor absoluto do produto de F ′(xi) para os m pontos no ciclo formaior que 1, o ciclo sera repulsivo; se o produto for menor que 1, o ciclo seraatractivo, e se o produto for igual a 1, o ciclo sera neutro.

Exemplo 10.5Encontre os ciclos de perıodo igual a 2 do sistema logıstico

xn+1 = 3.1 xn(1− xn)

e diga se sao atractivos, repulsivos ou neutros.

Comecamos por definir a funcao F (x) e o seu quadrado F 2(x)

(C12) flog: 3.1*x*(1-x)$

(C13) flog2: flog, x=flog;

(D13)

9.610000000000001 (1− x)x (1− 3.1 (1− x) x)


Os pontos periodicos de perıodo igual a dois sao as solucoes da equacaoF 2(x)− x = 0

(C14) periodicos: solve(flog2 - x);

(D14) [x = −

√41− 41

62, x =

√41 + 41

62, x =

2131

, x = 0

]

Os dois ultimos pontos, nomeadamente 0 e 21/31 sao tambem pontosfixos (isto ve-se pela solucao de F (x)−x = 0) e cada um constitui, isoladamente,um ciclo “trivial”, pois realmente dao origem a sequencias constantes, quetambem podem ser consideradas ciclos com um perıodo qualquer, mas nao saorealmente ciclos.

Os outros dois pontos fazem parte de um ciclo de perıodo dois, como sepode conferir facilmente: se calcularmos a orbita com valor inicial igual a umdesses pontos, obtem-se uma sequencia que oscila entre os dois valores.

Para saber se o ciclo e atractivo ou repulsivo, calculamos a derivada emcada um dos dois pontos do ciclo

(C15) dflog: diff(flog, x);

(D15)

3.1 (1− x)− 3.1x

(C16) dflog, periodicos[1], ratsimp, numer;

(D16)

−0.3596875787352

(C17) dflog, periodicos[1], ratsimp, numer;

(D17)

−0.3596875787352

O valor absoluto do produto das duas derivadas e menor que 1 e, portanto, o ciclo e atractivo.


10.5 Problemas

1. Considere a sequencia xn+1 = |xn − 2|

(a) Encontre todos os pontos fixos. Mostre os pontos fixos num graficocom as funcoes F (x) = |x− 2| e G(x) = x.

(b) Diga como sera a sequencia se x0 for um numero inteiro par ouımpar.

(c) Encontre a solucao para o valor inicial 8.2.

(d) Encontre todos os pontos periodicos de perıodo igual a dois. Mostreos pontos periodicos num grafico de F 2(x) e G(x).

2. Para cada funcao na lista que se segue, o ponto x = 0 faz parte de umaorbita periodica do sistema xn+1 = F (xn). Diga em cada caso qual e operıodo da orbita e calcule a respectiva derivada para determinar se aorbita e atractiva, repulsiva ou neutra. Desenhe um grafico de degrauspara uma orbita com valor inicial 0 e para valores iniciais perto de 0.

(a) F (x) = 1− x2

(b) F (x) =π

2cos x

(c) F (x) = −12x3 − 3

2x2 + 1

(d) F (x) = |x− 2| − 1

(e) F (x) = − 4π

arctan(x + 1)

3. Desenhe os diagramas de xn em funcao de n, com n entre 0 e 30, para asequencia xn+1 = x2

n− 2, nos casos x0 = 0, x0 = 0.0000001, x0 = 0.0001e x0 = 0.1. Desenhe o grafico de degraus em cada caso, com n = 30 ecom n = 100.

4. Considere a funcao

F (x) =

2x , se x ≤ 14− 2x , se x ≥ 1

(a) Mostre que F (x) e equivalente a 2− 2|x− 1|.

(b) Desenhe, no mesmo grafico, as funcoes F (x), F 2(x), F 3(x) e G(x) =x. Que pode concluir acerca dos pontos fixos e pontos periodicosde xn+1 = F (xn)?


(c) Faca uma tabela, ou desenhe um grafico, de n vs xn, entre n = 0 en = 100, com os valores iniciais 0.5, 0.6, 0.89, 0.893, 0.1111111111.Discuta e explique os resultados obtidos.

(d) Na alınea anterior, em todos os casos a sequencia permanece cons-tante a partir de n = 55. Obtenha novamente os resultados daalınea anterior usando o programa com uma maior precisao numerica:

evolution60(f, x0, n) :=

block([x:bfloat(x0),xlist:[0, x0],fpprec:60],

for i thru n

do (x:ev(f), xlist:append(xlist,[i,float(x)])),

openplot_curves([["plotpoints 1 nolines"],xlist])

)$

que pode concluir?


11. Sistemas em duas dimensoes e fractais

11.1 Equacoes de diferencas de segunda ordem

Em alguns sistemas dinamicos discretos o estado do sistema num perıododepende do estado nos ultimos dois perıodos. Por exemplo, a sucessao deFibonacci, 1, 1, 2, 3, 5, . . ., define-se a partir da relacao de recorrencia

xn+2 = xn + xn+1 (11.1)

esta equacao de segunda ordem e equivalente a um sistema de duas equacoesde primeira ordem: se definirmos uma sucessao yn igual a xn+1, obtemos osistema

xn+1 = yn (11.2)

yn+1 = xn + yn (11.3)

Com os valores iniciais x0 = 1 e y0 = 1 pode-se obter um termo da sucessao deFibonacci xn a partir dos termos anteriores das duas sequencias, xn−1 e yn−1.

Em geral um sistema de equacoes de segunda ordem define-se a partirde duas funcoes de duas variaveis F (x, y) e G(x, y)

xn+1 = F (xn, yn) (11.4)

yn+1 = G(xn, yn) (11.5)

A evolucao do sistema pode ser representada num grafico em duas di-mensoes por meio de um programa simples

(C1) evolution2d(f, g, x0, y0, n) :=

block

([x: x0, y: y0, xaux, xlist:[x0, y0], numer: true],

for i thru n do

(xaux: ev(f),


y: ev(g),

x: xaux,

xlist: append(xlist, [x, y])),

openplot_curves([["plotpoints 1 nolines"], xlist]))$

A figura 11.1 mostra a evolucao do sistema de equacoes 11.2 e 11.3.

xn

yn

2 4 6 8

2

4

6

8

Figura 11.1: Evolucao do sistema de Fibonacci

Um sistema discreto em duas dimensoes que apresenta uma solucaocaotica, e o chamado mapa de Henon

xn+1 = 1 + 0, 3 yn − 1, 4 x2n (11.6)

yn+1 = xn (11.7)

A figura 11.2 mostra o resultado obtido com 10000 iteracoes:

(C2) evolution2d(1+0.3*y-1.4*x^2, x, 0.2, 0.2, 10000);

O diagrama da figura 11.2 e uma fractal, isto e, uma figura que apre-senta partes que se repetem recursivamente em escalas menores e menores.Essa propriedade pode ser corroborada se obtivermos um grafico ampliado deuma parte da figura. O lado direito da figura 11.2 foi obtido ampliando aparte que aparece num quadro no lado esquerdo. Algumas linhas que pareci-am ser simples, resultam ser realmente um conjunto de 3 linhas paralelas. Secontinuarmos a ampliar a figura, cada linha resultara ser realmente outras 3linhas.

O programa usado para amplificar uma parte do diagrama no epacode fase, e uma versao modificada do programa evolution2d que so guarda

Sistemas em duas dimensoes e fractais 119

Figura 11.2: Mapa de Henon.

as coordenadas dos pontos que cairem dentro de uma caixa delimitada pelascoordenadas xmin, xmax, ymin e ymax. O programa completo e:

(C3) window2d(f, g, x0, y0, n, xmin, ymin, xmax, ymax) :=

block

([x: x0, y: y0, xaux, xlist:[], numer: true],

for i thru n do

(if x>xmin and y>ymin and x<xmax and y<ymax then

xlist: append(xlist, [x, y]),

xaux: ev(f),

y: ev(g),

x: xaux),

openplot_curves([["plotpoints 1 pointsize 0.7

nolines"], xlist]))$

11.2 Fractais

Um metodo que produz sistemas caoticos e atractores estranhos em duasdimensoes e o jogo do caos: escolhe-se um conjunto de pontos fixos no pla-no; a seguir desenha-se o ponto inicial (x0, y0) em qualquer parte e calcula-se


(x1, y1) deslocando o ponto inicial em direccao de um dos pontos fixos, escolhi-do aleatoriamente. O deslocamento em direccao do ponto fixo e uma fraccao β

da distancia; por exemplo, se β = 1/2, mede-se a distancia desde (x0, y0) ate oponto fixo, e desloca-se o ponto num meio da distancia. O processo repete-seindefinidamente, escolhendo em cada passo um dos pontos fixos em forma aaleatoria.

Em termos matematicos, se as coordenadas do ponto fixo forem (pi, qi),o ponto (xn+1, xn) calcula-se a partir do ponto (xn, yn) por meio de

xn+1 = pi + β(xn − pi) (11.8)

yn+1 = qi + β(yn − qi) (11.9)

Uma possıvel implementacao desse algoritmo em Maxima e:

(C4) chaosgame(point,p0,b,n) :=

block ([p:ev(p0,numer),j,m:length(point),game_dat,

numer:true, display2d:false],

with_stdout ("game.dat",

print("game_dat: [ ", p),

for i thru n do

(j: random(m) + 1,

p: ev(point[j] + b*(p-point[j]), numer),

print(", ",p)),

print("]$"),

batchload("game.dat")),


nolines"],game_dat]))$

Como primeiro exemplo, usaremos tres pontos atractores nos verticesde um triangulo equilatero e um factor de contraccao β = 1/2, com 25000iteracoes

(C5) chaosgame([[0,0],[0.5,sqrt(3)/2],[1,0]],[0.1,0.1],

1/2,25000);

Obtem-se uma figura fractal designada de triangulo de Sierpinski(lado esquerdo da figura 11.3). O lado direito da figura 11.3 mostra um atractorestranho com 5 pontos, nos vertices e no centro de um quadrado, com factorde contraccao igual a 1/2.


Figura 11.3: Triangulo de Sierpinski e quadrado fractal.

11.3 Sistemas iterativos de funcoes

As equacoes 11.8 e 11.9 podem ser reescritas em forma matricial:(xn+1

yn+1

)=(

β 00 β

)(xn − pi

yn − qi

)+(

pi

qi

)(11.10)

a matriz com β na diagonal e 0 fora da diagonal, corresponde a uma contraccaodas duas coordenada x e y, no mesmo factor β. Podemos generalizar o sistemaadmitindo diferentes factores de contraccao, βx e βy nos dois eixos. Podemosainda rodar o sistema de eixos em cada iteracao, multiplicando por uma matrizde rotacao:(

xn+1

yn+1

)=(

cos θ − sin θ

sin θ cos θ

)(βx 00 βy

)(xn − pi

yn − qi

)+(

pi

qi

)(11.11)

onde θ e o angulo de rotacao. Finalmente, se quisermos incluir tambem de-formacoes e reflexoes do sistema de coordenadas, podemos simplesmente usaruma matriz com 4 parametros a, b, c e d, e um vector com 2 parametros e e f(

xn+1

yn+1

)=(

a b

c d

)(xn

yn

)+(

e

f

)(11.12)

Podemos definir m matrizes e vectores e depois escolher um grupo emforma aleatoria e aplicar a transformacao 11.12. Se o processo for repetidomuitas vezes, o resultado final sera independente do ponto inicial escolhido eproduzira uma fractal. Para controlar melhor a densidade de pontos em di-ferentes partes da fractal, convem que algumas das m transformacoes sejamescolhidas com maior frequencia. Isso controla-se usando diferentes probabili-dades para a escolha de cada transformacao.


O programa ifs a seguir a este paragrafo implementa esse metodo, desig-nado de sistema iterativo de funcoes (em ingles Iterated Functions System,IFS). Os parametros de entrada do programa sao uma lista de probabilidadesacumuladas (cada elemento i na lista e a soma das probabilidades das trans-formacoes 1, 2, ate i), uma lista de m matrizes, uma lista de m pontos, umponto inicial, e o numero de iteracoes que deverao ser feitas.

(C6) ifs(prob, mat, point, p0, n) :=

block ([p:ev(p0,numer),s,r,m:last(prob),ifs_dat,

numer:true,display2d:false],

with_stdout ("ifs.dat",

print("ifs_dat: [ ",p),

for i thru n do

(r: random(m) + 1,

s: 0,

for j while (r-prob[j]) > 0 do s:j,

p: ev(mat[s+1].p + point[s+1], numer),

print(", ",first(transpose(p)))),

print("]$")),

batchload("ifs.dat"),


nolines"], ifs_dat]))$

Um exemplo proposto por Barnsley consiste em usar as quatro matrizese vectores seguintes

A1 =(

0.85 0.04−0.04 0.85

)b1 =

(0

1.6

)(11.13)

A2 =(

0.2 −0.260.23 0.22

)b2 =

(0

1.6

)(11.14)

A3 =(−0.15 0.280.26 0.24

)b3 =

(0

0.44

)(11.15)

A4 =(

0 00 0.16

)b4 =

(00

)(11.16)

com probabilidades

p1 = 0.85 p2 = 0.07 p3 = 0.07 p4 = 0.01 (11.17)


Os comandos do Maxima para definir as matrizes e executar o programaifs sao:

(C7) a1: matrix([0.85,0.04],[-0.04,0.85])$

(C8) a2: matrix([0.2,-0.26],[0.23,0.22])$

(C9) a3: matrix([-0.15,0.28],[0.26,0.24])$

(C10) a4: matrix([0,0],[0,0.16])$

(C11) p1: [0,1.6]$

(C12) p2: [0,1.6]$

(C13) p3: [0,0.44]$

(C14) p4: [0,0]$

(C15) prob: [85,92,99,100]$

(C16) ifs(prob,[a1,a2,a3,a4],[p1,p1,p3,p4],[5,0],25000);

O resultado e a figura 11.4. O parecido com uma folha de um feto esurpreendente, tendo em conta que toda a informacao do desenho esta contidaem apenas 4 matrizes e quatro vectores. De facto, Barnsley definiu as matrizesacima olhando para uma fotografia de uma folha real, e o resultado reproduzmuito fielmente a fotografia usada. Cada um dos quatro conjuntos de matriz evector transforma o desenho nos quatro desenhos menores que se mostram nolado direito da figura 11.4. Se cada um dos desenhos menores for transformadonovamente em forma recursiva, obtem-se a folha completa.

Figura 11.4: Fractal em forma de folha, gerada com 4 transformacoes.

Para construir fractais, identificam-se as partes necessarias para definira figura completa em forma recursiva. Os 6 parametros de cada transformacaocalculam-se a partir das 6 equacoes de transformacao para 3 pontos no desenho.


11.4 Problemas

1. Desenhe as duas fractais obtidas com quatro pontos atractivos, nosvertices de um quadrado, com factores de contracao iguais a 0, 45 e0, 55.

2. Os tres operadores A1, A2 e A3 que atraim um ponto qualquer para osvertices do triangulo formado por (0, 0), (1, 0) e (0, 1) sao

A1

(x

y

)=

12

(x

y

)A2

(x

y

)=

12

(x− 1

y

)+(

10

)A3

(x

y

)=

12

(x

y − 1

)+(

01

)onde o factor de contraccao e igual a 1/2.

(a) Desenhe os diagramas de evolucao dos 3 sistemas dinamicos(xn+1

yn+1

)= Ai

(xn

yn

)para cada um dos 3 operadores A1, A2 e A3, com quaisquer valoresiniciais (x0, y0).

(b) Desenhe o diagrama de evolucao do sistema(xn+1

yn+1

)= A3

(A2

(A1

(xn

yn

)))com valores iniciais (0.5, 0.5), e explique como e a solucao quandon e se aproxima de infinito.

(c) Desenhe a fractal obtida quando os tres operadores A1, A2 e A3

actuam alternadamente em forma aleatoria.

3. Considere o atractor estranho do problema anterior com os 3 vertices dotriangulo formado por (0, 0), (1, 0) e (0, 1), e com factor de contraccaoigual a 1/2. Para alem do factor de contraccao 1/2, incluia em cadapasso uma rotacao de 45, definida pela matriz de rotacao:(

cos(π/4) − sin(π/4)sin(π/4) cos(π/4)

)Desenhe a fractal obtida com esse sistema.


Indice remissivo

Aamortecedor, 73amortecimento crıtico, 74amplitude, 71analıtica, 30

Bbifurcacao, 88

Ccampo de direccoes, 37centro atractivo, 112centro repulsivo, 112ciclo, 114ciclos de fronteira, 98condicoes fronteira, 52condicoes iniciais, 52constante de fase, 71coordenadas normais, 94

Ddiagrama de degraus, 109

Eequacao autonoma, 39equacao de diferencas, 104equacao exacta, 45equacao inversa, 36equilıbrio estavel, 91equilıbrio instavel, 91espaco de fase, 2, 53estado, 1

Ffactor integrante, 45Fibonacci, 19fractal, 119frequencia angular, 72

JJacobiana, 93jogo do caos, 120

Llinear homogenea, 51

Mmetodo de Newton, 105metodo de iteracao, 107mapa de Henon, 119modelo logıstico, 41modo de oscilacao, 82

Nno atractivo, 112no repulsivo, 112

Oorbita, 104oscilacoes auto-excitadas, 98

Pponto de sela, 94ponto fixo, 93, 112pontos fixos, 88


pontos periodicos, 114

Rregime estacionario, 75regra da cadeia, 23relacao aurea, 19relacao de recorrencia, 105ressonancia, 75

Sserie binomial, 31serie de Taylor, 31sistema iterativo de funcoes, 123sobreamortecimento, 74solucao geral, 39solucao particular, 39

Ttriangulo de Sierpinski, 121

Vvalores proprios, 80variaveis de estado, 1, 53vectores proprios, 80

Indice remissivo 127

Documents

Introdu¸cão aos Sistemas Dinâmicos usando Maximamargarete/download/MET200-0/sistemasdinamicos.pdf · Introdu¸cão aos Sistemas Dinâmicos usando Maxima Jaime E. Villate Faculdade