Prévia do material em texto

Métodos Matemáticos
para Engenharia
Química
uma abordagem prática com o software livre Maxima
Adilson José de Assis
Métodos Matemáticos para Engenharia
Química: uma abordagem prática com o
software livre Maxima
Adilson José de Assis
Faculdade de Engenharia Química
Universidade Federal de Uberlândia
http://www.feq.ufu.br
Métodos Matemáticos para Engenharia Química: uma abordagem prática com o software livre
Maxima
Copyright © 2015 Adilson J. de Assis
E-mail: ajassis@ufu.br
Versão: 17-02-2015
ISBN:
Este livro pode ser copiado e reproduzido livremente, respeitando os termos da Licença Creative
Commons Atribuição-Partilha (versão 3.0). Para obter uma cópia desta licença, visite
http://creativecommons.org/licenses/by-sa/3.0/ ou envie uma carta para Creative
Commons, 559 Nathan Abbott Way, Stanford, California 94305, USA.
�A Matemática é o alfabeto com o qual Deus escreveu o Universo.�
�O livro da Natureza está escrito em caracteres matemáticos ...sem um conhecimento dos mesmos os homens
não poderão compreendê-lo�.
Galileu Galilei
�A matemática, corretamente vista, possui não apenas a verdade, mas suprema beleza - uma beleza fria e
austera, sem os enfeites lindos da pintura ou da música.�
Bertrand Russell
'
&
$
%
A natureza e a matemática são complementares e belas. Veja os exemplos nos vídeos a seguir:
https://vimeo.com/77330591: Beauty of mathematics
https://www.youtube.com/watch?v=CzqQyamt6Jw: A beleza da matemática
https://www.youtube.com/watch?v=QaWepnGWRs8: Sequência de Fibonacci e Número de Ouro
https://www.youtube.com/watch?v=BTiZD7p_oTc: Best fractals zoom ever
3
Sumário
Sumário 4
Lista de Figuras 7
1. Introdução: modelagem matemática, engenharia e equações diferenciais 9
2. Revisão de Matemática Básica e Cálculo Diferencial e Integral 12
3. Breve apresentação do software livre Maxima 19
3.1. Iniciando o uso do Maxima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.2. Processando um comando . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.3. Manipulação de variáveis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4. Funções matemáticas elementares no Maxima . . . . . . . . . . . . . . . . . . . . . . . 22
3.5. Operações básicas com matrizes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.6. Solução de sistemas lineares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.7. Solução analítica de uma equação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.8. Solução numérica de uma equação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.9. Manipulação algébrica de equações e funções . . . . . . . . . . . . . . . . . . . . . . . 26
3.10. Simplificações diversas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.11. Gráficos bi e tridimensionais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.12. Limites e continuidade de funções . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.13. Derivadas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.13.1. Derivadas simples e pontos críticos de uma função . . . . . . . . . . . . . . . . 31
3.13.2. Derivadas parciais de uma função . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.14. Integrais indefinidas, definidas e área abaixo de uma curva . . . . . . . . . . . . . . . . 32
3.15. Programação no Maxima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.15.1. Instrução: if . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.15.2. Instrução: for...end . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.15.3. Instrução block . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.16. Séries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.16.1. Série de Taylor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.16.2. Série de Fourier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.17. Regressão Linear . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4
SUMÁRIO
4. Campo de direções no Maxima 39
4.1. Exemplo 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.2. Exemplo 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.3. Exemplo 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.4. Exemplo 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.5. Função drawdf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5. Solução analítica de equações diferenciais ordinárias (EDOs) 45
5.1. Verificando se uma dada equação é solução de uma EDO . . . . . . . . . . . . . . . . . 45
5.2. Equações diferenciais de 1ª ordem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
5.2.1. Encontrando a solução geral de uma ED de 1ª ordem . . . . . . . . . . . . . . . 45
5.2.2. Encontrando a solução particular de uma ED de 1ª ordem . . . . . . . . . . . . 46
5.2.3. EDO com solução implícita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.2.4. Solução de EDOs exatas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
5.2.5. Exemplo 1: variação da massa de um soluto com o tempo a volume constante . 48
5.2.6. Exemplo 2: lei do resfriamento de Newton . . . . . . . . . . . . . . . . . . . . . 49
5.3. Equações diferenciais lineares de 2ª ordem . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.3.1. Solução geral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
5.3.2. Solução de PVIs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.3.3. Exemplo: difusão e reação de um componente gasoso em um líquido (Rice and
Do., 1995, exemplo 2.7) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.3.4. Solução de PVCs . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
5.3.5. EDOs não homogêneas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
5.3.6. Exemplo: Dinâmica de um tanque agitado com aquecimento elétrico . . . . . . 55
5.3.7. Wronkiano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.4. Pacote contrib_ode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
5.5. Solução usando série de potências para EDOs lineares de 2ª ordem . . . . . . . . . . . 59
5.6. Solução de EDOs por transformada de Laplace (EDOs simples e sistemas de EDOs) . 62
5.6.1. Transformada de Laplace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
5.6.2. Solução de EDOs por transformada de Laplace . . . . . . . . . . . . . . . . . . 63
5.7. Sistemas de EDOs lineares de 1ª ordem no plano . . . . . . . . . . . . . . . . . . . . . 65
5.7.1. Autovalores e autovetores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.7.2. Retrato de fase de sistemas lineares . . . . . . . . . . . . . . . . . . . . . . . . . 67
6. Solução numérica de equações diferenciais ordinárias usando o método de Runge-Kutta 72
7. Solução analítica de equações diferencias parciais (EDPs) 75
7.1. Exemplos de EDPs na Engenharia Química . . . . . . . . . . . . . . . . . . . . . . . . 75
7.2. Séries de Fourier . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
7.3. Solução de EDPs por separação de variáveis . . . . . . . . . . . . . . . . . . . . . . . . 82
5
SUMÁRIO
8. Solução numérica deequações diferencias parciais (EDPs) usando o método das diferenças
finitas 88
8.1. EDP parabólica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
8.2. EDP hiperbólica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
8.3. EDP elíptica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
Referências Bibliográficas 94
A. Lista de exercícios n. 1: Revisão de matemática básica e cálculo diferencial e integral 95
B. Lista de exercícios n. 2: Classificação, campo de direções e EDOs de 1ª ordem 98
C. Lista de exercícios n. 3: EDOs de 2ª ordem, solução em série de potências de EDOs 103
D. Lista de exercícios n. 4: solução de EDOs por transformada de Laplace, sistemas de EDOs107
E. Lista de exercícios n. 5: EDPs 112
F. Revisão de matemática básica e de cálculo diferencial e integral: uma abordagem prática
com o Maxima 115
6
Lista de Figuras
1.0.1.Matemática como literatura de suspense, segundo Calvin! . . . . . . . . . . . . . . . . 10
2.0.1.Definição das funções trigonométricas . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.0.2.Gráficos das funções exponencial e logarítmica . . . . . . . . . . . . . . . . . . . . . . . 13
2.0.3.Gráficos das funções hiperbólicas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3.1.1.Imagem de uma tela do Maxima, com a interface gráfica wxMaxima. . . . . . . . . . . 20
3.8.1.Função sin(x) exp(x) + x− 3 usada para o cálculo numérico de suas raízes. . . . . . . . 25
3.11.1.Gráficos bidimensionais no Maxima. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.11.2.Gráfico bidimensional com pontos discretos e curva contínua. . . . . . . . . . . . . . . 30
3.11.3.Gráficos tridimensionais no Maxima. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.12.1.Limites simples e laterais no Maxima. . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.13.1.Função usada no exemplo de derivadas e pontos críticos. . . . . . . . . . . . . . . . . . 32
3.14.1.Função usada no exemplo de integração. . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.16.1.Aproximações por série de Taylor para a função cos(x). . . . . . . . . . . . . . . . . . 36
3.16.2.Série de Fourier da função usada no exemplo. . . . . . . . . . . . . . . . . . . . . . . . 37
3.17.1.Exercício de aplicação de regressão linear. . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.0.1.Campo de direções como vetores tangentes à solução de uma EDO. . . . . . . . . . . . 40
4.1.1.Campo de direções: exemplo 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
4.2.1.Campo de direções: exemplo 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.3.1.Campo de direções: exemplo 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.4.1.Campo de direções: exemplo 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.5.1.Campo de direções com a função drawdf. . . . . . . . . . . . . . . . . . . . . . . . . . . 44
5.2.1.Solução de um PVI, com equação diferencial de 1ª ordem. . . . . . . . . . . . . . . . . 47
5.2.2.Variação da massa com o tempo para um tanque agitado (m0=1, c0=1, V= 2, r=1) . . 49
5.2.3.Solução do problema envolvendo a lei de resfriamento de Newton. . . . . . . . . . . . . 50
5.3.1.Funções de Bessel J0(x) e Y0(x). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
5.3.2.PVI envolvendo uma EDO de 2ª ordem. . . . . . . . . . . . . . . . . . . . . . . . . . . 52
5.3.3.Difusão e reação de um componente gasoso. . . . . . . . . . . . . . . . . . . . . . . . . 54
5.3.4.Solução de uma EDO de segunda ordem, não homogênea. . . . . . . . . . . . . . . . . 55
5.3.5.Tanque agitado com aquecimento elétrico. . . . . . . . . . . . . . . . . . . . . . . . . . 57
7
LISTA DE FIGURAS
5.6.1.Solução de um sistema de EDOs, por transformada de Laplace, sem usar o comando
desolve. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.7.1.Retrato de fase e evolução no tempo. O gráfico da direita é obtido clicando no último
ícone da barra superior da figura à esquerda. . . . . . . . . . . . . . . . . . . . . . . . 68
5.7.2.Retrato de fase: nó atrator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.7.3.Retrato de fase: foco atrator. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6.0.1.Solução numérica de EDOs: única EDO. . . . . . . . . . . . . . . . . . . . . . . . . . . 73
6.0.2.Solução numérica de EDOs: sistema de duas EDOs. . . . . . . . . . . . . . . . . . . . 73
7.1.1.Reator modelado no Exemplo 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
7.1.2.Solução, na forma gráfica, da EDP do Exemplo 1 (Rice and Do., 1995). . . . . . . . . . 77
7.1.3.Solução gráfica do Exemplo 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
7.1.4.Solução analítica na forma gráfica do Exemplo 4 (Rice and Do., 1995). . . . . . . . . . 79
7.1.5.Solução analítica na forma gráfica do Exemplo 5 (Rice and Do., 1995). . . . . . . . . . 81
7.2.1.Exemplo de série de Fourier no Maxima. . . . . . . . . . . . . . . . . . . . . . . . . . . 83
7.3.1.Resposta, na forma tridimensional, para o exemplo numérico da condução de calor em
uma barra, com extremidades a temperaturas fixas. . . . . . . . . . . . . . . . . . . . . 86
7.3.2.Variação da temperatura, com o tempo, no centro da barra. . . . . . . . . . . . . . . . 86
7.3.3.Variação da temperatura, com a posição, para vários tempos diferentes. . . . . . . . . 87
8.1.1.Solução numérica da equação da difusão, transiente, unidimensional. . . . . . . . . . . 89
8.2.1.Solução numérica da equação da onda, transiente, unidimensional. . . . . . . . . . . . 90
8.3.1.Solução numérica da equação de Laplace. . . . . . . . . . . . . . . . . . . . . . . . . . 92
8.3.2.Solução numérica da equação de Poisson. . . . . . . . . . . . . . . . . . . . . . . . . . . 93
B.0.1.Campo de direções de uma determinada equação diferencial. . . . . . . . . . . . . . . . 100
B.0.2.Tanque cilíndrico com entrada e saída de líquido. . . . . . . . . . . . . . . . . . . . . . 101
D.0.1.Dois tanques aquecidos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
D.0.2.Campo das direções de um circuito elétrico. . . . . . . . . . . . . . . . . . . . . . . . . 109
D.0.3.Sistema massa-mola, com dois corpos e três molas. . . . . . . . . . . . . . . . . . . . . 110
E.0.1.Esquema 2D de um aquífero. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
8
1. Introdução: modelagem matemática,
engenharia e equações diferenciais
�A modelagem matemática é a area do conhecimento que estuda a simulação de sistemas reais a fim
de prever o comportamento dos mesmos, sendo empregada em diversos campos de estudo, tais como
física, química, biologia, economia e engenharias. Ou seja, modelagem matemática consiste na arte
(ou tentativa) de se descrever matematicamente um fenômeno.
Dentre as diferentes formas e métodos de modelagem temos a modelagem via autômatos celulares e
equações diferenciais, parciais e/ou ordinárias. A modelagem matemática via equações diferenciais
tem um papel de enorme destaque, visto que tal técnica vem sendo utilizada para modelar fenômenos
desde o século XVII por Malthus e Verhulst,no final dos anos 1700. Pode-se, então, dizer que um
modelo matemático é desenvolvido para simular a realidade usando a linguagem matemática.
Os modelos matemáticos se subsidiam, por exemplo, das leis da física (como as leis de Kirchhoff para
sistemas eletricos e as leis de Newton para mecânicos) ou dados experimentais.
Frequentemente, os modelos atingem grau de sofisticação suficiente para justificar ferramentas compu-
tacionais, envolvendo sistemas deequações diferenciais.
A modelagem de um fenômeno via equações diferenciais é, normalmente, feita da seguinte forma:
através da simples observação conseguem-se informações sobre as taxas de variação do fenômeno (que
do ponto de vista matemático são derivadas), escreve-se a equação que relaciona as taxas de variação
e a função, isto é, a equação diferencial associada e, a partir da solução desta equação tem-se uma
possível descrição do fenômeno.
Então, tal modelo matemático será também composto por parâmetros (constantes), que são intrínsecas
ao sistema a ser estudado; variáveis que afetam o sistema, porém o modelo não foi designado para
estudar seu comportamento (variáveis independentes) e as variáveis as quais o modelo foi designado
para estudar (variáveis dependentes). Quando o sistema em questão busca retratar um fenômeno que
consiste na interação entre duas ou mais entidades, então a modelagem é feita através de um sistema
de equações diferenciais
9
1. Introdução: modelagem matemática, engenharia e equações diferenciais
Os modelos matemáticos apresentam uma série de aspectos úteis do ponto de vista científico. Além de
apresentar naturalmente uma linguagem concisa, que pode vir a facilitar sua manipulação, um modelo
matemático traz também aspectos como a possibilidade de confirmar ou rejeitar determinadas hipóteses
relacionadas a complexos sistemas, revelar contradições em dados obtidos e/ou hipóteses formuladas,
prever o comportamento de um sistema sob condições não testadas ou ainda não �testáveis�, dentre
outros.
Por outro lado, quanto maior é a proximidade do modelo com a realidade, mais complexo será o
modelo. Isto significa um maior numero de parâmetros e conseqüentemente uma maior dificuldade
tanto na obtenção de dados a partir do modelo quanto na interpretação desses dados gerados pelo
modelo em questão.
Figura 1.0.1.: Matemática como literatura de sus-
pense, segundo Calvin!
Modelos simples são mais fáceis de lidar, porém
modelos mais sofisticados são frequentemente ne-
cessários. É importante ressaltar que as previ-
sões do comportamento de um determinado mo-
delo matemático, caso se faça necessário depen-
dendo de sua complexidade, se dão através de si-
mulações computacionais do mesmo. Caso o mo-
delo seja suficientemente simples, teorias mate-
máticas são eficientes ferramentas para se obter
conclusões gerais. Então, pode-se dizer que ao
desenvolver um modelo matemático busca-se um
ponto ótimo entre a representação da realidade e
a complexidade do modelo, para que a obtenção
de resultados coerentes seja possível, bem como
sua interpretação. Segundo Howard Emmons, �o
desafio em modelagem matemática não é produ-
zir os modelos descritivos mais compreensíveis,
mas sim produzir modelos suficientemente sim-
ples que incorporam as principais características
do fenômeno em questão�. Portanto, a modela-
gem matemática ajuda a evitar ou reduzir a ne-
cessidade de gastos excessivos em experimentos,
ou até mesmo simular experimentos impossíveis
de serem realizados na prática.�Wikipédia (2013)
Particularmente, em Engenharia Química, os mo-
delos matemáticos na forma de equações diferencias tem um papel extremamente importante, pois
servem para descrever matematicamente:
1. Fenômenos físicos e químicos, tais como a transferência de quantidade de movimento (mecânica
dos fluidos), a transferência de calor e de massa, a velocidade de reações (bio)químicas;
2. Os balanços de massa, energia e quantidade de movimento em equipamentos (ou operações
unitárias);
3. A dinâmica de processos no estado transiente, base para o controle de processos;
Espero que este material contribua para que o universo das equações diferenciais deixe de ser um
�mistério�, como é para o Calvin, e se torne algo palpável em sua vida! LOL
No conteúdo que se encontra a seguir, o leitor encontrará uma abordagem prática, com o uso de um
�software� livre e gratuito, dos métodos matemáticos para a Engenharia Química, com ênfase na
solução analítica de equações diferenciais. A teoria dos métodos utilizados não será coberta aqui
e o leitor deve consultar os manuais clássicos de equações diferenciais, como por exemplo Kreyszig
et al. (2011) e Boyce and DiPrima (2010), ou materiais disponíveis gratuitamente na rede, como por
exemplo Villate (2011).
10
1. Introdução: modelagem matemática, engenharia e equações diferenciais
Este livro parte da premissa que conhecer o processo matemático de, por exemplo, resolver uma equação
diferencial ordinária é importante. Entretanto, com o advendo e o avanço da informática, permitindo
hoje o acesso a �softwares� que podem fazer boa parte do trabalho antes as vezes trabalhoso ou tedioso,
de manipulação de equações, o engenheiro pode concentrar seus esforços mais nas etapas de desenvol-
vimento do modelo matemático e na análise da solução, por economizar tempo na etapa de solução
do modelo. Com isto, sua atuação enfatiza mais a ênfase como engenheiro e não como matemático.
Com isto não se quer dizer que o lado matemático do engenheiro não seja importante, até porque
todos os softwares disponíveis são limitados e não livres de erros; porém, no estado da arte em que eles
se encontram, pelo menos os métodos clássicos estão relativamente bem implementados, permitindo
principalmente que a formação acadêmica em métodos matemáticos considere também os aspectos
práticos, permitindo ao estudante avançar na aplicação nos métodos estudados e o desenvolvimento de
projetos didáticos mais elaborados ou mais próximos dos modelos que serão encontrados nos ambientes
de atuação profissional onde tais modelos desempenham um papel relevante.
11
2. Revisão de Matemática Básica e Cálculo
Diferencial e Integral
A seguir será apresentada uma breve revisão de matemática básica e de cálculo diferencial e integral,
com o objetivo de melhor compreender os conteúdos que serão abordados posteriormente. No Apêndice
A encontra-se esta revisão com os comandos do Maxima apropriados.
ˆ (x± y)2 = x2 ± 2xy + y2
ˆ x2 − y2 = (x− y)(x+ y)
ˆ Fatorial de n: n! = 1.2.3. · · · .n
ˆ Para o triângulo retângulo mostrado na Fig. 2.0.1:
* Teorema de Pitágoras: OB
2
= OA
2
+AB
2
* sen2α+ cos2α = 1
Figura 2.0.1.: Definição das funções trigonométricas
ˆ ap.aq = ap+q
ˆ ap/aq = ap−q
ˆ (ap)q = apq
ˆ a0 = 1, a 6= 0
ˆ 0a = 0, a 6= 0
ˆ
0
5
= 0, mas
5
0
é indefinido
12
2. Revisão de Matemática Básica e Cálculo Diferencial e Integral
ˆ a−p = 1/ap, a 6= 0
ˆ (ab)p = apbp
ˆ
(a
b
)−n
=
(
b
a
)n
=
bn
an
ˆ a
n
m =
(
a
1
m
)n
= (an)
1
m
ˆ
n
√
a = a1/n
ˆ
n
√
am = am/n
ˆ
n
√
a/b = n
√
a/ n
√
b
ˆ Cuidado!!!
* (a+ b)2 6= a2 + b2
*
√
a2 + b2 6= a+ b
*
3a+ b
3a+ c
6= b
c
*
dpi3
dx
6= 3pi2
ˆ Função logarítmica −→ y = loga x. Ver Fig. 2.0.2
Figura 2.0.2.: Gráficos das funções exponencial e logarítmica
ˆ logaMN = logaM + logaN
ˆ logaM/N = logaM − logaN
ˆ logaM
p = p logaM
ˆ blogb x = x
ˆ i =
√−1
ˆ i2 = −1
ˆ
√−a = i√a se a ≥ 0
ˆ eiθ = cos θ + i senθ
ˆ e−iθ = cos θ − i senθ
ˆ sen θ =
eiθ − e−iθ
2i
13
2. Revisão de Matemática Básica e Cálculo Diferencial e Integral
ˆ cos θ =
eiθ + e−iθ
2
ˆ Seno hiperbólico de x : senh x =
ex − e−x
2
Ver Fig. 2.0.3.
Figura 2.0.3.: Gráficos das funções hiperbólicas
ˆ Cosseno hiperbólico de x : cosh x =
ex + e−x
2
ˆ Equação quadrática ou do 2º grau: ax2 + bx + c = 0; raízes: x =
−b±√b2 − 4ac
2a
. Se a, b e c
são reais e se D = b2 − 4ac é o discriminante, as raízes da equação são:
i) reais e desiguais se D > 0
ii) reais e iguais se D = 0
iii) conjugadas e complexas se D < 0
ˆ lim
x→a bf(x) = b limx→a f(x)
ˆ lim
x→a[f(x) + g(x)] =limx→a f(x) + limx→a g(x)
ˆ lim
x→a[f(x)g(x)] = limx→a f(x) · limx→a g(x)
ˆ lim
x→a
f(x)
g(x)
=
limx→a f(x)
limx→a g(x)
, se lim
x→a g(x) 6= 0
ˆ lim
x→a[f(x)]
n = [lim
x→a f(x)]
n
ˆ lim
x→∞ e
x =∞ e lim
x→−∞ e
x = 0
ˆ lim
x→∞ ln(x) =∞ e limx→0+ ln(x) = −∞
ˆ Descontinuidades:
* Removível: y = sen x/x, quando x = 0
* Infinita: y = 1/x, quando x = 0
* Salto: y = 10/(1 + e1/x), quando x = 0+ y = 0+ ou x = 0 y = 0 ou x = 0− y = 10
14
2. Revisão de Matemática Básica e Cálculo Diferencial e Integral
ˆ Se y = P (x), a derivada de y ou de P(x) em relação à x é definida como
dy
dx
= lim
h→0
f(x+ h)− f(x)
h
=
lim
∆x→0
f(x+ ∆x)− f(x)
∆x
ˆ A derivada é também designada por y′, df/dx ou f ′(x). O processo de derivar/obter uma derivada
é chamado diferenciação.
ˆ
d
dx
(c) = 0
ˆ
d
dx
(cx) = c
ˆ
d
dx
(cxn) = ncxn−1
ˆ
d
dx
(u± v ± w ± . . .) = du
dx
± dv
dx
± dw
dx
± . . .
ˆ
d
dx
(cv) = c
dv
dx
ˆ
d
dx
(uv) = u
dv
dx
+ v
du
dx
ˆ
d
dx
(u
v
)
=
v dudx − u dvdx
v2
ˆ
d
dx
(un) = nun−1
du
dx
ˆ Regra da cadeia:
dy
dx
=
dy
du
du
dx
ˆ
du
dx
=
1
dx/du
ˆ
dy
dx
=
dy/du
dx/du
ˆ
d
dx
senu = cosu
du
dx
ˆ
d
dx
cosu = −senudu
dx
ˆ
d
dx
loga u =
loga e
u
du
dx
, a 6= 0, 1
ˆ
d
dx
lnu =
d
dx
loge u =
1
u
du
dx
ˆ
d
dx
au = au ln a
du
dx
ˆ
d
dx
eu = eu
du
dx
ˆ
d
dx
uv =
d
dx
ev lnu = ev lnu
d
dx
[v lnu] = vuv−1
du
dx
+ uv lnu
dv
dx
ˆ
d
dx
senh u = coshu
du
dx
ˆ
d
dx
cosh u = senh u
du
dx
ˆ Derivadas superiores −→ segunda derivada: d
dx
(
dy
dx
)
=
d2y
dx2
= f”(x) = y”
ˆ Diferencial: d(u± v ± w ± . . .) = du± dv ± dw ± . . .
15
2. Revisão de Matemática Básica e Cálculo Diferencial e Integral
ˆ Diferencial: d(uv) = udv + vdu
ˆ Seja f(x, y), uma função de duas variáveis x e y. Então definimos a derivada parcial de f(x, y)
com relação a x, considerando y constante, por
∂f
∂x
= lim
∆x→0
f(x+ ∆x, y)− f(x, y)
∆x
Igualmente a derivada parcial de f(x, y) com relação a y, considerando x constante é definida
por:
∂f
∂y
= lim
∆y→0
f(x, y + ∆y)− f(x, y)
∆y
ˆ A diferencial de f(x, y) (ou derivada total de f) é definida como: df =
(
∂f
∂x
)
y
dx+
(
∂f
∂y
)
x
dy
ˆ Teorema de L'Hospital: Formas como 0/0,∞/∞, 0×∞ etc são chamadas de indeterminadas;
se lim
x→a
f(x)
g(x)
=
0
0
ou
±∞
±∞ então limx→a
f(x)
g(x)
= lim
x→a
f ′(x)
g′(x)
. Exemplo: lim
x→0
sen(x)
x
= lim
x→0
d sen(x)
dx
=
lim
x→0
cos(x)
1
= 1
ˆ Se
dy
dx
= f(x), então y é a função cuja derivada é f(x) e é chamada de anti-derivada de f(x) ou
integal indefinida de f(x), designada por y =
∫
f(x)dx
ˆ
∫
a dx = ax+ C [as constantes de integração serão omitidas, mas estão implícitas]
ˆ
∫
af(x) dx = a
∫
f(x) dx
ˆ
∫
(u± v ± w ± . . .) dx =
∫
u dx±
∫
v dx±
∫
w dx± . . .
ˆ Integração por partes:∫
u dv = uv −
∫
v du
ˆ
∫
un du =
un+1
n+ 1
, com n 6= −1
ˆ
∫
du
u
= lnu se u > 0 ou ln(−u) se u < 0 = ln |u|
ˆ
∫
eu du = eu
ˆ
∫
sen u du = − cosu
ˆ
∫
cosu du = sen u
ˆ Transformações importantes:
*
∫
F (ax+ b) dx =
1
a
∫
F (u) du onde u = ax+ b
*
∫
F (
√
ax+ b) dx =
2
a
∫
uF (u) du onde u =
√
ax+ b
*
∫
F (
n
√
ax+ b) dx =
n
a
∫
un−1F (u) du onde u = n
√
ax+ b
*
∫
F (eax) dx =
1
a
∫
F (u)
u
du onde u = eax
16
2. Revisão de Matemática Básica e Cálculo Diferencial e Integral
*
∫
F (lnx) dx =
∫
F (u)eu du onde u = lnx
ˆ Métodos de integração
1
* Fórmula direta: quando o integrando pode ser transformado numa forma facilmente
integrável;
* Substituição trigonomética: aplica-se em geral quando o integrando está na forma de
um radical;
* Substituição algébrica: aplica-se em geral quando há funções da forma (a+ bx)1/n;
* Frações parciais: redução de funções racionais na forma f(x)/g(x) para frações mais
facilmente integráveis;
* Integração por partes
* Expansão em séries
ˆ
∫ b
a
f(x) dx = g(x)|ba = g(b)− g(a)
ˆ
∫ a
a
f(x) dx = 0
ˆ
∫ b
a
f(x) dx = −
∫ a
b
f(x) dx
ˆ
∫ b
a
f(x) dx =
∫ c
a
f(x) dx+
∫ b
c
f(x) dx com a < c < b
ˆ Função Gama: Γ(n) =
∫ ∞
0
tn−1e−t dt com n > 0
ˆ Função Erro: erf(x) =
2√
pi
∫ x
0
e−u
2
du
ˆ erf(x) =
2√
pi
(
x− x
3
3.1!
+
x5
5.2!
− x
7
7.3!
+ . . .
)
ˆ erf(−x) = −erf(x), erf(0) = 0, erf(∞) = 1
ˆ Série de Taylor [expansão de uma função f(x) em torno de um ponto a]: f(x) = f(a) +
f ′(a)(x− a) + f
′′(a)(x− a)2
2!
+ . . .+
f (n−1)(a)(x− a)n−1
(n− 1)! +Rn onde Rn é o resto após n termos.
Se a = 0 a série é frequentemente chamada série de Maclaurin
ˆ ex = 1 + x+
x2
2!
+
x3
3!
+ . . . com −∞ < x <∞
ˆ ln(1 + x) = x− x
2
2
+
x3
3
− x
4
4
+ . . . com −1 < x < 1
ˆ lnx =
(
x− 1
x
)
+
1
2
(
x− 1
x
)2
+
1
3
(
x− 1
x
)3
+ . . . com x ≥ 12
ˆ sen x = x− x
3
3!
+
x5
5!
− x
7
7!
+ . . . com −∞ < x <∞
ˆ cosx = 1− x
2
2!
+
x4
4!
− x
6
6!
+ . . . com −∞ < x <∞
ˆ Laplaciano: U = ∇2U = ∇ · (∇U) = ∂
2U
∂x2
+
∂2U
∂y2
+
∂2U
∂z2
1
Consultar a literatura de Cálculo Diferencial e Integral para detalhes dos métodos.
17
2. Revisão de Matemática Básica e Cálculo Diferencial e Integral
ˆ Equação diferencial de Bessel: x2y′′ + xy′ + (x2 − n2)y = 0 com n > 0. As soluções desta
equação são chamadas funções de Bessel de ordem n.
ˆ Funções de Bessel de primeira classe de ordem n:
* Jn(x) =
∞∑
k=0
(−1)k(x/2)n+2k
k!Γ(n+ k + 1)
* J−n(x) =
∞∑
k=0
(−1)k(x/2)2k−n
k!Γ(k + 1− n)
* J−n(x) = (−1)nJn(x) n = 0, 1, 2, . . .
* J0(x) = 1− x
2
22
+
x4
22 · 42 −
x6
22 · 42 · 62 + . . .
* J1(x) =
x
2
− x
3
22 · 4 +
x5
22 · 42 · 6 −
x7
22 · 42 · 62 · 8 + . . .
* J ′0(x) = −J1(x)
ˆ Funções de Bessel de segunda classe de ordem n: Yn(x)
18
3. Breve apresentação do software livre
Maxima
3.1. Iniciando o uso do Maxima
O Maxima é um sistema de computação algébrica baseado em uma versão de 1982 do Macsyma. Ele
é escrito em Common Lisp e funciona em todas as plataformas POSIX, tais como Mac OS X, Unix,
BSD, e GNU/Linux bem como no Microsoft Windows e mesmo em dispositivos que rodam o Android.
Trata-se de um software livre
1
cuja licença é a GNU General Public License.
O Maxima é um sistema para a manipulação de expressões simbólicas, semelhantemente ao Maple e ao
Mathematica, e numéricas, semelhante ao Matlab e ao Scilab, incluindo diferenciação, integração, série
de Taylor, transformadas de Laplace, equações diferenciais ordinárias, sistemas de equações lineares,
polinômios, conjuntos, listas, vetores, matrizes e tensores. O Maxima produz resultados numéricos de
alta precisão usando frações exatas, inteiros de precisão arbitrária e números de ponto flutuante de
precisão variável. O Maxima pode plotar funções e dados em duas e três dimensões.
O Maxima é um sistema de propósito geral, e cálculos de casos especiais tais como a fatoração de
números grandes, a manipulação de polinômios extremamente grandes, etc, algumas vezes são melhor
desempenhados com sistemas especializados.
O Maxima é constituído de um núcleo, que recebe e interpreta os comandos e devolve os resulta-
dos. Entretanto, não é prático usar diretamente o núcleo, pois ele carece de opções de salvamento,
menus e configurações. Há diversas interfaces gráficas para o núcleo,sendo a mais amigável delas a
wxMaxima, mostrada na Figura 3.1.1. Link para baixar o Maxima, já com a inteface wxMaxima:
http://maxima.sourceforge.net/
Quando se instala o Maxima no computador, há na aba �Ajuda� o manual do software. Entretanto, é
a versão em Inglês que está disponível. No link a seguir, há uma tradução para o Português do manual
do Maxima: http://maxima.sourceforge.net/docs/manual/pt/maxima.html
Há várias boas apostilas na internet abordando o uso do Maxima para fins diversos. Os principais
materiais que foram usados para este tutorial foram: Santos (2009), do Prado et al. (2008), de Oliveira
(2014); Torres et al. (2006), Woollett (2009) e Nodarse (2014).
3.2. Processando um comando
Após abrir o wxMaxima, comece digitando os comandos a seguir; para o Maxima processar as entradas,
pressione simultaneamente SHIFT+ENTER. Caso seja pressionado apenas ENTER, o Maxima apenas
cria uma nova linha de comando. Cada linha em uma célula deve ser finalizada por ponto e vírgula
(;):
1
�Os objetivos do Software Livre e (controle na própria computação e cooperação livre) são atingidas por concessão do
seguinte-direitos de liberdade: os usuários são livres para executar, copiar, distribuir, estudar, mudar e melhorar o
software, estas liberdades são explicitamente concedidos e não suprimidas (como é o caso do software proprietário).
Assim, o software livre é uma questão de liberdade, não de preço (os usuários são livres - o que inclui a liberdade
de redistribuir o software, que pode ser feito gratuitamente ou por uma taxa). Software livre garante as liberdades
dos usuários: estudar e modificar software, pela disponibilidade do código fonte, bem como a liberdade de copiar e
distribuir.� Fonte: http://pt.wikipedia.org/wiki/Software_livre
19
3. Breve apresentação do software livre Maxima
Figura 3.1.1.: Imagem de uma tela do Maxima, com a interface gráfica wxMaxima.
20
3. Breve apresentação do software livre Maxima
(%i1)
5+2;
%pi;
float(%pi);
%pi, numer;
%e, numer;
(%o1) 7
(%o2) pi
(%o3) 3.141592653589793
(%o4) 3.141592653589793
(%o5) 2.718281828459045
3.3. Manipulação de variáveis
Uma ferramenta importante no Maxima é a capacidade de atribuir e manipular variáveis. Uma variável,
em programação, é um identificador ao qual se pode atribuir valores. No Maxima a instrução de
atribuição é feita com dois pontos (:), como pode ser visto no exemplo a seguir:
(%i6)
a:7;
b:3;
c:a+b /* isto é um comentário */;
eq1:a*x^2+b*x+c=0;
kill(all) /* limpa as variáveis existentes */;
eq1:a*x^2+b*x+c=0;
eq2:A*x;
eq3:eq1*eq2;
eq4:expand(%);
(%o6) 7
(%o7) 3
(%o8) 10
(%o9) 7x2 + 3x+ 10 = 0
(%o0) done
(%o1) a x2 + b x+ c = 0
(%o2) xA
(%o3) x
(
a x2 + b x+ c
)
A = 0
(%o4) a x3A+ b x2A+ c xA = 0
21
3. Breve apresentação do software livre Maxima
3.4. Funções matemáticas elementares no Maxima
(%i5) abs(-1);
4!;
sqrt(x);
x^(a/b);
log(x) /* logaritmo neperiano */;
log(%e);
log(5)/log(10), numer /* log de 5 na base 10 */;
exp(x);
sin(%pi/2);
cos(%pi/2);
tan(%pi/4);
sinh(x);
asin(1/2);
(%o5) 1
(%o6) 24
(%o7)
√
x
(%o8) x
a
b
(%o9) log (x)
(%o10) 1
(%o11) 0.69897000433601
(%o12) ex
(%o13) 1
(%o14) 0
(%o15) 1
(%o16) sinh (x)
(%o17)
pi
6
22
3. Breve apresentação do software livre Maxima
3.5. Operações básicas com matrizes
(%i18) A:matrix([5, -3, -4],
[1, -7, 12],
[10, 4, 6]);
B:matrix([1, -3, -6],
[2, -7, 12],
[10, 1, -6]);
A+B;
invert(A);
A^^(-1);
3*A /* multiplicacao por escalar */;
A.B /* multiplicacao matricial */;
invert(A).A;
determinant(A);
transpose(A);
(%o18)
 5 −3 −41 −7 12
10 4 6

(%o19)
 1 −3 −62 −7 12
10 1 −6

(%o20)
 6 −6 −103 −14 24
20 5 0

(%o21)
 45544 − 1544 117− 57544 − 35544 117
− 37544 25544 134

(%o22)
 45544 − 1544 117− 57544 − 35544 117
− 37544 25544 134

(%o23)
15 −9 −123 −21 36
30 12 18

(%o24)
−41 2 −42107 58 −162
78 −52 −48

(%o25)
1 0 00 1 0
0 0 1

(%o26) − 1088
(%o27)
 5 1 10−3 −7 4
−4 12 6

23
3. Breve apresentação do software livre Maxima
(%i1) A: matrix(
[a,b],
[c,d]
)$
determinant(A);
invert(A).A$ ratsimp(%);
A^(-1);
A.A^(-1)$ ratsimp(%);
A^^(-1).A$ ratsimp(%);
adjoint(A);
(%o2) a d− b c
(%o4)
(
1 0
0 1
)
(%o5)
(
1
a
1
b
1
c
1
d
)
(%o7)
(
c+b
c
a d+b2
b d
a d+c2
a c
c+b
b
)
(%o9)
(
1 0
0 1
)
(%o10)
(
d −b
−c a
)
3.6. Solução de sistemas lineares
(%i28) linsolve([2*x+y+z=3, x+2*y+3*z=6, 2*x+3*y+4*z=1],[x,y,z]);
2*x+y+z, x:8,y:-37,z:24 /* verificando a solucao na 1a equacao */;
(%o28) [x = 8, y = −37, z = 24]
(%o29) 3
3.7. Solução analítica de uma equação
(%i30)
eq1;
eq5:solve(eq1,x);
eq5[1] /* pegando a 1a raiz */;
eq5[2] /* pegando a 2a raiz */;
rhs(eq5[1]) /* pegando somente o resultado da 1a raiz */;
(%o30) a x2 + b x+ c = 0
(%o31) [x = −
√
b2 − 4 a c+ b
2 a
, x =
√
b2 − 4 a c− b
2 a
]
(%o32) x = −
√
b2 − 4 a c+ b
2 a
(%o33) x =
√
b2 − 4 a c− b
2 a
(%o34) −
√
b2 − 4 a c+ b
2 a
24
3. Breve apresentação do software livre Maxima
Figura 3.8.1.: Função sin(x) exp(x) + x− 3 usada para o cálculo numérico de suas raízes.
3.8. Solução numérica de uma equação
Serão usados dois métodos diferentes. Os valores aproximados das raízes da função podem ser vistos
na Figura 3.8.1.
(%i35) eq6:sin(x)*%e^x+x-3;
plot2d(eq6,[x,-4,4]) /* visualizando os valores aproximados
das raizes */;
find_root(eq6,x,-1,2) /* necessita um intervalo */;
find_root(eq6,x,2,4);
load(mnewton) /* carrega a funcao mnewton */;
mnewton(eq6,x,1) /* necessita de uma aproximacao */;
(%o35) ex sin (x) + x− 3
(%o36) /home/adilson/maxout.gnuplot_pipes
(%o37) 0.93839691426479
(%o38) 3.147945514659651
(%o39) /usr/share/maxima/5.32.1/share/mnewton/mnewton.mac
(%o40) [[x = 0.93839691426479]]
25
3. Breve apresentação do software livre Maxima
3.9. Manipulação algébrica de equações e funções
(%i41) eq7:x^3+4=(6+x^2)/x;
lhs(eq7) /* left hand side */;
rhs(eq7) /* right hand side */;
(x/4)*eq7;
expand(%);
eq7, x:2 /* atribuindo momentaneamente um valor para x */;
f(x,y):=x+(x/y) /* definindo uma funcao */;
f(9,3);
subst(y+z,x,eq7) /* subst y+z em x na eq7 */;
ev(eq7,x=y+z) /* avalia e eq7 com x = y+z */;
eq8:x^2+y=1;
eq9:x-2*y=2;
solve([eq8,eq9],[x,y]) /* sol analitica duas equacoes */;
float(%);
(%o41) x3 + 4 =
x2 + 6
x
(%o42) x3 + 4
(%o43)
x2 + 6
x
(%o44)
x
(
x3 + 4
)
4
=
x2 + 6
4
(%o45)
x4
4
+ x =
x2
4
+
3
2
(%o46) 12 = 5
(%o47) f (x, y) := x+
x
y
(%o48) 12
(%o49) (z + y)3 + 4 =
(z + y)2 + 6
z + y
(%o50) (z + y)3 + 4 =
(z + y)2 + 6
z + y
(%o51) y + x2 = 1
(%o52) x− 2 y = 2
(%o53) [[x = −
√
33 + 1
4
, y = −
√
3
√
11 + 9
8
], [x =
√
33− 1
4
, y =
√
3
√
11− 9
8
]]
(%o54) [[x = −1.686140661634507, y = −1.843070330817253], [x = 1.186140661634507, y = −0.40692966918274]]
26
3. Breve apresentação do software livre Maxima
3.10. Simplificações diversas
(%i55) (eq8+1/x)*eq9;
ratsimp(%);
eq10:(%e^x-1)/(1+%e^(x/2));
radcan(%) /* Simplifica expr que contem logs, exponenciais e radicais */;
radcan((log(x+x^2)-log(x))^a/log(1+x)^(a/2));
eq11:x^2+2*x*y+y^2;
factor(eq11) /* fatoriza uma expressao */;
expand(%);
log(a^b), logexpand=super;
log(a*b), logexpand=super;
log(a/b), logexpand=super;
logcontract(%);
eq12:sin(x)^2+cos(x)^2;
trigsimp(eq12) /* simplificacao trigonometrica */;
trigexpand(sin(2*x));
trigreduce(-sin(x)^2+3*cos(x)^2+x);(%o55) (x− 2 y)
(
y + x2 +
1
x
)
= 2
(
1
x
+ 1
)
(%o56) − 2x y
2 +
(
2x3 − x2 + 2) y − x4 − x
x
=
2x+ 2
x
(%o57)
ex − 1
e
x
2 + 1
(%o58) e
x
2 − 1
(%o59) log (x+ 1)
a
2
(%o60) y2 + 2x y + x2
(%o61) (y + x)2
(%o62) y2 + 2x y + x2
(%o63) log (a) b
(%o64) log (b) + log (a)
(%o65) log (a)− log (b)
(%o66) log
(a
b
)
(%o67) sin (x)2 + cos (x)2
(%o68) 1
(%o69) 2 cos (x) sin (x)
(%o70)
cos (2x)
2
+ 3
(
cos (2x)
2
+
1
2
)
+ x− 1
2
3.11. Gráficos bi e tridimensionais
Há duas maneiras de plotar funções no Maxima. A primeira, com a função plot2d, o gráfico é mostrado
numa janela separada. Caso o usuário queira que o gráfico fique embutido na planilha do Maxima,
deve usar a função wxplot2d. A sintaxe de ambas é idêntica:
27
3. Breve apresentação do software livre Maxima
ˆ plot2d(fun(x),[x,a,b]) - plota o gráfico de fun(x), sendo �x� a variável independente, de �a�
a �b�
ˆ wxplot2d(fun(x),[x,a,b])
Há argumentos opcionais, como mudar o nome dos eixos, legenda etc. Os resultados estão mostrados
nas Figuras 3.11.1, 3.11.2 e 3.11.3.
Os gráficos gerados no Maxima podem ser salvos diretamente como arquivos, nos formato png (opção
png_file), eps (opção ps_file) e pdf (opção pdf_file). O caminho onde o arquivo foi gravado
aparece na saída do comando:
(%i14) plot2d (3*x^3+5*x^2-x+6, [x,-3,1], [png_file, "funcao1.png"]);
(%o14) /home/adilson/maxout.gnuplot/home/adilson/funcao1.png
3.12. Limites e continuidade de funções
Em um comando do Maxima, ao se colocar um apóstrofo (') antes de tal comando, isto inibe o Maxima
de processar o comando, retornando a própria entrada. ATENÇÃO: ao se gerar o documento pdf, o
apóstrofo aparece como uma aspa simples. O caracter correto é o mostrado no código a seguir, colado
aqui como figura:
.
Comparar com:
(%i1) 'limit(sin(x)/x,x,0)=limit(sin(x)/x,x,0);
(%o1) lim
x→0
sin (x)
x
= 1
Tal funcionalidade é útil para se visualizar o que se deseja calcular, antes do resultado do cálculo, como
mostrado nos comandos a seguir. O gráfico da função
1
x−1 está mostrado na Figura 3.12.1.
(%i54) 'limit((x^2+x+1)/(x+1),x,1)=limit((x^2+x+1)/(x+1),x,1);
'limit(1/(x-1),x,1,plus)=limit(1/(x-1),x,1,plus);
'limit(1/(x-1),x,1,minus)=limit(1/(x-1),x,1,minus);
plot2d(1/(x-1),[x,0.5,1.5],[y,-100,100]);
'limit(1/(x-1),x,minf)=limit(1/(x-1),x,minf);
'limit(1/(x-1),x,inf)=limit(1/(x-1),x,inf);
'limit(sin(x)/x,x,0)=limit(sin(x)/x,x,0);
(%o54) lim
x→1
x2 + x+ 1
x+ 1
=
3
2
(%o55) lim
x→1+
1
x− 1 =∞
(%o56) lim
x→1−
1
x− 1 = −∞
(%o58) lim
x→−∞
1
x− 1 = 0
(%o59) lim
x→∞
1
x− 1 = 0
(%o60) lim
x→0
sin (x)
x
= 1
28
3. Breve apresentação do software livre Maxima
Figura 3.11.1.: Gráficos bidimensionais no Maxima.
29
3. Breve apresentação do software livre Maxima
(%i8) xy: [[10, .6], [20, .9], [30, 1.1], [40, 1.3], [50, 1.4]]$
plot2d([[discrete, xy], 2*%pi*sqrt(l/980)], [l,0,50],
[style, points, lines], [color, red, blue],
[point_type, asterisk],
[legend, "experimento", "teórico"],
[xlabel, "comprimento do pêndulo (cm)"],
[ylabel, "período (s)"])$
Figura 3.11.2.: Gráfico bidimensional com pontos discretos e curva contínua.
(%i6) plot3d (u^2 - v^2, [u, -2, 2], [v, -3, 3], [grid, 100, 100],
[mesh_lines_color,false])$
(%i7) plot3d (log (x^2*y^2), [x, -2, 2], [y, -2, 2],[grid, 29, 29],
[palette, [gradient, red, orange, yellow, green]], color_bar,
[xtics, 1], [ytics, 1], [ztics, 4],[color_bar_tics, 4])$
Figura 3.11.3.: Gráficos tridimensionais no Maxima.
30
3. Breve apresentação do software livre Maxima
Figura 3.12.1.: Limites simples e laterais no Maxima.
3.13. Derivadas
3.13.1. Derivadas simples e pontos críticos de uma função
O gráfico da função usada no exemplo a seguir está mostrado na Figura 3.13.1. Lembrar que os pontos
críticos de uma função são os zeros da derivada primeira. O sinal da derivada segunda, avaliada nos
pontos críticos, determinam se tais pontos são máximo, mínimo ou de inflexão.
(%i15)
f1:x^3-2.5*x^2+0.5*x+1.0;
f2:'diff(f1,x)=diff(f1,x) /* derivada primeira */;
f2a:diff(f1,x)$
f2b:f2a,x:2 /* derivada no ponto x=2 */;
f3:'diff(f1,x,2)=diff(f1,x,2) /* derivada segunda */;
f4:float(solve(f2a,x)) /* pontos criticos da funcao f1 */;
f4a:diff(f1,x,2)$
f4b:f4a,x:0.10685017607655 /* analisando se é ponto de max ou mín */;
f4c:f4a,x:1.559816490590112 /* analisando se é ponto de max ou mín */;
f5:solve(f1,x) /* raizes da funcao f1 */;
plot2d(f1,[x,-1,2.5]);
(%o15) x3 − 2.5x2 + 0.5x+ 1.0
(%o16)
d
d x
(
x3 − 2.5x2 + 0.5x+ 1.0) = 3x2 − 5.0x+ 0.5
(%o18) 2.5
(%o19)
d2
d x2
(
x3 − 2.5x2 + 0.5x+ 1.0) = 6x− 5.0
(%o20) [x = 0.10685017607655, x = 1.559816490590112]
(%o22) − 4.3588989435407
(%o23) 4.358898943540671
(%o24) [x = 2, x = −1
2
, x = 1]
3.13.2. Derivadas parciais de uma função
diff(f(x1 , x2 · · ·), x1, n1, x2, n2, · · ·): derivada parcial da função f (x1, x2, · · ·) em
relação à x1 , n1 vezes, em relação a x2 , n2 vezes,...
31
3. Breve apresentação do software livre Maxima
Figura 3.13.1.: Função usada no exemplo de derivadas e pontos críticos.
Como exemplo, deseja-se calcular as derivadas
∂3f
∂2x∂y
e
∂2f
∂x∂y
, assim como a matriz jacobiana e a
matriz hessiana, para a função f(x, y) = sen(x) cos(y2).
(%i1) f(x,y):=sin(x)*cos(y^2);
diff(f(x,y),x,2,y,1);
diff(f(x,y),x,1,y,1);
jacobian([f(x,y)],[x,y]) /* matriz jacobiana */;
jacobian([exp(x+y),x^2-y^2],[x,y])
/* matriz jacobiana de um campo vetorial */;
hessian(f(x,y),[x,y]) /* matriz hessiana */;
(%o1) f (x, y) := sin (x) cos
(
y2
)
(%o2) 2 sin (x) y sin
(
y2
)
(%o3) − 2 cos (x) y sin (y2)
(%o4)
(
cos (x) cos
(
y2
) −2 sin (x) y sin (y2))
(%o5)
(
ey+x ey+x
2x −2 y
)
(%o6)
( −sin (x) cos (y2) −2 cos (x) y sin (y2)
−2 cos (x) y sin (y2) −2 sin (x) sin (y2)− 4 sin (x) y2 cos (y2)
)
3.14. Integrais indefinidas, definidas e área abaixo de uma curva
A função f(x) = sen5(x) cos(x), usada no exemplo a seguir, cujas raízes são 0 e pi2 , está mostrada na
Figura 3.14.1.
(%i1) f6:sin(x)^5*cos(x);
f7:'integrate(f6,x)=integrate(f6,x) /* integral indefinida */;
diff(f7,x) /* checando a solução */;
solve(f6,x) /* calculando as raízes de f6 */;
f8:'integrate(f6,x,0,%pi/2)=integrate(f6,x,0,%pi/2)
/* integral definida ou a área abaixo da curva f6, entre as raízes */;
plot2d(f6,[x,0,%pi/2])$
(%o1) cos (x) sin (x)5
32
3. Breve apresentação do software livre Maxima
Figura 3.14.1.: Função usada no exemplo de integração.
(%o2)
∫
cos (x) sin (x)5dx =
sin (x)6
6
(%o3) cos (x) sin (x)5 = cos (x) sin (x)5
(%o4) [x = 0, x =
pi
2
]
(%o5)
∫ pi
2
0
cos (x) sin (x)5dx =
1
6
O comando assume serve para informar o Maxima a condição de uma constante, quando o sinal desta
afeta o cálculo a ser feito; caso contrário, o Maxima perguntará ao usuário se a constante é positiva,
negativa ou zero. Neste caso, basta digitar �p�, �n� ou �z� e em seguida SHIFT+ENTER.
(%i1) kill(all);
'integrate (x^a*exp(-x), x, 0, inf)=integrate (x^a*exp(-x), x, 0, inf);
(%o0) done
Is a+1 positive, negative or zero? p;
(%o1)
∫ ∞
0
xa e−xdx = γ (a+ 1)
(%i8) assume (a > 1)$
'integrate (x^a*exp(-x), x, 0, inf)=integrate (x^a*exp(-x), x, 0, inf);
(%o9)
∫ ∞
0
xa e−xdx = γ (a+ 1)
Quando não é possível calcular a integral definida analiticamente
2
, pode-se utilizar métodos numéricos
de integração(quadpack e romberg). O método quadpack retorna uma lista com quatro elementos:
1. Valor da integral;
2. Erro absoluto estimado;
3. Número de avaliações da função integranda;
4. Código de erro (ver manual ou o Help do Maxima para a descrição doserros, sendo 0 se não foi
encontrado nenhum problema);
2
Quando isto acontece o Maxima retorna o mesmo comando que foi fornecido.
33
3. Breve apresentação do software livre Maxima
(%i7) f9:sin(1/(x+1/10))$
f10:integrate(f9,x,0,1);
quad_qags(f9, x, 0, 1);
romberg(f9, x, 0, 1);
(%o8)
∫ 1
0
sin
(
1
x+ 110
)
dx
(%o9) [0.59450627843597, 1.9863236970333475 10−10, 147, 0]
(%o10) 0.59450623879743
3.15. Programação no Maxima
3.15.1. Instrução: if
Tomada de uma decisão numa situação condicional:
if condição then instrução;
end
if condição then instrução1 else instrução2;
end
(%i1) a:5$
if a>0 then print("número positivo")$
b:-4$
if b>0 then print("número positivo") else print("número negativo")$
número positivo
número negativo
3.15.2. Instrução: for...end
Há duas possibilidades para o uso do for :
for variável:valor1 step valor2 thru valor3 do instrução;
ou
for variável:valor1 step valor2 while condição do instrução;
(%i1) /* calcula o quadrado dos números ímpares no intervalo de 1 até 10 */
for i:1 step 2 thru 10 do print(i^2);
1
9
25
49
81
(%o1) done
(%i2) for i:1 step 2 while i<10 do print(i^2);
1
34
3. Breve apresentação do software livre Maxima
9
25
49
81
(%o2) done
3.15.3. Instrução block
Para definir procedimentos mais elaborados que englobem sequências de instruções, deve-se utilizar a
instrução block, cuja estrutura é:
nome_do_comando(argumento1, argumento2,...,argumentok):=block(
[variavel_local1, variavel_local2,...,variavel_localk],
instrução 1,
instrução 2,
...
instrução n);
(%i3) circulo(r):=block([area, perimetro],
area:%pi*r^2,
perimetro:2*%pi*r,
return([area, perimetro])
);
circulo(3);
(%o3) circulo (r) := block
(
[area, perimetro], area : pi r2, perimetro : 2pi r, return ([area, perimetro])
)
(%o4) [9pi, 6pi]
3.16. Séries
3.16.1. Série de Taylor
A função taylor(expr, x, a, n) expande a expressão expr em uma série de Taylor ou de Laurent,
truncada, na variável x em torno do ponto a, contendo termos até (x - a)^n. No exemplo a seguir,
encontram-se séries para a função sen(x) em torno do ponto zero, a forma geral da série de potência
sem truncamento e novamente a série truncada expandida em torno do ponto pi.
(%i1) f1:sin(x);
f2:taylor(f1, x, 0, 8);
f3:niceindices(powerseries(f1, x, 0));
f4:taylor(f1, x, %pi, 8);
(%o1) sin (x)
(%o2)/T/x− x
3
6
+
x5
120
− x
7
5040
+ ...
(%o3)
∞∑
i=0
(−1)i x2 i+1
(2 i+ 1)!
(%o4)/T/ − (x− pi) + (x− pi)
3
6
− (x− pi)
5
120
+
(x− pi)7
5040
+ ...
35
3. Breve apresentação do software livre Maxima
Figura 3.16.1.: Aproximações por série de Taylor para a função cos(x).
Dada a função cos(x), a Figura 3.16.1 mostra a função e os polinômios numa vizinhança de 0 (ponto
de aproximação). Percebe-se que a precisão da aproximação depende da distância do ponto de apro-
ximação (quanto mais próximo deste, melhor) e do grau do polinômio usado (quanto maior o grau,
melhor).
(%i1) f1:cos(x);
f5:taylor(f1,x,0,1);
f6:taylor(f1,x,0,4);
f7:taylor(f1,x,0,6);
f8:taylor(f1,x,0,8);
plot2d([f1,f5,f6,f7,f8],[x,-4,4],[style,[lines,1,2],[lines,1,3],
[lines,1,4],[lines,1,5],[lines,1,6]],[legend,"Função cos(x)",
"Polinômio de Taylor de grau 1", "Polinômio de Taylor de grau 4",
"Polinômio de Taylor de grau 6", "Polinômio de Taylor de grau 8"]);
(%o1) cos (x)
(%o2)/T/1 + ...
(%o3)/T/1− x
2
2
+
x4
24
+ ...
(%o4)/T/1− x
2
2
+
x4
24
− x
6
720
+ ...
(%o5)/T/1− x
2
2
+
x4
24
− x
6
720
+
x8
40320
+ ...
3.16.2. Série de Fourier
Seja f uma função periódica, com período 2L, e contínua no intervalo [-L,L]. Então a série de Fourier
de f é a série de funções:
a0 +
∞∑
n=1
[
an cos
(npix
L
)
+ bnsen
(npix
L
)]
36
3. Breve apresentação do software livre Maxima
Figura 3.16.2.: Série de Fourier da função usada no exemplo.
Os coeficientes an e bn são chamados de coeficientes de Fourier e são calculados do seguinte modo:
a0 =
1
2L
∫ L
−L
f(x)dx
an =
1
L
∫ L
−L
f(x)cos
(npix
L
)
dx, n = 1, 2, 3, ...
bn =
1
L
∫ L
−L
f(x)sen
(npix
L
)
dx, n = 1, 2, 3, ...
As funções que serão usadas no Maxima são:
ˆ fourier (f, x, p): retorna uma lista dos coeficientes de Fourier de f(x) definidos no intervalo
[-p, p].
ˆ fourexpand (l, x, p, limit): constroi e retorna a série de Fourier partindo da lista l dos
coeficientes de Fourier até ao termo limit (limit pode ser inf). x e p tem o mesmo significado
que na função fourier.
(%i1) load(fourie)$
f(x):=abs(x)-2$
g:fourexpand(fourier(f(x),x,1),x,1,3);
plot2d([f(x),g],[x,-1.5,1.5]);
(%t3) a0 = −3
2
(%t4) an = 2
(
−sin (pi n)
pi n
+
cos (pi n)
pi2 n2
− 1
pi2 n2
)
(%t5) bn = 0
(%o5) − 4 cos (3pi x)
9pi2
− 4 cos (pi x)
pi2
− 3
2
3.17. Regressão Linear
Pretende-se estudar a relação da produção de uma variedade de trigo, em kg, e a quantidade de
adubo fertilizante utilizado, em kg, conforme mostrado na Tabela 3.1. O gráfico comparando o modelo
estimado e os dados experimentais está mostrado na Figura 3.17.1.
37
3. Breve apresentação do software livre Maxima
Tabela 3.1.: Dados exeperimentais usados no exemplo de regressão linear.
Fertilizante, kg 16 24 32 40 48 56 64 72 80
Quantidade de trigo, kg 199 214 230 248 255 305 298 323 359
Figura 3.17.1.: Exercício de aplicação de regressão linear.
(%i11) dados:[[16,199],[24,214],[32,230],[40,248],[48,255],[56,305],
[64,298],[72,323],[80,359]]$
load(stats)$
simple_linear_regression(dados);
f(x):= 2.416666666666668*x + 154.111111111111;
plot2d([[discrete,dados],f(x)],[x,16,84],[style,[points,1,5],[lines,1,4]],
[legend,"Dados experimentais","Modelo calculado"],
[xlabel,"Fertilizante, kg"],[ylabel,"Produção, kg"]);
(%o13)

SIMPLE LINEAR REGRESSION
model = 2.416666666666668x+ 154.111111111111
correlation = 0.98266512147706
v_estimation = 114.0317460317458
b_conf_int = [2.009183744425169, 2.824149588908167]
hypotheses = H0 : b = 0, H1 : b#0
statistic = 14.0239216561971
distribution = [student_t, 7]
p_value = 2.220330900915002 10−6

(%o14) f (x) := 2.416666666666668x+ 154.111111111111
38
4. Campo de direções no Maxima
O campo de direções é uma ferramenta importante para a avaliação de equações diferenciais pois
dão fornecem uma solução qualitativa da equação diferencial (ou de um sistema de duas equações
diferenciais de primeira ordem autônomas) através dos valores para os quais a solução de uma dada
ED tende. O campo das direções é composto por vetores tangentes às soluções nos pontos do plano
xy.
Considere a EDO a seguir:
dv
dt
= g −
( γ
m
)
v
com: g = 9,8 m/s
2
; γ = 2 kg/s; m = 10 kg. Fixando-se valores de t e v, podemos calcular dv/dt; como
a derivada é a inclinação da reta tangente à curva v(t), os valores de dv/dt calculados são vetores que
calculados ao longo de um dado intervalo, indicam o comportamento da solução v(t). Neste exemplo,
como a EDO não depende da variável �t�, independente do valor desta, os valores da derivada somente
variam com �v�, como mostrado na Figura 4.0.1.
(%i2) kill(all);
dvdt(t,v):=g-gamma/m*v;
v:[40, 44, 48, 52, 56, 60]$
t:0$
dvdt(t,v), g:9.8, gamma:2, m:10;
t:2$
dvdt(t,v), g:9.8, gamma:2, m:10;
t:4$
dvdt(t,v), g:9.8, gamma:2, m:10;
kill(all);
plotdf(9.8-2/10*v,[t,v],[t,-0.9,10],[v,-0.9,65]);
(%o0) done
(%o1) dvdt (t, v) := g − γ
m
v
(%o4) [1.8, 1.0, 0.2,−0.59999999999999,−1.399999999999998,−2.199999999999999]
(%o6) [1.8, 1.0, 0.2,−0.59999999999999,−1.399999999999998,−2.199999999999999]
(%o8) [1.8, 1.0, 0.2,−0.59999999999999,−1.399999999999998,−2.199999999999999](%o0) done
A função que desenha o campo de direções no Maxima é plotdf, sendo que tal função só desenha
campos de direções para equações diferenciais ordinárias de primeira ordem e que seja possível explicitar
a derivada:
dy
dx
= f(x, y)
Possibilidades de uso da função plotdf:
ˆ plotdf (dydx, ...options...)
ˆ plotdf (dvdu, [u,v], ...options...)
39
4. Campo de direções no Maxima
Figura 4.0.1.: Campo de direções como vetores tangentes à solução de uma EDO.
ˆ plotdf ([dxdt,dydt], ...options...)
ˆ plotdf ([dudt,dvdt], [u,v], ...options...)
dydx, dxdt e dydt são expressões que dependem de x e y. dvdu, dudt e dvdt são expressões que
dependem de u e v. Além dessas duas variáveis, as expressões podem também depender de um conjunto
de parâmetros, com valores numéricos fornecidos com a opção de parâmetros (a sintaxe de opção é
dada abaixo), ou com um intervalo de valores permitidos especificados por uma opção �sliders�. Muitas
outras opções podem ser fornecidas dentro do comando, ou selecionadas no menu. Curvas integrais
podem ser obtidas por meio de um clique no gráfico, ou com a opção trajectory_at. A direção da
integração pode ser controlada com a opção direction, que pode ter valores de �forward� (adiante),
�backward� (para trás) ou �both� (ambos). O número de passos de integração é fornecido por meio de
nsteps e o intervalo de tempo entre eles é escolhido com a opção tstep. O método numérico utilizado
é Runge-Kutta de 4ª ordem com intervalos de tempo variáveis.
4.1. Exemplo 1
Campo de direções de
dy
dx = exp(−x) + y, mostrando a solução particular que passa pelo ponto (2;
-0,1), na Figura 4.1.1.
(%i1) plotdf(exp(-x)+y,[trajectory_at,2,-0.1])$
4.2. Exemplo 2
O seguinte exemplo mostra o campo de direção de um oscilador harmônico, definido pelas duas equações
40
4. Campo de direções no Maxima
Figura 4.1.1.: Campo de direções: exemplo 1.
dx
dt
= y e
dy
dt
=
(−k
m
)
x
e a curva integral em todo o intervalo (x,y) = (6,0), com um botão de deslizamento que irá permitir a
você mudar o valor de m interativamente (k está fixado em 2). O resultado está mostrado na Figura
4.2.1.
(%i11) plotdf([v,-k*z/m], [z,v], [parameters,"m=2,k=2"],
[sliders,"m=1:5"], [trajectory_at,6,0])$
4.3. Exemplo 3
Para montar o gráfico do campo de direção da equação de Duffing,
m
d2x
dt2
+ c
dx
dt
+ kx+ bx3 = 0
introduzimos a variável y=x' e usamos o código a seguir. O resultado é o que está mostrado na Figura
4.3.1.
(%i13) plotdf([y,-(k*x + c*y + b*x^3)/m],
[parameters,"k=-1,m=1.0,c=0,b=1"],
[sliders,"k=-2:2,m=-1:1"],[tstep,0.1]);
41
4. Campo de direções no Maxima
Figura 4.2.1.: Campo de direções: exemplo 2.
Figura 4.3.1.: Campo de direções: exemplo 3.
42
4. Campo de direções no Maxima
Figura 4.4.1.: Campo de direções: exemplo 4.
4.4. Exemplo 4
O campo de direção para um pêndulo amortecido, descrito por:
da
dt
= w
dw
dt
=
(−g
l
)
sen(a)−
(
b
ml
)
w
incluindo a solução para as condições iniciais fornecidas, com um botão de deslizamento que pode ser
usado para mudar o valor da massa m, e com um gráfico das duas variáveis de estado como uma função
do tempo, está no exemplo a seguir. Os resultados estão mostrados na Figura 4.4.1.
(%i17) plotdf([w,-g*sin(a)/l - b*w/m/l], [a,w],
[parameters,"g=9.8,l=0.5,m=0.3,b=0.05"],
[trajectory_at,1.05,-9],[tstep,0.01],
[a,-10,2], [w,-14,14], [direction,forward],
[nsteps,300], [sliders,"m=0.1:1"], [versus_t,1])$
4.5. Função drawdf
Existe também no Maxima a função drawdf que é quase idêntica ao plotdf, porém usa uma biblioteca
diferente para plotar o gráfico resultante. Usuários de wxMaxima ou Imaxima podem opcionalmente
usar wxdrawdf, que é idêntica ao drawdf, exceto que os gráficos são desenhados embutidos no �works-
pace� do Maxima, usando wxdraw. Para detalhes completos e exemplos de sua utilização, consultar o
Help do Maxima. Antes de utilizar a função, ela deve ser carregada, como mostrado na Figura 4.5.1.
43
4. Campo de direções no Maxima
Figura 4.5.1.: Campo de direções com a função drawdf.
44
5. Solução analítica de equações diferenciais
ordinárias (EDOs)
5.1. Verificando se uma dada equação é solução de uma EDO
Muitas vezes nos deparamos com uma equação e queremos simplesmente verificar se ela satisfaz uma
dada EDO, sem nos preocuparmos com o método de solução que levou à equação dada. O comando
ev serve para avaliar uma expressão na qual parte da mesma vai ser substituída por uma outra; já
o comando diff, colocando após a vírgula, instrui ao Maxima para realizar a derivada da expressão
montada com o comando anterior. Veja o exemplo a seguir.
Considere a equação de Clairaut:
y′(x)2 − xy′(x) + y(x) = 0
Verifique se y(x) = Cx− C2 e y(x) = x24 são soluções da equação de Clairaut.
(%i1) eq4: 'diff(y,x)^2 - x*'diff(y,x) + y = 0;
eq5: y = C*x - C^2;
ev(eq4,eq5);
ev(eq4,eq5),diff;
eq6: y = x^2/4;
ev(eq4,eq6);
ev(eq4,eq6),diff;
(%o1)
(
d
d x
y
)2
− x
(
d
d x
y
)
+ y = 0
(%o2) y = xC − C2
(%o3)
(
d
d x
(
xC − C2))2 − x ( d
d x
(
xC − C2))− C2 + xC = 0
(%o4) 0 = 0
(%o5) y =
x2
4
(%o6)
(
d
d x
x2
4
)2
− x
(
d
d x
x2
4
)
+
x2
4
= 0
(%o7) 0 = 0
Percebemos que ambas as soluções satisfazem a equação de Clairaut, já que após a substituição de
cada solução na EDO, econtrou-se uma identidade matemática.
5.2. Equações diferenciais de 1ª ordem
5.2.1. Encontrando a solução geral de uma ED de 1ª ordem
A solução geral de uma ED de 1ª ordem é encontrada usando as seguintes etapas:
45
5. Solução analítica de equações diferenciais ordinárias (EDOs)
1. Entrar com a EDO; nesta etapa, usar o apóstrofo antes do comando diff para inibir o Maxima
de proceder à derivada e assim permitir montar a equação diferencial;
2. Resolver a EDO usando a função ode2(edo, variável dependente, variável independente);
Se por alguma razão ode2 não conseguir encontrar a solução, retornará false, após talvez mostrar
uma mensagem de erro. Os métodos implementados para equações diferenciais de primeira ordem, na
ordem em que serão testados, são: linear, separável, exata - talvez requerendo um fator de integração,
homogênea, equação de Bernoulli, homogênea generalizada.
Durante o processo de resolução da EDO, serão dados valores a várias variáveis locais, com fins pu-
ramente informativos: method denota o método de solução usado (por exemplo, linear), intfactor
denota qualquer fator integrante utilizado, odeindex denota o índice para o método de Bernoulli ou
para o método homogêneo generalizado, e yp denota a solução particular no método de variação dos
parâmetros (solução de EDOs de 2ª ordem, não homogêneas).
Na solução geral aparece a constante de integração, que no caso de ED de 1ª ordem é nomeada
automaticamente pelo Maxima como %c. No exemplo a seguir, resolve-se a EDO de 1ª ordem, linear,
não homogênea, de coeficientes variáveis, ty′ + 2y = 4t2.
(%i4) kill(all);
eq1:t*'diff(y,t)+2*y=4*t^2;
eq2:ode2(eq1,y,t);
method;
intfactor;
(%o0) done
(%o1) t
(
d
d t
y
)
+ 2 y = 4 t2
(%o2) y =
t4 + %c
t2
(%o3) linear
(%o4) false
ATENÇÃO: Na montagem da EDO no Maxima, o comando diff é precedido por um apóstrofo (').
Entretanto, ao se gerar o documento pdf a partir do L
A
T
E
X, o apóstrofo aparece como uma aspa simples.
Compare o mesmo código anterior com o mostrado a seguir, colado como figura do Maxima, a fim de
não se alterar a formatação!
5.2.2. Encontrando a solução particular de uma ED de 1ª ordem
Considere o PVI (Problema de Valor Inicial) a seguir:
46
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.2.1.: Solução de um PVI, com equaçãodiferencial de 1ª ordem.
ty′ + 2y = 4t2
y(1) = 2
A solução deste PVI no Maxima consiste nas seguintes etapas:
1. Entrar com a EDO;
2. Resolver a EDO usando a função ode2(edo, variável dependente, variável independente);
3. Aplicar a condição inicial usando a função ic1(solução da edo, condição inicial da variável
independente, condição inicial da variável dependente);
4. Fazer um gráfico da solução (opcional, mas auxilia a visualizar o tipo de solução encontrado);
As etapas anteriores estão no código do Maxima a seguir e a solução particular encontrada está mos-
trada na Figura 5.2.1:
(%i6) kill(all);
eq1:t*'diff(y,t)+2*y=4*t^2;
eq2:ode2(eq1,y,t);
eq3:ic1(eq2,t=1,y=2);
method;
plot2d(rhs(eq3),[t,0.1,10]);
(%o0) done
(%o1) t
(
d
d t
y
)
+ 2 y = 4 t2
(%o2) y =
t4 + %c
t2
(%o3) y =
t4 + 1
t2
(%o4) linear
5.2.3. EDO com solução implícita
Nem sempre na solução de uma EDO é possível encontrar uma resposta explícita na variável depen-
dente, como é o caso da solução desta equação diferencial separável:
dy
dx
=
x2
1− y2 . A solução no
47
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Maxima, que não difere em nada dos exemplos anteriores, está mostrada a seguir.
(%i16) kill(all);
eq1:'diff(y,x)=x^2/(1-y^2);
eq2:ode2(eq1,y,x);
eq3:ratsimp(%)*3;
method;
(%o0) done
(%o1)
d
d x
y =
x2
1− y2
(%o2) − y
3 − 3 y
3
=
x3
3
+ %c
(%o3) 3 y − y3 = x3 + 3 %c
(%o4) separable
Exercício: Como você faria um esboço da solução geral do exemplo anterior usando o Maxima?
5.2.4. Solução de EDOs exatas
Considere a EDO exata a seguir: (y cos(x) + 2xey) + (sen(x) + x2ey − 1)y′ = 0. Usando o Maxima,
mostre que se trata de uma EDO exata e em seguida resolva-a.
(%i7) kill(all);
eq1:(y*cos(x)+2*x*%e^y)+(sin(x)+x^2*%e^y-1)*'diff(y,x);
M:(y*cos(x)+2*x*%e^y);
N:(sin(x)+x^2*%e^y-1);
ratsimp(diff(M,y)-diff(N,x)) /* verificando se a EDO é exata */;
eq2:ode2(eq1,y,x);
method;
(%o0) done
(%o1)
(
x2 ey + sin (x)− 1) ( d
d x
y
)
+ 2x ey + cos (x) y
(%o2) 2x ey + cos (x) y
(%o3) x2 ey + sin (x)− 1
(%o4) 0
(%o5) x2 ey + (sin (x)− 1) y = %c
(%o6) exact
5.2.5. Exemplo 1: variação da massa de um soluto com o tempo a volume constante
Considere um tanque agitado, a volume constante, no qual há a variação da massa de um soluto de
uma solução com o tempo, conforme descrito pelo esquema mostrado na Figura 5.2.2 e no modelo a
seguir:
dm
dt
+
( r
V
)
m = rc0
t = 0; m = m0
48
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.2.2.: Variação da massa com o tempo para um tanque agitado (m0=1, c0=1, V= 2, r=1)
No modelo, m é massa dentro do tanque; r é a vazão volumétrica [volume/tempo] e c concentração
mássica [massa/volume].
A solução deste modelo, no Maxima, é (ver a Figura 5.2.2 para a solução na forma gráfica para um
conjunto específico das constantes do modelo):
(%i1) edo1:'diff(m,t)+(r/V)*m=r*c0;
sol1:ode2(edo1,m,t)$ expand(%);
sol2:ic1(sol1,t=0,m=m0)$ expand(%);
sol3:sol2, m0:1, c0:1, V: 2, r:1$ expand(%);
'limit(rhs(sol3),t,inf)=limit(rhs(sol3),t,inf);
plot2d(rhs(sol3),[t,0,15],[ylabel,"m"])$
(%o1)
mr
V
+
d
d t
m = c0 r
(%o3) m = %c e−
r t
V + c0V
(%o5) m = −c0V e− r tV +m0 e− r tV + c0V
(%o7) m = 2− e− t2
(%o8) lim
t→∞ e
− t
2
(
2 e
t
2 − 1
)
= 2
5.2.6. Exemplo 2: lei do resfriamento de Newton
Um corpo que está em uma temperatura T, exposto a um ambiente de temperatura Ta, troca calor
com tal ambiente segundo a lei do resfriamento de Newton, descrita por:
dT
dt
= −k [T − Ta]. Ta pode variar no tempo, como por exemplo: Ta(t) = T0 + T1 cos(ωt). A solução
gráfica está na Figura 5.2.3.
49
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.2.3.: Solução do problema envolvendo a lei de resfriamento de Newton.
(%i1) kill(all);
edo1:'diff(T,t)=-k*(T-Ta(t));
sol1:ode2(edo1,T,t)$ expand(%);
edo2:'diff(T,t)=-k*(T-T0+T1*cos(\omega*t));
sol2:ode2(edo2,T,t)$ expand(%);
sol3:ic1(sol2,t=0,T=20)$
sol4:sol3,k:0.2, \omega:0.5, T0:25, T1:1;
plot2d(rhs(sol4),[t,0,50],[ylabel,"T"])$
(%o0) done
(%o1)
d
d t
T = −k (T − Ta (t))
(%o3) T = k e−k t
∫
ek t Ta (t) dt+ %c e−k t
(%o4)
d
d t
T = −k (cos (ω t) T1− T0 + T )
(%o6) T = −k ω sin (ω t) T1
ω2 + k2
− k
2 cos (ω t) T1
ω2 + k2
+ T0 + %c e−k t
(%o8) T =− 3.448275862068965 e−0.2 t(
0.1 e0.2 t sin (0.5 t) + 0.04 e0.2 t cos (0.5 t) + 25
(
0.29− 0.29 e0.2 t)− 5.84)
5.3. Equações diferenciais lineares de 2ª ordem
5.3.1. Solução geral
Os tipos de equações de segunda ordem que podem ser resolvidas são: coeficientes constantes, exatas,
linear homogêneas com coeficientes não constantes que possam ser transformados para constantes,
equação de Euler ou equi-dimensional, equações que possam ser resolvidas pelo método de variação
dos parâmetros, e equações que não dependam ou da variável independente ou da variável dependente,
de modo que possam ser reduzidas a duas equações lineares de primeira ordem a serem resolvidas
sequencialmente. Algumas equações especiais também são reconhecidas. Caso o Maxima não consiga
resolver a EDO, é retornado false.
50
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.3.1.: Funções de Bessel J0(x) e Y0(x).
Considere a equação de Bessel de ordem zero:
x2y′′ + xy′ + x2y = 0
Cuja solução geral é dada nos manuais como sendo:
y(x) = c1J0(x) + c2Y0(x)
(%i1) eq1:x^2*'diff(y,x,2)+x*'diff(y,x)+x^2*y=0;
sol1:ode2(eq1,y,x);
plot2d([bessel_y(0,x), bessel_j(0,x)],[x,0,15],[y,-0.5,1])$
limit(bessel_y(0,x),x,0);
(%o1) x2
(
d2
d x2
y
)
+ x
(
d
d x
y
)
+ x2 y = 0
(%o2) y = bessel_y (0, x) %k2 + bessel_j (0, x) %k1
(%o4) −∞
Na solução anterior, J0(x) e Y0(x) são chamadas de funções de Bessel, do primeiro tipo e ordem zero.
Percebe-se na Figura 5.3.1 o comportamento oscilatório de tais funções, passando por valores negativos.
Isto leva à necessidade de adimensionalizar modelos descritos por equações de Bessel, que aparecem
em geometrias cilíndricas e esféricas, para o intervalo [0, 1], a fim de garantir que a solução obtida
seja fisicamente compatível, por exemplo, evitando de serem obtidas concentrações negativas. Outra
observação, é que a função Y0(x) tende para menos infinito quando x tende a zero. Isto faz com que
a constante de integração associada a esta função seja arbitrada como sendo nula a fim de evitar que
concentração ou temperatura vá para menos infinito no centro do corpo, ou em x = 0, em geometrias
cilíndricas e esféricas. Seus valores encontram-se em tabelas ou podem ser calculados por �softwares� e
calculadoras científicas. O Maxima dispõe de todas as funções de Bessel na classe de funções especiais.
(%i5) [bessel_j(0,1), bessel_y(0,1), bessel_k(0,1), bessel_i(0,1)], numer;
(%o5) [0.76519768655796, 0.088256964215676, 0.4210244382407, 1.266065877752008]
5.3.2. Solução de PVIs
Considere o PVI a seguir:
51
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.3.2.: PVI envolvendo uma EDO de 2ª ordem.
d2y
dx2
− 3dy
dx
+ 2y = 0
y(0) = 0 y′(0) = 0, 5
A solução encontra-se a seguir, onde a função ic2 é usada para inserir a condição inicial envolvendo
a variável dependente e sua derivada, conhecidas para um mesmo ponto da variável independente. A
solução geral encontrada é comparada com as raízes da equação característica. A solução particular
do problema está mostrada na Figura 5.3.2.
(%i10) kill(all);
eq1:'diff(y,x,2)-3*'diff(y,x)+2*y /* edo homogenea */;
eq2:ode2(eq1,y,x);
eq3:ic2(eq2,x=0,y=0,'diff(y,x)=0.5);
plot2d(rhs(eq3),[x,-5,1]);
method;
eq_caracteristica:r^2-3*r+2;
solve(eq_caracteristica,r);(%o0) done
(%o1)
d2
d x2
y − 3
(
d
d x
y
)
+ 2 y
(%o2) y = %k1 e2x + %k2 ex
rat: replaced 0.5 by 1/2 = 0.5
(%o3) y =
e2x
2
− e
x
2
(%o5) constcoeff
(%o6) r2 − 3 r + 2
(%o7) [r = 1, r = 2]
5.3.3. Exemplo: difusão e reação de um componente gasoso em um líquido (Rice and
Do., 1995, exemplo 2.7)
52
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Considere uma torre de absorção reativa, usada
na lavagem de um gás tóxico, na qual o gás solúvel
A, reativo, se dissolve numa interface plana de um
corpo longo de líquido reagente, com o qual ele
reage através da equação da taxa irreversível e
não linear:
rA = knC
n
A
Assumindo que a reação seja de primeira ordem
(n = 1) e que em z = 0, CA = 0 e
dCA
dz = 1. O ba-
lanço de massa no sistema, no estado estacinário,
conduz a:
D
d2CA
dz2
− knCA = 0. A solução no Maxima é:
(%i1) kill(all)$
assume(D>0,k_n>0);
edo1:D*'diff(C_A,z,2)-k_n*C_A^1;
solve(D*r^2-k_n,r) /* equação característica */;
sol3:ode2(edo1,C_A,z);
sol4:ic2(sol3, z=0, C_A=0,'diff(C_A,z)=1);
sol5:rhs(sol4),D=0.1,k_n=0.5;
plot2d(sol5,[z,0,1],[style,[lines,5]],[ylabel,"C_A"],
[xlabel,"z"])$
(%o1) [D > 0, k_n > 0]
(%o2)
(
d2
d z2
C_A
)
D − k_nC_A
(%o3) [r = −
√
k_n√
D
, r =
√
k_n√
D
]
(%o4) C_A = %k1 e
√
k_n z√
D + %k2 e
−
√
k_n z√
D
(%o5) C_A =
√
De
√
k_n z√
D
2
√
k_n
−
√
De
−
√
k_n z√
D
2
√
k_n
(%o6) 0.22360679774997 e2.23606797749979 z − 0.22360679774997 e−2.23606797749979 z
Adotando D = 0, 1 e kn = 0, 5; o gráfico da concentração em função do comprimento está mostrado na
Figura 5.3.3, onde percebe-se o comportamento exponencial do sistema, ditado pelas raízes da equação
característica. Este comportamento é genérico para este tipo de sistema ou é um caso particular para
os valores das constantes adotadas?
5.3.4. Solução de PVCs
O PVC:
d2y
dx2
− 3dy
dx
+ 2y = 0
y(0) = 0 y(1) = 2
tem a solução mostrada a seguir, onde a função bc2 é usada para inserir a condição de contorno
envolvendo a variável dependente e a variável independente, conhecidas para dois pontos diferentes da
variável independente:
53
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.3.3.: Difusão e reação de um componente gasoso.
(%i10) kill(all);
eq1:'diff(y,x,2)-3*'diff(y,x)+2*y /* edo homogenea */;
eq2:ode2(eq1,y,x);
eq3:bc2(eq2,x=0,y=0,x=1,y=2);
eq3,x:0 /* checando as condicoes de contorno */;
eq3,x:1, numer;
method;
eq_caracteristica:r^2-3*r+2;
solve(eq_caracteristica,r);
(%o0) done
(%o1)
d2
d x2
y − 3
(
d
d x
y
)
+ 2 y
(%o2) y = %k1 e2x + %k2 ex
(%o3) y =
2 e2x
e2 − e −
2 ex
e2 − e
(%o4) y = 0
(%o5) y = 1.999999999999999
(%o6) constcoeff
(%o7) r2 − 3 r + 2
(%o8) [r = 1, r = 2]
5.3.5. EDOs não homogêneas
Considere o PVI a seguir:
d2y
dx2
− 3dy
dx
+ 2y = 5 exp(x)− sin(x)
y(0) = 0 y′(0) = 0, 5
54
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.3.4.: Solução de uma EDO de segunda ordem, não homogênea.
A solução segue os mesmos passos já apresentados; a única novidade aqui é que o método de solução
usado pelo Maxima é a variação dos parâmetros, que fica armazenada na variável yp. A solução
particular está mostrada na Figura 5.3.4.
(%i9) kill(all);
eq1:'diff(y,x,2)-3*'diff(y,x)+2*y-5*%e^x+sin(x) /* edo nao homogenea */;
eq2:ode2(eq1,y,x);
eq3:ic2(eq2,x=0,y=0,'diff(y,x)=0.5);
plot2d(rhs(eq3),[x,-5,1],[ylabel,"y"]);
method;
eq_caracteristica:r^2-3*r+2;
solve(eq_caracteristica,r);
yp;
(%o0) done
(%o1)
d2
d x2
y − 3
(
d
d x
y
)
+ 2 y + sin (x)− 5 ex
(%o2) y = −sin (x) + 3 cos (x) + (50x+ 50) e
x
10
+ %k1 e2x + %k2 ex
rat: replaced 10.6 by 53/5 = 10.6
(%o3) y =
53 e2x
10
− sin (x) + 3 cos (x) + (50x+ 50) e
x
10
(%o5) variationofparameters
(%o6) r2 − 3 r + 2
(%o7) [r = 1, r = 2]
(%o8) − sin (x) + 3 cos (x) + (50x+ 50) e
x
10
5.3.6. Exemplo: Dinâmica de um tanque agitado com aquecimento elétrico
Um fenômeno que é amplamente estudado no Controle de Processos é o arranjo em série de sistemas,
resultando numa dinâmica global de ordem superior. Por exemplo, dois sistemas de primeira
55
5. Solução analítica de equações diferenciais ordinárias (EDOs)
ordem em série, resultam num sistema de segunda ordem. Isto pode ser facilmente visualizado
quando se toma as equações diferenciais individuais de cada sistema, de primeira ordem, e combina-se
tais equações numa única, resultando numa EDO de segunda ordem.
Considere um tanque agitado com aquecimento elétrico, como o mostrado na Figura 5.3.5. A dinâmica
da temperatura dentro do tanque e no sistema de aquecimento pode ser representada por:
mC
dT
dt
= wC(Ti − T ) + heAe(Te − T )
meCe
dTe
dt
= Q− heAe(Te − T )
A variável Te pode ser isolada na primeira EDO e substituída na segunda, resultando numa EDO de
2ª ordem, linear, de coeficientes constantes (com os parâmetros constantes), não homogênea, apenas
com variável dependente T e variável independente t. Em seguida, a EDO de 2ª ordem resultante é
resolvida, considerando m = 2, C = 1, me = 1, Ce = 0, 5; w = 1, he = 1, Ae = 1, Ti = 10 + 1, 5sent,
Te = 15, Qe = 0, 5. Como as raízes da equação característica são r1 = �2,6 e r2 = �0,38, portanto reais,
diferentes e de sinal negativo, parte da resposta deve ter um comportamento exponencial e assintótico,
o que de fato se verifica na solução encontrada (ver Figura 5.3.5c). Como a temperatura de entrada
Ti = 10 + 1, 5sent é função do tempo e oscilatória (influenciando na parte não homogênea da EDO),
isto faz com que o sistema seja continuamente excitado, embora resulte numa dinâmica final estável,
para o conjunto de parâmetros adotados.
Neste exemplo, mostra-se que um sistema de duas EDOs pode ser combinado numa única EDO de
ordem superior. A recíproca também é verdadeira, ou seja, uma EDO de ordem n pode ser decomposta
em um sistema de n EDOs de 1ª ordem.
(%i1) depends([T,Te,Ti],t);
eqn1:m*C*'diff(T,t)=w*C*(Ti-T)+he*Ae*(Te-T);
eqn2:me*Ce*'diff(Te,t)=Q-he*Ae*(Te-T);
eqn3:solve(eqn1,Te)[1];
eqn4:(Ae*he*ev(eqn2,eqn3,diff)+Ae*he*m*C*('diff(T,t,1))+
Ce*me*('diff(Ti,t,1))*w*C+Ae*he*w*C*T)/(C*w*he*Ae)$
eqn4:expand(%);
eqn5:eqn4,m:2,C:1,me:1,Ce:0.5,w:1,he:1,Ae:1,Ti:10+1.5*sin(t),Te:15,Q:0.5,diff;
eqn6:ode2(eqn5,T,t)$
eqn7:ic2(eqn6,t=0,T=5,'diff(T,t)=0);
plot2d(rhs(eqn7),[t,0,30],[xlabel,"t"],[ylabel,"T"]);
(%o1) [T (t) ,Te (t) ,Ti (t)]
(%o2) mC
(
d
d t
T
)
= wC (Ti− T ) +Aehe (Te− T )
(%o3) Ceme
(
d
d t
Te
)
= Q−Aehe (Te− T )
(%o4) Te =
mC
(
d
d t T
)
+ (wC +Aehe) T − TiwC
Aehe
(%o6)
Cemme
(
d2
d t2
T
)
Aehew
+
Ceme
(
d
d t T
)
wC
+
m
(
d
d t T
)
w
+
Ceme
(
d
d t T
)
Aehe
+T =
Q
wC
+
Ceme
(
d
d t Ti
)
Aehe
+Ti
(%o7) 1.0
(
d2
d t2
T
)
+ 3.0
(
d
d t
T
)
+ T = 1.5 sin (t) + 0.75 cos (t) + 10.5
(%o9) T =
sin (t)− 2 cos (t) + 42
4
−
(
31
√
5 + 50
)
e
(
√
5−3) t
2
20
+
(
31
√
5− 50) e (−√5−3) t2
20
56
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.3.5.: Tanque agitado com aquecimento elétrico.
5.3.7. Wronkiano
O wronkiano de uma base de soluções pode ser calculado com o comando wronskian([fun1, fun2],
var). Na verdade o comando wronskian apenas monta a matriz das soluções e suas respectivas
derivadas primeiras; depois, deve ser seguido do comando determinant para calcular o determinante
da matriz anterior:
(%i1) load(functs)$
wronskian([f(x), g(x)],x);
w:determinant(wronskian([%e^x, %e^(2*x)],x));
(%o2)
(
f (x) g (x)
d
d x f (x)
d
d x g (x))
(%o3) e3x
5.4. Pacote contrib_ode
A função ode2 do Maxima, resolve uma equação diferencial ordinária (EDO) linear, de primeira ou de
segunda ordem, como mostrado nos exemplos anteriores. O pacote contrib_ode contém, dentre outras
função, a função de mesmo nome do pacote que estende as capacidades de ode2 com métodos adicionais
para EDOs de primeira ordem lineares e não lineares e EDOs de segunda ordem lineares homogêneas.
O código ainda está em desenvolvimento e a seqüência de chamamento pode mudar em versões futuras.
Uma vez que o código tenha se estabilizado, pode ser movido a partir do diretório �contrib� e integrado
ao Maxima. Este pacote deve ser carregado com o comando load(contrib_ode) antes do seu uso.
A convenção de chamada para contrib_ode é idêntico ao ode2. Leva três argumentos: uma EDO
(só precisa de ser dado o lado esquerdo, se o lado direito é 0), a variável dependente e as variáveis
independentes. Quando bem sucedido, ele retorna uma lista de soluções. A forma da solução difere
de ode2. Como equações não lineares podem ter várias soluções, contrib_ode retorna uma lista de
soluções. Cada solução pode ter um número de formas:
ˆ uma solução explícita para a variável dependente;
ˆ uma solução implícita para a variável dependente;
ˆ uma solução paramétrica em termos de variável %t, ou uma transformação em outra EDO na
variável %u;
57
5. Solução analítica de equações diferenciais ordinárias (EDOs)
%c é usado para representar a constante de integração para equações de primeira ordem. %k1 e %k2
são as constantes para equações de segunda ordem. Se contrib_ode não pode obter uma solução
por qualquer motivo, ele retorna false, após talvez mostrar uma mensagem de erro. É necessário
retornar uma lista de soluções, já que até mesmo EDOs de primeira ordem não lineares podem ter
várias soluções. Por exemplo:
(%i1) load(contrib_ode)$
eq1:x*'diff(y,x)^2-(1+x*y)*'diff(y,x)+y=0$
eq2:contrib_ode(eq1,y,x);
(%t3) x
(
d
d x
y
)2
− (x y + 1)
(
d
d x
y
)
+ y = 0
first order equation not linear in y'
(%o3) [y = log (x) + %c, y = %c ex]
EDOs não lineares podem ter soluções singulares sem constantes de integração, como na segunda
solução do exemplo a seguir:
(%i10) load(contrib_ode)$
eq3:'diff(y,x)^2+x*'diff(y,x)-y=0$
eq4:contrib_ode(eq3,y,x);
method;
(%t11)
(
d
d x
y
)2
+ x
(
d
d x
y
)
− y = 0
first order equation not linear in y'
(%o11) [y = %c x+ %c2, y = −x
2
4
]
(%o12) clairault
A EDO seguinte tem duas soluções paramétricas em termos da variável muda %t. Neste caso, as
soluções paramétricas podem ser manipuladas para dar soluções explícitas:
(%i1) load(contrib_ode)$
eq5:'diff(y,x)=(x+y)^2$
eq6:contrib_ode(eq5,y,x);
method;
(%o3) [[x = %c− atan
(√
%t
)
, y = −x−
√
%t], [x = atan
(√
%t
)
+ %c, y =
√
%t− x]]
(%o4) lagrange
A função ode_check retorna o valor da EDO eqn após a substituição de uma possível solução soln.
O valor é equivalente a zero se soln é uma solução de eqn:
(%i1) load(contrib_ode)$
eqn:'diff(y,x,2)+(a*x+b)*y;
ans:[y = bessel_y(1/3,2*(a*x+b)^(3/2)/(3*a))*%k2*sqrt(a*x+b)
+bessel_j(1/3,2*(a*x+b)^(3/2)/(3*a))*%k1*sqrt(a*x+b)];
ode_check(eqn,ans[1]);
(%o2)
d2
d x2
y + (a x+ b) y
(%o3) [y = bessel_y
(
1
3
,
2 (a x+ b)
3
2
3 a
)
%k2
√
a x+ b+ bessel_j
(
1
3
,
2 (a x+ b)
3
2
3 a
)
%k1
√
a x+ b]
58
5. Solução analítica de equações diferenciais ordinárias (EDOs)
(%o4) 0
5.5. Solução usando série de potências para EDOs lineares de 2ª
ordem
Não há no Maxima nenhum comando que permite resolver diretamente uma EDO por série de potên-
cias. O software apenas auxilia na manipulação algébrica, incluindo os somatórios derivadas e a lei geral
da relação de recorrência, embora usuários mais avançados podem usar os recursos de programação
existentes no Maxima para automatizar o procedimento de solução.
Considere-se a solução da seguinte EDO de 2ª ordem, com coeficientes variáveis, por série de potências
(Rice and Do., 1995, exemplo 3.1):
4x y′′ + 6 y′ − y = 0
Será proposta uma solução do tipo:
y =
∞∑
n=0
anx
n+c
(%i1) kill(all);
y:sum(a[n]*x^(n+c),n,0,inf);
dy:diff(y,x);
d2y:diff(y,x,2);
resp1:4*x*d2y + 6*dy - y = 0;
(%o0) done
(%o1)
∞∑
n=0
an x
n+c
(%o2)
∞∑
n=0
(n+ c) an x
n+c−1
(%o3)
∞∑
n=0
(n+ c− 1) (n+ c) an xn+c−2
(%o4) −
( ∞∑
n=0
an x
n+c
)
+ 6
( ∞∑
n=0
(n+ c) an x
n+c−1
)
+ 4x
∞∑
n=0
(n+ c− 1) (n+ c) an xn+c−2 = 0
Ajustando as potências de x para (n+c):
(%i5) resp2: -sum(a[n-1]*x^(n+c-1),n,1,inf)+
6*dy+4*sum((n+c-1)*(n+c)*a[n]*x^(n+c-1),n,0,inf)=0;
(%o5) 4
( ∞∑
n=0
(n+ c− 1) (n+ c) an xn+c−1
)
+ 6
( ∞∑
n=0
(n+ c) an x
n+c−1
)
−
∞∑
n=1
an−1 xn+c−1 = 0
Abrindo os dois primeiros somatórios:
(%i6) 4*(c-1)*c*a[0]*x^(c-1)+6*c*a[0]*x^(c-1)-
sum(a[n-1]*x^(n+c-1),n,1,inf)+6*sum((n+c)*a[n]*x^(n+c-1),n,1,inf)+
4*sum((n+c-1)*(n+c)*a[n]*x^(n+c-1),n,1,inf)=0;
(%o6) 4
( ∞∑
n=1
(n+ c− 1) (n+ c) an xn+c−1
)
+ 6
( ∞∑
n=1
(n+ c) an x
n+c−1
)
−
( ∞∑
n=1
an−1 xn+c−1
)
+
4 a0 (c− 1) c xc−1 + 6 a0 c xc−1 = 0
59
5. Solução analítica de equações diferenciais ordinárias (EDOs)
A equação indicial é:
(%i7) eqind0:4*(c-1)*c+6*c=0;
eqind1:expand(%);
eqind2:solve(eqind1,c);
(%o7) 4 (c− 1) c+ 6 c = 0
(%o8) 4 c2 + 2 c = 0
(%o9) [c = −1
2
, c = 0]
Relação de recorrência 1, com c = −12 , válida para n ≥ 1:
(%i10) relrec1a:4*(n+c-1)*(n+c)*a[n]+6*(n+c)*a[n]-a[n-1]=0, c=-1/2;
load(solve_rec);
y1:solve_rec(relrec1a,a[n]);
(%o10) 4
(
n− 3
2
) (
n− 1
2
)
an + 6
(
n− 1
2
)
an − an−1 = 0
(%o11) /usr/share/maxima/5.32.1/share/solve_rec/solve_rec.mac
(%o12) an =
√
pi%k1
22n n! Γ
(
n+ 12
)
(%i13) y1:y1,%k[1]=a[0];
(%o13) an =
√
pi a0
22n n! Γ
(
n+ 12
)
Na solução anterior, Γ(x) é a função gama. Relação de recorrência 2 com c = 0, válida para n ≥ 1:
(%i14) relrec1b:4*(n+c-1)*(n+c)*b[n]+6*(n+c)*b[n]-b[n-1]=0, c=0;
y2:solve_rec(relrec1b,b[n]);
y2:y2,%k[1]=b[0];
(%o14) 4 (n− 1) n bn + 6n bn − bn−1 = 0
(%o15) bn =
√
pi%k1 2
−2n−1
n! Γ
(
n+ 32
)
(%o16) bn =
√
pi b0 2
−2n−1
n! Γ
(
n+ 32
)
Solução geral, com a0 e b0 sendo as constantes de integração válidas para n ≥ 1 e verificando a solução,
substituindo na EDO:
(%i17) sol:sum(rhs(y1)*x^(n-1/2),n,0,inf)+sum(rhs(y2)*x^(n),n,0,inf);
dsol:diff(sol,x)$
d2sol:diff(sol,x,2)$
resp2:4*x*d2sol + 6*dsol - sol = 0;
(%o17)
√
pi b0
( ∞∑
n=0
2−2n−1 xn
n! Γ
(
n+ 32
))+√pi a0 ∞∑
n=0
xn−
1
2
22n n! Γ
(
n+ 12
)
60
5. Solução analítica de equações diferenciais ordinárias (EDOs)
(%o20) −√pi b0
( ∞∑
n=0
2−2n−1 xn
n! Γ
(
n+ 32
))−√pi a0 ( ∞∑
n=0
xn−
1
2
22n n! Γ
(
n+ 12
))+
6
(
√
pi b0
( ∞∑
n=0
n 2−2n−1 xn−1
n! Γ
(
n+ 32
) )+√pi a0 ∞∑
n=0
(
n− 12
)
xn−
3
2
22n n! Γ
(
n+ 12
))+
4x
(
√
pi b0
( ∞∑
n=0
(n− 1) n 2−2n−1 xn−2
n! Γ
(
n+ 32
) )+√pi a0 ∞∑
n=0
(
n− 32
) (
n− 12
)
xn−
5
2
22n n! Γ
(
n+ 12
) ) = 0
Expandindo a solução para 10 termos e verificando a solução expandida na EDO:
(%i21) Y:sum(rhs(y1)*x^(n-1/2),n,1,10)+sum(rhs(y2)*x^(n),n,1,10);
dY:diff(Y,x)$
d2Y:diff(Y,x,2)$
resp3:4*x*d2Y + 6*dY - Y = 0$
ratsimp(expand(resp3)),x=1,a[0]=-1,b[0]=1,numer;
(%o21)
b0 x
10
51090942171709440000
+
a0 x
19
2
2432902008176640000
+
b0 x
9
121645100408832000
+
a0 x
17
2
6402373705728000
+
b0 x
8
355687428096000
+
a0 x
15
2
20922789888000
+
b0 x
7
1307674368000
+
a0 x
13
2
87178291200
+
b0 x
6
6227020800
+a0 x
11
2
479001600
+
b0 x
5
39916800
+
a0 x
9
2
3628800
+
b0 x
4
362880
+
a0 x
7
2
40320
+
b0 x
3
5040
+
a0 x
5
2
720
+
b0 x
2
120
+
a0 x
3
2
24
+
b0 x
6
+
a0
√
x
2
rat: replaced -7.829176425356505E-18 by -1/127727355429273600 = -7.829176425356505E-18
(%o25) − 7.829176425356505 10−18 = 0
Como pode-se considerar que ambos os lados da última igualdade são iguais a zero, a solução consi-
derada satisfaz a EDO. Lembrar que foram usados apenas 10 termos no somatório. Comparando a
solução obtida no Maxima com a solução apresentada pelo Rice and Do. (1995):
(%i26) sol_rice:a[0]*cosh(sqrt(x))/sqrt(x)+b[0]*sinh(sqrt(x))/sqrt(x);
Y, x=1, a[0]=-1, b[0]=1, numer;
sol_rice, x=1, a[0]=-1, b[0]=1, numer;
(%o26)
b0 sinh (
√
x)√
x
+
a0 cosh (
√
x)√
x
(%o27) − 0.36787944117144
(%o28) − 0.36787944117144
Como pode ser percebido, ambas as soluções tem o mesmo valor, para x = 1. Fazendo o teste da
indepêndia linear da solução:
(%i29) f:a[0]*cosh(sqrt(x))/sqrt(x);
g:b[0]*sinh(sqrt(x))/sqrt(x);
load(functs);
w1:determinant(wronskian([f,g],x));
w2:ratsimp(w1);
w3:trigsimp(w2);
(%o29)
a0 cosh (
√
x)√
x
(%o30)
b0 sinh (
√
x)√
x
61
5. Solução analítica de equações diferenciais ordinárias (EDOs)
(%o31) /usr/share/maxima/5.32.1/share/simplification/functs.mac
(%o32)
a0 cosh (
√
x)
(
b0 cosh(
√
x)
2x −
b0 sinh(
√
x)
2x
3
2
)
√
x
−
b0 sinh (
√
x)
(
a0 sinh(
√
x)
2x −
a0 cosh(
√
x)
2x
3
2
)
√
x
(%o33) − a0 b0 sinh (
√
x)
2 − a0 b0 cosh (
√
x)
2
2x
3
2
(%o34)
a0 b0
2x
3
2
Como a0, b0 e x (ponto singular em x = 0) são todos diferentes de zero, o wronskiano é também
diferente de zero; portanto, as soluções são independentes.
5.6. Solução de EDOs por transformada de Laplace (EDOs simples e
sistemas de EDOs)
5.6.1. Transformada de Laplace
A função laplace (expr, t, s) tenta calcular a transformada de Laplace de expr com relação a uma
variável t e parâmetro de transformação s. Se laplace não pode achar uma solução, um substantivo
laplace é retornado. laplace reconhece em expr as funções delta, exp, log, sin, cos, sinh, cosh, e erf,
também derivative, integrate, sum, e ilt. Se algumas outras funções estiverem presente, laplace pode
não ser habilitada a calcular a tranformada.
Exemplos:
(%i1) laplace (exp (2*t + a) * sin(t) * t, t, s);
laplace ('diff (f (x), x), x, s);
diff (diff (delta (t), t), t);
laplace (%, t, s);
(%o1)
ea (2 s− 4)
(s2 − 4 s+ 5)2
(%o2) s laplace (f (x) , x, s)− f (0)
(%o3)
d2
d t2
δ (t)
(%o4) − d
d t
δ (t)
∣∣∣∣
t=0
+ s2 − δ (0) s
(%i1) assume(a>0)$
laplace(gamma_incomplete(a,t),t,s),gamma_expand:true;
Is a an integer? y;
(%o2)
Γ (a)
s
− Γ (a) s
−a−1(
1
s + 1
)a
A função ilt(expr, s, t) calcula a transformada inversa de Laplace de expr em relação a s e o
parâmetro t. expr deve ser uma razão de polinômios cujo denominador tem apenas fatores lineares
ou quadráticos. Exemplo:
(%i1) ilt( 1/s^2 ,s , t );
(%o1) t
62
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Um exemplo mais complexo é mostrado a seguir:
(%i2) 'integrate (sinh(a*x)*f(t-x), x, 0, t) + b*f(t) = t**2;
laplace (%, t, s);
linsolve ([%], ['laplace(f(t), t, s)]);
ilt (rhs (first (%)), s, t);
(%o2)
∫ t
0
f (t− x) sinh (a x) dx+ b f (t) = t2
(%o3)
a laplace (f (t) , t, s)
s2 − a2 + b laplace (f (t) , t, s) =
2
s3
(%o4) [laplace (f (t) , t, s) =
2 s2 − 2 a2
b s5 + (a− a2 b) s3 ]
Is a b (a b− 1) positive, negative or zero?p;
(%o5) −
2 cosh
(√
a b (a b−1) t
b
)
a3 b2 − 2 a2 b+ a +
a t2
a b− 1 +
2
a3 b2 − 2 a2 b+ a
Frações parciais também podem ser obtidas no Maxima com o comando partfrac(expr, var):
(%i1) expr1:-x/(x^3+4*x^2+5*x+2);
expr2:partfrac(expr1,x);
ratsimp(%) /* verificação da solução */;
(%o1) − x
x3 + 4x2 + 5x+ 2
(%o2)
2
x+ 2
− 2
x+ 1
+
1
(x+ 1)2
(%o3) − x
x3 + 4x2 + 5x+ 2
5.6.2. Solução de EDOs por transformada de Laplace
A função desolve resolve tanto uma única EDO, quanto sistemas de equações diferenciais ordinárias
lineares usando transformada de Laplace. Aqui as expressões eqn são equações diferenciais nas variáveis
dependentes x_1, ..., x_n. A relação funcional de x_1, ..., x_n na variável independente deve ser
indicada explicitamente nas variáveis e nas suas derivadas. Por exemplo, esta forma de definir as
equações não é correta:
(%i1) eqn_1: 'diff(f,x,2) = sin(x) + 'diff(g,x)$
eqn_2: 'diff(f,x) + x^2 - f = 2*'diff(g,x,2)$
A forma correta é:
(%i3) eqn_1: 'diff(f(x),x,2) = sin(x) + 'diff(g(x),x);
eqn_2: 'diff(f(x),x) + x^2 - f(x) = 2*'diff(g(x),x,2);
(%o3)
d2
d x2
f (x) =
d
d x
g (x) + sin (x)
(%o4)
d
d x
f (x)− f (x) + x2 = 2
(
d2
d x2
g (x)
)
Assim, a chamada à função desolve é:
63
5. Solução analítica de equações diferenciais ordinárias (EDOs)
(%i3) desolve([eqn_1, eqn_2], [f(x),g(x)])$
Se as condições iniciais em x=0 forem conhecidas, poderão ser fornecidas antes de usar desolve, através
de atvalue. Por exemplo, considere o sistema linear a seguir:
df(x)
dx =
dg(x)
dx + sen(x)
d2g(x)
x2
= df(x)dx − cos(x)
dg(x)
dx |x=0 = a
f(0) = 1
Sua resolução no Maxima, com o comando desolve, é:
(%i1) eqn_1:'diff(f(x),x)='diff(g(x),x)+sin(x);
eqn_2:'diff(g(x),x,2)='diff(f(x),x)-cos(x);
atvalue('diff(g(x),x),x=0,a);
atvalue(f(x),x=0,1);
eqn_3:desolve([eqn_1,eqn_2],[f(x),g(x)]);
[eqn_1,eqn_2],eqn_3,diff;
(%o1)
d
d x
f (x) =
d
d x
g (x) + sin (x)
(%o2)
d2
d x2
g (x) =
d
d x
f (x)− cos (x)
(%o3) a
(%o4) 1
(%o5) [f (x) = a ex − a+ 1, g (x) = cos (x) + a ex − a+ g (0)− 1]
(%o6) [a ex = a ex, a ex − cos (x) = a ex − cos (x)]
A solução de uma única EDO se dá de modo análogo:
(%i1) edo1: 'diff(y(t), t, 2) + 5*'diff(y(t), t) + 4*y(t) = t;
atvalue(y(t), t=0, 0);
atvalue('diff(y(t), t), t= 0, 0);
desolve(edo1,y(t));
(%o1)
d2
d t2
y (t) + 5
(
d
d t
y (t)
)
+ 4 y (t) = t
(%o2) 0
(%o3) 0
(%o4) y (t) =
e−t
3
− e
−4 t
48
+
t
4
− 5
16
Se desolve não pode obter uma solução, retorna false.
Embora o comando desolve seja muito útil, pois ele realiza a transformada de Laplace, resolve a(s)
equação(ões) no domínio de Laplace e realiza a transformada inversa, pode ser que em alguns casos
o usuário queira realizar cada uma destas etapas separadamente. Isto é o que ocorre no exemplo a
seguir, cuja solução final na forma gráfica está mostrada na Figura 5.6.1:
64
5. Solução analítica de equações diferenciais ordinárias (EDOs)
(%i1) eq1:'diff(y1(t),t,2)=-k*y1(t)+k*(y2(t)-y1(t));
eq2:'diff(y2(t),t,2)=-k*(y2(t)-y1(t))-k*y2(t);
atvalue(y1(t),t=0,1)$
atvalue(y2(t),t=0,1)$
atvalue('diff(y1(t),t),t=0,sqrt(3*k))$
atvalue('diff(y2(t),t),t=0,-sqrt(3*k))$
eq3:laplace(eq1,t,s);
eq4:laplace(eq2,t,s);
eq5:solve([eq3,eq4],[laplace(y1(t),t,s),laplace(y2(t),t,s)])[1];
assume(k>0);
eq6:expand(ilt(eq5[1],s,t));
eq7:expand(ilt(eq5[2],s,t));
plot2d([ev(rhs(eq6),k=1),ev(rhs(eq7),k=1)],[t,0,15],
[style, [lines,3,1], [lines,2,2]],[legend, "eq6", "eq7"],
[xlabel, "t"], [ylabel, "y1,y2"]);
(%o1)
d2
d t2
y1 (t) = k (y2 (t)− y1 (t))− k y1 (t)
(%o2)
d2
d t2
y2 (t) = −k (y2 (t)− y1 (t))− k y2 (t)
(%o7) s2 laplace (y1 (t) , t, s)− s−
√
3
√
k = k (laplace (y2 (t) , t, s)− laplace (y1 (t) , t, s))−
k laplace (y1 (t) , t, s)
(%o8) s2 laplace (y2 (t) , t, s)− s+
√
3
√
k = −k (laplace (y2 (t) , t, s)− laplace (y1 (t) , t, s))−
k laplace (y2 (t) , t, s)
(%o9) [laplace (y1 (t) , t, s) =
s3 +
√
3
√
k s2 + 3 k s+
√
3 k
3
2
s4 + 4 k s2 + 3 k2,
laplace (y2 (t) , t, s) =
s3 −√3√k s2 + 3 k s−√3 k 32
s4 + 4 k s2 + 3 k2
]
(%o10) [k > 0]
(%o11) y1 (t) = sin
(√
3
√
k t
)
+ cos
(√
k t
)
(%o12) y2 (t) = cos
(√
k t
)
− sin
(√
3
√
k t
)
5.7. Sistemas de EDOs lineares de 1ª ordem no plano
5.7.1. Autovalores e autovetores
Considere o seguinte sistema de EDOs lineares:
dy1
dx
= 11y1 − y2
dy2
dx
= y1 + 7y2
65
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.6.1.: Solução de um sistema de EDOs, por transformada de Laplace, sem usar o comando
desolve.
Que pode ser escrito na forma matricial como sendo:
dy
dx
=
(
11 −1
1 7
)
y
e: y = [y1 y2]
T
Os autovalores (ou valores próprios) e os autovetores (ou vetores próprios), associados ao sistema de
EDOs lineares, podem ser calculados usando as funções eigenvalues e eigenvectors.
(%i1) M1 : matrix ([11, -1], [1, 7]);
eigenvalues(M1);
eigenvectors(M1);
(%o1)
(
11 −1
1 7
)
(%o2) [[9−
√
3,
√
3 + 9], [1, 1]]
(%o3) [[[9−
√
3,
√
3 + 9], [1, 1]], [[[1,
√
3 + 2]], [[1, 2−
√
3]]]]
Perceber que a função eigenvalues retorna, na forma de uma lista, além dos autovalores (a primeira
sublista), a multiplicidade de cada um deles (a segunda sublista). Já o comando eigenvectors, além
de uma base de autovetores, retorna inicialmente os autovalores. Estas informações podem ser melhor
explicitadas na sequência de comandos a seguir:
(%i4) M1 : matrix ([11, -1], [1, 7]);
[vals, vecs] : eigenvectors (M1);
for i thru length (vals[1]) do disp (val[i] = vals[1][i],
mult[i] = vals[2][i], vec[i] = vecs[i]);
(%o4)
(
11 −1
1 7
)
(%o5) [[[9−
√
3,
√
3 + 9], [1, 1]], [[[1,
√
3 + 2]], [[1, 2−
√
3]]]]
val1 = 9−
√
3
mult1 = 1
66
5. Solução analítica de equações diferenciais ordinárias (EDOs)
vec1 = [[1,
√
3 + 2]]
val2 =
√
3 + 9
mult2 = 1
vec2 = [[1, 2−
√
3]]
(%o6) done
No exemplo seguinte, tem-se um sistema com multiplicidade maior do que um, a título de ilustração:
(%i7) M1 : matrix ([0, 1, 0, 0], [0, 0, 0, 0], [0, 0, 2, 0],
[0, 0, 0, 2]);
[vals, vecs] : eigenvectors (M1);
for i thru length (vals[1]) do disp (val[i] = vals[1][i],
mult[i] = vals[2][i], vec[i] = vecs[i]);
(%o7)

0 1 0 0
0 0 0 0
0 0 2 0
0 0 0 2

(%o8) [[[0, 2], [2, 2]], [[[1, 0, 0, 0]], [[0, 0, 1, 0], [0, 0, 0, 1]]]]
val1 = 0
mult1 = 2
vec1 = [[1, 0, 0, 0]]
val2 = 2
mult2 = 2
vec2 = [[0, 0, 1, 0], [0, 0, 0, 1]]
(%o9) done
5.7.2. Retrato de fase de sistemas lineares
Inicialmente será feito o plano de fases com o comando plotdf, já visto anteriormente. O retrato de
fase obtido está mostrado na Figura 5.7.1. Será considerado o seguinte sistema de EDOs:
dx1
dt
= −3x1 + x2
dx2
dt
= −2x1
x1(0) = 2; x2(0) = 0
(%i1) kill(all)$
/* entrada do sistema no estado estacionário */
eq1:-3*x1+x2$
eq2:-2*x1$
/* cálculo dos pontos de equilíbrio */
pontos_equilibrio:solve([eq1,eq2]);
/* extrai a matriz dos coeficientes */
A: coefmatrix([eq1,eq2],[x1,x2]);
/* calcula os autovalores e autovetores */
eigenvectors(A);
/* plota o retrato de fases */
vars: [x1, x2]$
plotdf([eq1,eq2],vars,[x1,-0.5,2],[x2,-1,0],[trajectory_at,0,-1])$
(%o3) [[x2 = 0, x1 = 0]]
67
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.7.1.: Retrato de fase e evolução no tempo. O gráfico da direita é obtido clicando no último
ícone da barra superior da figura à esquerda.
(%o4)
(−3 1
−2 0
)
(%o5) [[[−2,−1], [1, 1]], [[[1, 1]], [[1, 2]]]]
Outra maneira de obter o retrato de fase no Maxima, caso não haja sucesso com o procedimento
anterior, é efetuar o que seria feito manualmente. Neste caso o Maxima serve para automatizar os
cálculos. Este segundo procedimento permite entender melhor o processo de construir o retrato de
fase, sendo, portanto, mais didático.
Serão considerados dois exemplos: um sistema com nó (o mesmo considerado no exemplo anterior)
e outro com foco, ambos atratores. Inicialmente, resolve-se o sistema de EDOs por transformada
de Laplace com o comando desolve, com as condições iniciais pertinentes. Após serem calculados
os limites, para se verificar o comportamento do sistema no estado estacionário, o retrato de fase é
construído, usando um vetor da variável tempo e as soluções previamente encontradas. Os gráficos,
mostrando os retratos de fase, construídos a partir de pontos discretos, estão mostrados nas Figuras
5.7.2 e 5.7.3.
(%i1) kill(all)$
atvalue(x1(t),t=0,2)$
atvalue(x2(t),t=0,0)$
eq1:'diff(x1(t),t)=-3*x1(t)+x2(t);
eq2:'diff(x2(t),t)=-2*x1(t);
eq3:desolve([eq1,eq2],[x1(t),x2(t)]);
eq3a:eq3[1];
eq3b:eq3[2];
limit(rhs(eq3a),t,inf);
limit(rhs(eq3b),t,inf);
plot2d([rhs(eq3a),rhs(eq3b)],[t,0,20],[legend, "x1", "x2"],[ylabel, "x1,x2"])$
68
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.7.2.: Retrato de fase: nó atrator.
(%o3)
d
d t
x1 (t) = x2 (t)− 3 x1 (t)
(%o4)
d
d t
x2 (t) = −2 x1 (t)
(%o5) [x1 (t) = 4 e−2 t − 2 e−t, x2 (t) = 4 e−2 t − 4 e−t]
(%o6) x1 (t) = 4 e−2 t − 2 e−t
(%o7) x2 (t) = 4 e−2 t − 4 e−t
(%o8) 0
(%o9) 0
(%i11) tt:[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,1,1.25,1.5,1.75,2,
2.5,3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5]$
eq3a1:rhs(eq3a),t:tt,numer$
eq3a2:rhs(eq3b),t:tt,numer$
plot2d([discrete,eq3a1,eq3a2],[ylabel, "x2"],[xlabel, "x1"])$
69
5. Solução analítica de equações diferenciais ordinárias (EDOs)
(%i24) kill(all)$
A: matrix(
[-0.7,1],
[-1,0]);
lambda:eigenvalues(A);
eigenvectors(A);
expand(lambda),numer;
atvalue(x1(t),t=0,2)$
atvalue(x2(t),t=0,0)$
eq_a:A.matrix([x1(t)],[x2(t)]);
eq_a[1][1];
eq1:'diff(x1(t),t)=%;
eq_a[2][1];
eq2:'diff(x2(t),t)=%;
eq3:desolve([eq1,eq2],[x1(t),x2(t)]);
eq3a:eq3[1]$
eq3b:eq3[2]$
limit(rhs(eq3a),t,inf);
limit(rhs(eq3b),t,inf);
plot2d([rhs(eq3a),rhs(eq3b)],[t,0,20],[legend, "x1", "x2"],[ylabel, "x1,x2"])$
l1:create_list(i, i, 0,20)$
l2:create_list(i, i, 0.5,20.5)$
tt:join(l1,l2)$
eq3a1:rhs(eq3a),t:tt,numer$
eq3a2:rhs(eq3b),t:tt,numer$
plot2d([discrete,eq3a1,eq3a2],[ylabel, "x2"],[xlabel, "x1"])$
(%o1)
(−0.7 1
−1 0
)
(%o2) [[−3
√
39 i+ 7
20
,
3
√
39 i− 7
20
], [1, 1]]
(%o3) [[[−3
√
39 i+ 7
20
,
3
√
39 i− 7
20
], [1, 1]], [[[1,−3
√
39 i− 7
20
]], [[1,
3
√
39 i+ 7
20
]]]]
(%o4) [[−0.93674969975975 i− 0.35, 0.93674969975975 i− 0.35], [1, 1]]
(%o7)
(
x2 (t)− 0.7 x1 (t)
−x1 (t)
)
(%o8) x2 (t)− 0.7 x1 (t)
(%o9)
d
d t
x1 (t) = x2 (t)− 0.7 x1 (t)
(%o10) − x1 (t)
(%o11)
d
d t
x2 (t) = −x1 (t)
(%o12) [x1 (t) =
e−
7 t
20
(
20 cos
(
3
√
39 t
20
)
− 140 sin
(
3
√
39 t
20
)
3
√
39
)
10
, x2 (t) = −
40 e−
7 t
20 sin
(
3
√
39 t
20
)
3
√
39
]
(%o15) 0
(%o16) 0
70
5. Solução analítica de equações diferenciais ordinárias (EDOs)
Figura 5.7.3.: Retrato de fase: foco atrator.
Nos exemplos anteriores, ambos os sistemas eram lineares. Caso o sistema seja não linear, a análise é
feita de modo análogo, mas após a linearização do sistema. Não será abordado aqui este item, pois há
excelentes livros gratuitos sobre o assunto disponíveis para �download� na rede (por exemplo, Villate
(2007, 2013), disponível em: http://www.villate.org/docs.html).
71
6. Solução numérica de equações diferenciais
ordinárias usando o método de Runge-Kutta
A solução de EDOs numericamente, usando o Maxima, deve antes carregar o pacote dynamics que
inclui comandos para explorar sistemas dinâmicos e a implementação do método de Runge-Kutta de
4ª ordem.
O uso da função rk, para se resolver uma única EDO, é: rk(ODE, var, initial, domain)- EDO
a ser resolvida, variável dependente, valor inicial da variável dependente, lista contendo a variável
independente, valor inicial e final da variável independente, incremento a ser usado no intervalo anterior.
ODE deve ser uma expressão que depende somente das variáveis independente e dependentes e define
a derivada da variável dependente com relação à variável independente.
Veja o exemplo que se segue, juntamente com a solução, mostrada na Figura 6.0.1. Considere a solução
do Problema de Valor Inicial (PVI) dado por:
dy
dx
= 0, 2 exp(x) + xy
y(0) = 0; x ∈ [0, 4]
(%i9) kill(all)$
load(dynamics)$
edo:0.2*exp(x)+x*y;
solpvi:rk([edo],[y],[0],[x,0,4,0.05]) /* h=0.05, n=80 */$
plot2d([discrete,solpvi],[style,lines],[x,-0.1,4],[y,-100,2000]);
(%o2) x y + 0.2 ex
Para se resolver numericamente um sistema de EDOs, a função rk deve ser chamada, como segue: rk
([ODE1,...,ODEm], [v1,...,vm], [init1,...,initm], domain) - lista com as EDOs, lista com as
variáveis dependentes, lista com os valores iniciais das variáveis dependentes, lista contendo a variável
independente, valor inicial e final da variável independente, incremento a ser usado no intervalo anterior.
Deve haver apenas uma única variável independente. Veja o exemplo que se segue, juntamente com a
solução, mostrada na Figura 6.0.2.
dx
dt
= 4− x2 − 4y2
dy
dt
= y2 − x2 + 1
x(0) = −1, 25; y(0) = 0, 75; t ∈ [0, 4]
(%i7) sol: rk([4-x^2-4*y^2, y^2-x^2+1],[x,y],[-1.25,0.75],[t,0,4,0.02])$
plot2d ([[discrete,makelist([p[1],p[3]],p,sol)],
[discrete,makelist([q[1],q[2]],q,sol)]],
[xlabel,"t"],[ylabel,"x e y"],[legend,"y","x"])$
Outra função disponível no Maxima para solução numérica de EDOs é runge1, que pertence ao pacote
diffeq, cuja sintaxe é:
runge1(f, x0, x1, h, y0): f é a função f(x, y) da equação y' = f (x, y), x0 e x1 os valores iniciais x0,
e final, x1, da variável independente, respectivamente, h é o comprimento (ou passo) dos sub-intervalos
e y0 é o valor inicial y0 que toma y em x0.
72
6. Solução numérica de equações diferenciais ordinárias usando o método de Runge-Kutta
Figura 6.0.1.: Solução numérica de EDOs: única EDO.
Figura 6.0.2.: Solução numérica de EDOs: sistema de duas EDOs.
73
6. Solução numérica de equações diferenciais ordinárias usando o método de Runge-Kutta
O resultado é uma lista que por sua vez contém três listas: a primeira, contem a abcissa x; a segunda,
as ordenadas y; a terceira, as derivadas y'. O uso da função está demonstrado no exemplo que se segue,
para a mesma EDO resolvida no exemplo anterior:
(%i6) kill(all);
load(diffeq)$
f(x,y):=0.2*exp(x)+x*y$
h:1/20$
solnum:runge1(f,0,4,h,0)$
plot2d([discrete,solnum[1],solnum[2]])$
74
7. Solução analítica de equações diferencias
parciais (EDPs)
7.1. Exemplos de EDPs na Engenharia Química
Para recordar, uma EDP é formada de derivadas parciais, nas quais:
lim
∆x→0
f(x+ ∆x, t)− f(x, t)
∆x
=
∂f
∂x
, a t cte (7.1.1)
lim
∆t→0
f(x, t+ ∆t)− f(x, t)
∆t
=
∂f
∂t
, a x cte (7.1.2)
Se
dy
dx
= 0→ y(x) = cte (7.1.3)
Se
∂f
∂x
= 0→ f = g(t) (7.1.4)
As EDPs são classificadas do seguinte modo:
P
∂2z
∂x2
+Q
∂2z
∂x∂y
+R
∂2z
∂y2
= S (7.1.5)
ax2 + 2bxy + cy2 = d Portanto: ∆ = b2 − 4ac (7.1.6)
∆ < 0→ equação elíptica (7.1.7)
∆ = 0→ equação parabólica (7.1.8)
∆ > 0→ equação hiperbólica (7.1.9)
EXEMPLO 1: Reator tubular, oco, com paredes revestidas de catalisador, processando a reação A →
produtos (Rice and Do., 1995, item 10.4.1). O esquema do reator está mostrado na Figura 7.1.1.
75
7. Solução analítica de equações diferencias parciais (EDPs)
Figura 7.1.1.: Reator modelado no Exemplo 1.
DA
1
r
∂
∂r
(
r
∂CA
∂r
)
= v0
∂CA
∂z
(7.1.10)
com as condições de contorno : z = 0, CA(r, 0) = CA0 (7.1.11)
r = 0,
∂CA
∂r
= 0 ou CA é finito (7.1.12)
r = R, −DA
[
∂CA
∂r
]
= kCA (7.1.13)
(7.1.14)
A solução analítica do modelo anterior é (ver também Figura 7.1.2):
CA(r, z)
CA0
=
∞∑
n=1
2Ha
(λ2n +Ha
2)
J0(λnξ)
J0(λn)
exp(−λ2nζ) (7.1.15)
Ha = kR/DA (7.1.16)
ξ = r/R (7.1.17)
ζ =
z
v0
DA
R2
(7.1.18)
J1(λn) = Ha
J0(λn)
λn
(7.1.19)
Os valores de λn são calculados de modo a obedecer a última equação. No código que se segue,
calculam-se os 10 primeiros valores, para Ha = 10:
(%i1) eq1:bessel_j(1,%lambda)-Ha*bessel_j(0,%lambda)/%lambda;
load(newton1);
newton(eq1,%lambda,7,0.0001),Ha=10;
eq2:eq1,Ha=10;
for n:1 thru 10 step 1 do res[n]:newton(eq2,%lambda,n,0.0001);
for n:1 thru 10 step 1 do display(res[n])$
(%o1) bessel_j (1, λ)− bessel_j (0, λ) Ha
λ
(%o2) /usr/share/maxima/5.32.1/share/numeric/newton1.mac
(%o3) 7.956844234713721
(%o4) bessel_j (1, λ)− 10 bessel_j (0, λ)
λ
76
7. Solução analítica de equações diferencias parciais (EDPs)
Figura 7.1.2.: Solução, na forma gráfica, da EDP do Exemplo 1 (Rice and Do., 1995).
(%o5) done
res1 = 2.179496572558621
res2 = 2.17949656007887
res3 = 2.179495180034295
res4 = 5.033203763150477
res5 = 5.033211956888001
res6 = 13.95803019757838
res7 = 7.956844234713721
res8 = 7.956883405578219
res9 = 10.93633017715493
res10 = 10.93629509424527
Exemplo 2: Resfriamento de um fluido num tubo com escoamento laminar, ignorando condução axial
(Rice and Do., 1995, exemplo 10.3):
2v0
[
1−
( r
R
)2] ∂T
∂z
=
k
ρCp
1
r
∂
∂r
(
r
∂T
∂r
)
(7.1.20)
com as condições de contorno : z = 0, T = T0 (7.1.21)
r = 0,
∂T
∂r
= 0 (7.1.22)
r = R, T = Tw (7.1.23)
Cuja solução é:
T − Tw
T0 − Tw =
∞∑
n=1
AnGzn(ξ)exp(−λ2nζ) (7.1.24)
Gz(1, λn) = 0 ≡ funções de Graetz (7.1.25)
ξ =
r
R
, ζ =
zα
R22v0
, α =
k
ρCp
(7.1.26)
Exemplo 3: Condução de calor bidimensional no estado estacionário (eq. de Laplace):
77
7. Solução analítica de equações diferencias parciais (EDPs)
Figura 7.1.3.: Solução gráfica do Exemplo 3.
∂2T
∂x2
+
∂2T
∂y2
= 0 (7.1.27)
0 < x < a; T (x, 0) = 0; T (x, b) = 0 (7.1.28)
0 < y < b; T (0, y) = 0; T (a, y) = f(y) (7.1.29)
Exemplo 4: Determinação da difusividade térmica em corpos esféricos:
∂T
∂t
= α
1
r2
∂
∂r
(
r2
∂T
∂r
)
(7.1.30)
CI: t = 0, T (0, r) = Ts (7.1.31)
CC1: r = 0,
∂T
∂r
= 0 (7.1.32)
CC2: r = R, T (t, R) = Tf (7.1.33)
ξ =
r
R
, τ =
αt
R2
(7.1.34)
T (r, t) =
2(Tf − Ts)R
r
∞∑
n=1
(−1)n
npi
exp
[
−1(npi)
2αt
R2
]
sen
(npir
R
)
+ Tf (7.1.35)
Esfera de raio R = 1 cm, com temp. Inicial de Ts = 30
o
C, colocada num meio de temp. constante
de Tf = 20
o
C. Após 6,05 s, a temperatura no centro atingiu 20,86
o
C. Daí se calculou a difusividade
térmica do material como sendo 1/19 cm
2
/s.
78
7. Solução analítica de equações diferencias parciais (EDPs)
Figura 7.1.4.: Solução analítica na forma gráfica do Exemplo 4 (Rice and Do., 1995).
Problema para pensar (...e calcular): Quanto tempo uma almôndega congelada de soja, na qual se
usou ovo na sua fabricação, necessita ficar em óleo fervente para garantir que uma possível contaminação
com salmonela seja eliminada?
79
7. Solução analítica de equações diferencias parciais (EDPs)
Exemplo 5: Hemodiálise com microtubos poliméricos ocos, com a finalidade de retirar solutos indese-
jáveis (ex., ácido úrico)
Cada fibra oca pode ser modelada com um tubo reto com paredes permeáveis a certos solutos. Balanço
de massa para a espécie A:
v0
∂CA
∂z
= DA
1
r
∂
∂r
(
r
∂CA
∂r
)
(7.1.36)
com as condições de contorno : z = 0, CA = C0 (7.1.37)
r = 0,
∂CA
∂r
= 0 ou CA é finito (7.1.38)
r = R, −DA
[
∂CA
∂r
]
= K0L[CA − CD] (7.1.39)
1
K0L
=
tw
Dw
+
1
k0
(7.1.40)
CA
C0
= 4
∞∑
n=1
exp
[
−λ2nLDAv0R2
]
λ2n
[
1 + λ
2
n
Bi
2
]
(7.1.41)
Bi = K0LR/DA (7.1.42)
λnJ1(λn) = BiJ0(λn) (7.1.43)
Com a solução acima pode-se determinar a velocidade do sangue num único tubo para se diminuir a
concentração de ácido úrico de 20 mg/dL para 7,5 mg/dL.
7.2. Séries de Fourier
Em 1807 Fourier apresentou um trabalho à Academia de Paris sobre condução de calor no qual defendeu
a tese de que qualquer função arbitrária poderia ser representada através de um somatório envolvendo
as funções seno e cosseno. Embora o excesso de ambição por parte de Fourier, muitas funções de fato
podem ser representadas através de uma série convergente da forma:
80
7. Solução analítica de equações diferencias parciais (EDPs)
Figura 7.1.5.: Solução analítica na forma gráfica do Exemplo 5 (Rice and Do., 1995).
f(x) =
a0
2
+
∞∑
m=1
(
amcos
mpix
l
+ bmsen
mpix
l
)
, na qual : (7.2.1)
am =
1
l
∫ l
−l
f(x)cos
mpix
l
dx, m = 0, 1, 2, ... (7.2.2)
bm =
1
l
∫ l
−l
f(x)sen
mpix
l
dx, m = 1, 2, ... (7.2.3)
f(x) é uma função periódica e tem período 2l. Por exemplo, representar a função f a seguir por série
de Fourier:
f(x) =
{ −x , −l ≤ x < 0
x , 0 ≤ x < l
}
O cálculo dos coeficientes da série e a própria série de Fourier está mostrado a seguir:
m = 0
a0 =
1
l
∫ 0
−l(−x)dx+ 1l
∫ l
0(x)dx =
1
l
l2
2 +
1
l
l2
2 = l
m > 0
am =
1
l
∫ 0
−l(−x) cos mpixl dx+ 1l
∫ l
0(x) cos
mpix
l dx =
2l
(mpi)2
[cos(mpi)− 1], m = 1, 2, 3,...
= −4l/(mpi)2, se m for ímpar e , 0, se m for par
bm = 0,m = 1, 2, 3, ...
f(x) = l2 +
∑∞
m=1{−4l/(mpi)2} cos mpixl
= l2 − 4lpi2
(
cos pixl +
1
32
cos 3pixl +
1
52
cos 5pixl + . . .
)
= l2 − 4lpi2
∑∞
m=1,3,5...
cos(pimx/l)
m2
= l2 − 4lpi2
∑∞
n=1
cos{(2n−1)pix/l}
(2n−1)2
O código do Maxima está mostrado a seguir, juntamente com a Figura 7.2.1, que mostra a série
aberta para 1 a 4 termos no somatório (funções f4 a f7). Percebe-se claramente que a medida que
81
7. Solução analítica de equações diferencias parciais (EDPs)
se aumenta o número de termos da série, esta se aproxima mais da função original (ou função a
ser aproximada), mesmo que esta apresente descontinuidades (ou ausência de derivadas) nos pontos
x = (−1; −0, 5; 0; 0, 5; 1).
(%i1) kill(all);
f1:abs(x-1);
f2:abs(x);
f3:abs(x+1);
L:0.5 /* funcao de periodo 2*L = 2*0.5 = 1 */;
L/2 - 4*L/(%pi)^2*sum(cos((2*n-1)*%pi*x/L)/(2*n-1)^2,n,1,inf);
f4:L/2 - 4*L/(%pi)^2*sum(cos((2*n-1)*%pi*x/L)/(2*n-1)^2,n,1,1);
f5:L/2 - 4*L/(%pi)^2*sum(cos((2*n-1)*%pi*x/L)/(2*n-1)^2,n,1,2);
f6:L/2 - 4*L/(%pi)^2*sum(cos((2*n-1)*%pi*x/L)/(2*n-1)^2,n,1,3);
f7:L/2 - 4*L/(%pi)^2*sum(cos((2*n-1)*%pi*x/L)/(2*n-1)^2,n,1,4);
fig1:plot2d([f1,f2,f3,f4,f5,f6,f7],[x,-1,1],[y,0,0.5],
[legend,"f1","f2","f3","f4","f5","f6","f7"])$
(%o0) done
(%o1) |x− 1|
(%o2) |x|
(%o3) |x+ 1|
(%o4) 0.5
(%o5) 0.25−
2.0
∑∞
n=1
cos(2.0pi (2n−1)x)
(2n−1)2
pi2
(%o6) 0.25− 2.0 cos (2.0pi x)
pi2
(%o7) 0.25−
2.0
(
cos(6.0pi x)
9 + cos (2.0pi x)
)
pi2
(%o8) 0.25−
2.0
(
cos(10.0pi x)
25 +
cos(6.0pi x)
9 + cos (2.0pi x)
)
pi2
(%o9) 0.25−
2.0
(
cos(14.0pi x)
49 +
cos(10.0pi x)
25 +
cos(6.0pi x)
9 + cos (2.0pi x)
)
pi2
7.3. Solução de EDPs por separação de variáveis
O Maxima não dispõe de nenhuma função específica ou mesmo um pacote para a resolução analítica de
EDPs, dada a especificidade de cada caso. Entretanto, pode-se usar o Maxima nas etapas intermediárias
dos principais métodos, como será visto a seguir.
Será considerado o caso da condução de calor em uma barra com as extremidades a temperaturas
fixas. Supondo que a condução de calor se realiza numa geometria cilíndrica e que as variações de
temperatura nas direções radial e angular sejam desprezíveis, tem-se:
α
d2T
dz2
=
dT
dt
sendo 0 ≤ z ≤ L e t ≥ 0
Com as seguintes condições de contorno e inicial:
CI: t = 0, 0 ≤ z ≤ L, T (z, 0) = f(z)
CC1: t > 0, z = 0, T (0, t) = T1
82
7. Solução analítica de equações diferencias parciais (EDPs)
Figura 7.2.1.: Exemplo de série de Fourier no Maxima.
83
7. Solução analítica de equações diferencias parciais (EDPs)
CC2: t > 0, z = L, T (L, t) = T2
(%i1) depends(T,[z,t],zp,z,theta,t);
eq1:alpha*'diff(T,z,2)='diff(T,t);
(%o1) [T (z, t) , zp (z) , θ (t)]
(%o2) α
(
d2
d z2
T
)
=
d
d t
T
Definindo as variáveis separáveis, diferenciando, substituindo na EDP e separando as variáveis:
(%i3)
eq2:T=zp*theta;
eq3:diff(eq2,z,2);
eq4:diff(eq2,t);
eq5:ev(eq1,[eq3,eq4],diff);
eq6:eq5/(theta*zp*alpha);
(%o3) T = θ zp
(%o4)
d2
d z2
T = θ
(
d2
d z2
zp
)
(%o5)
d
d t
T =
(
d
d t
θ
)
zp
(%o6) α θ
(
d2
d z2
zp
)
=
(
d
d t
θ
)
zp
(%o7)
d2
d z2
zp
zp
=
d
d t θ
α θ
Gerar em seguir três conjuntos de EDOs e resolvê-las, um para cada caso de λ:
λ = 0
(%i8) eq7a:ode2(lhs(eq6),zp,z);
eq7b:ode2(rhs(eq6),theta,t);
(%o8) zp = %k2 z + %k1
(%o9) θ = %c
λ > 0 e checando o limite quando t tende ao infinito, a fim de verificar se a solução tem validade física
(%i10) assume(%lambda>0);
assume(alpha>0);
eq8a:ode2(lhs(eq6)+%lambda,zp,z);
eq8b:ode2(rhs(eq6)+%lambda,theta,t);
limit(eq8b,t,inf);
(%o10) [λ > 0]
(%o11) [α > 0]
(%o12) zp = %k1 sin
(√
λ z
)
+ %k2 cos
(√
λ z
)
(%o13) θ = %c e−λα t
(%o14) θ = 0
λ < 0, calculando o limite considerando a constante de integração sendo positiva ou negativa. Perceber
que neste caso a solução não é fisicamente viável, pois a temperatura no estado estacionário não é finita.
84
7. Solução analítica de equações diferencias parciais (EDPs)
(%i15) eq9a:ode2(lhs(eq6)-%lambda,zp,z);
eq9b:ode2(rhs(eq6)-%lambda,theta,t);
limit(eq9b,t,inf),%c=1;
limit(eq9b,t,inf),%c=-1;
(%o15) zp = %k1 e
√
λ z + %k2 e−
√
λ z
(%o16) θ = %c eλα t
(%o17) θ =∞
(%o18) θ = −∞
A seguir, usaremos valores numéricos para obter resposta na forma gráfica para uma barra de cobre de
1 m de comprimento, com extremidades mantidas a 25 C, distribuição inicial da temperatura T(z,0) =
75
o
C. O gráfico, mostrado na Figura 7.3.1, será obtido com a resposta com 100 termos no somatório:
(%i1) fz:75$
L:1$
%alpha:11.27*10^(-5)$
T1:25$
T2:25$
An1:2/L*integrate((fz-(T2-T1)/L-T1)*sin(n*%pi*z/L),z,0,L);
eq2:(T2-T1)*z/L + T1 +
sum(An1*sin(n*%pi*z/L)*%e^(-(n*%pi/L)^2*%alpha*t),n,1,100)$
eq2,z=0.5,t=10,numer;
eq2,z=0.5,t=30,numer;
eq2,z=0,t=0,numer;
eq2,z=1,t=0,numer;
plot3d(eq2,[z,0,1],[t,0,3600],[zlabel,"T"]);
(%o6) 100
(
1
pi n
− cos (pi n)
pi n
)
(%o8) 75.00000000000001
(%o9) 74.99999988011971
(%o10) 25.0
(%o11) 25.00000000000039
(%o12) /home/adilson/maxout.gnuplot_pipes
Como a temperatura varia no centro da barra em função do tempo? Ver Figura 7.3.2. Perceber que
no centro demoram-se alguns instantes para que a mudança nas extremidades seja sentida. Após, o
decaimento é exponencial, até a estabilização no novo estado estacionário.
(%i13) eq3:eq2,z=0.5$
plot2d([eq3],[t,0,3600],[xlabel,"t"],[ylabel,"T (z = 0,5)"])$
(%o14)
Como a temperatura varia com a posição em vários instantes de tempo? Ver Figura 7.3.3. Percebe-se
que tempos pequenos a maior parte da barra ainda está na temperatura inicial, embora as extremidades
da barra tenham suas temperaturas alteradas �instantaneamente�. Isto mostra a capacidade das séries
de Fourier em representar funções complexas, com grandes variações em pequenos intervalos. Para um
tempo muito grande, o novo estado estacionário é atingido.
85
7. Solução analítica de equações diferencias parciais (EDPs)
Figura 7.3.1.: Resposta, na forma tridimensional, para o exemplonumérico da condução de calor em
uma barra, com extremidades a temperaturas fixas.
Figura 7.3.2.: Variação da temperatura, com o tempo, no centro da barra.
86
7. Solução analítica de equações diferencias parciais (EDPs)
Figura 7.3.3.: Variação da temperatura, com a posição, para vários tempos diferentes.
(%i15) eq4:eq2,t=10$
eq5:eq2,t=50$
eq6:eq2,t=300$
eq7:eq2,t=700$
eq8:eq2,t=1500$
eq9:eq2,t=3600$
plot2d([eq4,eq5,eq6,eq7,eq8,eq9],[z,0,1],[ylabel,"T"],
[legend,"t=10","t=50","t=300","t=700","t=1500","t=3600"])$
87
8. Solução numérica de equações diferencias
parciais (EDPs) usando o método das
diferenças finitas
8.1. EDP parabólica
Todos os exemplos desta seção, juntamente com o código do Maxima, foram retirados de Sirin (2013).
Considere-se a equação da difusão transiente, unidimensional:
k
∂f
∂t
=
∂2f
dx2
na qual k é uma constante. Uma aproximação apropriada por diferenças finitas é:
k
f(i,j+1) − f(i,j)
∆t
=
f(i+1,j) − 2f(i,j) + f(i−1,j)
∆x2
na qual x = i∆x, i = 0, 1, 2,...,n, t = j∆t, j = 0,1,2,...,n. Na discretização foi utilizada a aproximação
à frente no lugar da derivada em relação a t e aproximação por diferença central para a derivada em
relação a x. Fazendo:
r =
∆t
k∆2
a equação discretizada por ser escrita como:
f(i,j+1) = r · f(i+1,j) + (1− r) · f(i,j) + r · f(i−1,j)
A solução está mostrada na Figura 8.1.1. As figuras 3D geradas pelo Maxima podem ser giradas em
qualquer ângulo e direção, bastando para isso clicar com o mouse da esquerda sobre ela e mantendo o
mouse pressionado, arrastá-la no ângulo ou direção desejado. Este recurso é muito útil quando o ângulo
padrão de exibição da figura não permite uma visualização conveniente dos resultados ou mesmo para
explorar outras perspectivas de visão em relação à solução obtida!
(%i17) declare([r,iend,jend,fini],constant)$
r:0.5$
iend:20.0$
jend:20$
fini:100.0$
for j:1.0 thru jend do for i:1.0 thru iend do f[i,j]:fini$
for j:1.0 thru jend do f[1.0,j]:0.0$
for i:iend thru iend do for j:1.0 thru jend do f[iend,j]:0.0$
f[1.0,1.0]:fini/2.0$
f[iend,1.0]:fini/2.0$
f[i,j]:for j:1.0 thru jend-1.0 do for i:2.0 thru iend-1.0 do
(f[i,j+1.0]:r*(f[i+1.0,j]+f[i-1.0,j])+(1-2*r)*f[i,j])$
load(draw)$
M:apply(matrix,makelist(makelist(f[i,j],j,1.0,19.0),i,1.0,19.0 ))$
wxdraw3d(contour_levels={10, 20,30,40,50,60,70,80,90,100},
contour = both,color = blue,elevation_grid(M,0,0,1,1),xlabel
= "x",ylabel = "y",surface_hide = true)$
im: apply(matrix,makelist(makelist (f[i,j],j,1.0,19.0),i,1.0,19.0))$
wxdraw2d(palette=gray,image(im,0,0,100,100))$
88
8. Solução numérica de equações diferencias parciais (EDPs) usando o método das diferenças finitas
Figura 8.1.1.: Solução numérica da equação da difusão, transiente, unidimensional.
8.2. EDP hiperbólica
A equação da onda é um exemplo clássico de EDP hiperbólica, dada por:
u2
∂2f
∂x2
=
∂2f
∂t2
na qual u é a velocidade da onda. Se forem aplicadas diferenças finitas centrais para as variáveis x e
t, obtem-se:
u2
f(i+1,j) − 2 · f(i,j) + f(i−1,j)
∆x2
=
f(i,j+1) − 2 · f(i,j) + f(i,j−1)
∆t2
Se for escolhido x = i∆x, t = j∆t, para (i,j) = 0,1,2,3,..., a equação anterior, na forma discretizada,
pode ser escrita como:
f(i,j+1) = 2(1− r) · f(i,j) + r
{
f(i+1,j) + f(i−1,j)
}− f(i,j−1)
na qual f(i,j)é uma aproximação para f(x,t) e r é a �razão de aspecto� dada por:
r =
(
u
∆t
∆x
)2
Para as condições de contorno; t ≥ 0 e 0 < x < 1: f(0, t) = 0 e f(1, t) = 0
Com as condições iniciais:
f(x, 0) = senpix e
(
∂f
∂t
)
= ft(x, 0) = 0
Considerando as condições iniciais e tomando arbitrariamente o valor de r como sendo 1 (r = 1):
f(i,j+1) = f(i−1,j) + f(i+1,j) − f(i,j−1) j ≥ 1
ft =
f(i,1) − f(i,−1)
2 ·∆t j = 0
Obtem-se:
f(i,1) = f(i,−1) e f(i,1) =
1
2
{
f(i−1,1) + f(i,1)
}
89
8. Solução numérica de equações diferencias parciais (EDPs) usando o método das diferenças finitas
Figura 8.2.1.: Solução numérica da equação da onda, transiente, unidimensional.
A solução no Maxima da equação da onda, juntamente com as condições iniciais e de contorno consi-
deradas, está a seguir. Na Figura 8.2.1, mostra-se o resultado na forma gráfica.
(%i16) kill(all)$
declare([iend,jend,dx,dt],constant)$
iend:20$
jend:20$
dx:.15$
dt:.5$
for i from 1 thru 20 do for j from 1 thru 20 do (f[i,j]:0)$
for i:1 thru iend do (xx[i]:dx*i)$
for i:1 thru iend do (f[i,1]:sin(3.14*xx[i]))$
for i:2 thru iend-1 do (f[i,2]:0.5*(f[i-1,1]+f[i+1,1]))$
f[i,j]:for j:2 thru jend-1 do for i:2 thru iend-1 do
(f[i,j+1]:(f[i+1,j]+f[i-1,j]-f[i,j-1]))$
M:apply(matrix,makelist(makelist(f[i,j+1],j,2,19),i,2,19))$
for i from 1 thru iend do (y[i]:dx*i)$
for j from 1 thru jend do (x[j]:dt*j)$
load(draw)$
draw3d (contour_levels=
{x[1],x[2],x[3],
x[4],x[5],x[6],x[7],x[8],x[9],x[10],x[11],
x[12], x[13],x[14],x[15], x[16], x[17],x[18],x[19],x[20]},
contour = both, color = blue, elevation_grid(M,0,0,8,24),
xlabel = "Time", ylabel = "x",zlabel = "Function
Value",surface_hide =false)$
90
8. Solução numérica de equações diferencias parciais (EDPs) usando o método das diferenças finitas
8.3. EDP elíptica
As equações de Poisson:
∇2f = ∂
2f
∂x2
+
∂2f
∂y2
= g(x, y)
e de Laplace:
∇2f = ∂
2f
∂x2
+
∂2f
∂y2
= 0
são típicas EDPs elípticas. Quando se usa aproximação por diferenças centrais para as variáveis x e y,
obtem-se:
∂2f
∂x2
≈ f(i,j+1) − 2 · f(i,j) + f(i,j−1)
(∆x)2
∂2f
∂y2
≈ f(i,j+1) − 2 · f(i,j) + f(i,j−1)
(∆y)2
usando x = i∆x, i = 0,1,2,...,n e y = j∆y, j = 0,1,2,...,m. Ao se escolher uma malha quadrada de
discretização, ou ∆x = ∆y = h, obtem-se:
f(i,j) =
1
4
{f(i+1,j) + f(i,j) + f(i,j+1) + f(i,j−1)} −
h2
4
g(i, j)
A solução da equação de Laplace está mostrada a seguir, com os resultados na Figura 8.3.1.
(%i52) /* solução numérica da equação de Laplace */
nx:10$ ny:10$ a:1.0$ b:1.0$ v1:10.0$ v2:60.0$ v3:0.0$ v4:30.0$
for i:1 thru nx-1 do for j:1 thru ny-1
do(v[i,j]:(v1+v2+v3+v4)/4)$
fpprintprec:4$
for i:1 thru nx-1 do (v[i,1]:v1)$
for i:1 thru nx-1 do(v[i,ny]:v3)$
for j:1 thru ny-1 do (v[1,j]:v4)$
for j:1 thru ny-1 do(v[nx,j]:v2)$
for k:1 thru 23 do for i:2 thru nx-1 do for j:2 thru ny-1 do
(v[i,j]:0.25*(v[i+1,j]+v[i-1,j]+v[i,j+1]+v[i,j-1]))$
load(draw)$
M:apply(matrix,makelist(makelist(v[i,j],j,1,ny-1),i,1,nx))$
draw3d(contour_levels = {10,15,20, 25,30,35, 40,45,50,60,
70,80,90,100},contour = both,color = blue,
elevation_grid(M,0,0,100,100),xlabel = "x",
ylabel = "y",surface_hide = true)$
91
8. Solução numérica de equações diferencias parciais (EDPs) usando o método das diferenças finitas
Figura 8.3.1.: Solução numérica da equação de Laplace.
(%i70) /* solução numérica da equação de Poisson */
nx:8$ ny:8$ a:1.0$ h:a/nx$ v1:10.0$ v2:-30.0$ v3:100.0$ v4:-10.0$
for i:1 thru nx-1 do for j:1 thru ny-1 do
(v[i,j]:(v1+v2+v3+v4)/4)$
for i:1 thru nx-1 do (v[i,1]:v1,v[i,ny]:v3)$
for j:1 thru ny-1 do (v[1,j]:v4,v[nx,j]:v2)$
t:cos(3.14/nx)+cos(3.14/ny)$
w: (8.0-sqrt(64.0-16.0*t^2))/(t^2)$
w4:w/4.0$
for k:1 thru 10 do for i:2 thru nx-1 do for j:2
thru ny-1 do (v[i,j]:(x:h*i,y:h*j,g:-36.0*3.14*x*(y-1.0),
r:w4*(v[i+1,j]+v[i-1,j]+v[i,j+1]+v[i,j-1]-4*v[i,j]-
g*h^2),v[i,j]:v[i,j]+(w4)*r))$
load(draw)$
M:apply(matrix,makelist(makelist (v[i,j],j,1,7),i,1,7))$
draw3d(contour_levels = {10,20,30, 40,50, 60,
70,80,90, 100}, contour= both,color = blue, elevation_grid
(M,0,0,1,1), xlabel = "x",ylabel = "y",surface_hide = true)$
92
8. Solução numérica de equações diferencias parciais (EDPs) usando o método das diferenças finitas
Figura 8.3.2.: Solução numérica da equação de Poisson.
93
Referências Bibliográficas
William E. Boyceand Richard C. DiPrima. Equações diferenciais elementares e problemas de valores
no contorno. LTC, 9 edition, 2010.
Fabio Luiz de Oliveira. Abordagem de conceitos de funções de duas variáveis com o uso do software
maxima. Produto Educacional do Mestrado Profissional em Educação Matemática da Universidade
Federal de Ouro Preto, UFOP, 2014. URL http://www.ppgedmat.ufop.br/arquivos/produtos_
2014/Produto_Fabio%20Luiz%20de%20Oliveira.pdf.
Naimara Vieira do Prado, Pétterson Vinícius Pramiu, and Rogério Luis Rizzi. O Emprego do software
Maxima no apoio ao ensino de Matemática. UNIOESTE, Cascavel, 2008. URL ftp://ftp.feq.
ufu.br/adilson/MM/apostilas_Maxima/MaximaApoioEnsinoMatematica.pdf.
Erwin Kreyszig, Herbert Kreyszig, and Edward J. Norminton. Advanced engineering mathematics.
John Wiley, 10 edition, 2011.
Renato A Nodarse. Introduccíon al MAXIMA CAS con aplicaciones a la resolucíon de ecuaciones
diferenciales ordinarias. "on line", 2014. URL http://euler.us.es/~renato.
R. G. Rice and D. D. Do. Applied mathematics and modeling for chemical engineers. John Wiley,
1995.
Bruna Santos. Introdução ao software Maxima. Centro de Matemática da Universidade do Porto,
2009. URL http://www.tellau.com.br/ufes/2014-2/LabMat/Aulas/Aula6/maxima/apostilas/
document.pdf.
Tufan Sirin. Solutions of partial differential equations with computer algebra. In 2013 International
Conference on Electronics, Computer and Computation (ICECCO), pages 351�354. IEEE, 2013.
doi: 10.1109/ICECCO.2013.6718300. URL http://ieeexplore.ieee.org/xpl/articleDetails.
jsp?tp=&arnumber=6718300&url=http%3A%2F%2Fieeexplore.ieee.org%2Fxpls%2Fabs_all.jsp%
3Farnumber%3D6718300.
Cristina M. Torres, Filipe R. Polizel, Letícia A. B. Silva, Luís Fernando L. Sato, Natália do Carmo Car-
valho, Thiago M. Sfredo, and Tiago B. de Lima. Utilização do software Maxima. UNICAMP, 2006.
URL http://www.ime.unicamp.br/~marcio/ss2006/grupo10.pdf.
Jaime E. Villate. Introdução aos Sistemas Dinâmicos: uma abordagem prática com Maxima. author's
edition, Porto, 2007. URL http://www.villate.org/doc/sistemasdinamicos/sistdinam-1_2.
pdf.
Jaime E. Villate. Equações diferenciais e equações de diferenças. "on line", 2011. URL http://www.
villate.org/doc/eqdiferenciais/eqdif_20110426.pdf.
Jaime E. Villate. Dinâmica e Sistemas Dinâmicos. author's edition, Porto, 2013. URL http://def.
fe.up.pt/pt/Din%C3%A2mica_e_Sistemas_Din%C3%A2micos.
Wikipédia. Modelagem matemática � wikipédia, a enciclopédia livre, 2013. URL http://pt.
wikipedia.org/w/index.php?title=Modelagem_matem%C3%A1tica&oldid=37475222. [Online; ac-
cessado em 15 de fevereiro de 2015].
Edwin L. Woollett. Maxima by example, 2009. URL http://www.csulb.edu/~woollett/.
94
A. Lista de exercícios n. 1: Revisão de
matemática básica e cálculo diferencial e
integral
�
�
�
�
Universidade Federal de Uberlândia � Faculdade de Engenharia Química
Programa de Pós-Graduação em Engenharia Química � Métodos Matemáticos em Engenharia Química
Prof. Adilson J. de Assis
1. (a) Converta 45º em radianos;
(b) Converta pi/6 radianos em graus;
(c) Avalie as expressões sem uma calculadora: (−3)4, −34, 3−4, 5
23
521
,
(
2
3
)−2
, 16−3/4;
(d) Simplifique a expressão:
(
3x3/2y3
x2y−1/2
)−2
;
(e) Resolva as equações:
(i)
2x
x+ 1
=
2x− 1
x
;
(ii) x4 − 3x2 + 2 = 0;
2. Encontre os limites:
(a) lim
x→1
ex
3−x
;
(b) lim
x→3
x2 − 9
x2 + 2x− 3 ;
(c) lim
x→−3
x2 − 9
x2 + 2x− 3 ;
(d) lim
x→1+
x2 − 9
x2 + 2x− 3 ;
(e) lim
h→0
(h− 1)3 + 1
h
;
(f) lim
x→pi−
ln(sen(x));
(g) lim
x→−∞
1− 2x2 − x4
5 + x− 3x4 ;
(h) lim
x→∞ e
x−x2
;
3. Derive as funções:
(a) f(t) = sen(3t+ 1);
(b) f(t) = cos2t;
(c) f(t) = tan(t2) + sec(1− 3t);
(d) f(t) = e−x sen x;
(e) f(u) = e−u/2 cos(2piu);
95
A. Lista de exercícios n. 1: Revisão de matemática básica e cálculo diferencial e integral
(f) f(u) =
cos(u)
1− cos(u) ;
(g) f(x) = ln
x+ 1√
x− 2 ;
4. Usando séries de potências, mostre que:
(a) ex =
∞∑
n=0
xn
n!
, para todo x; qual foi o ponto de expansão considerado neste caso? Usando a
série de potências, calcule e2,5 com 5 termos no somatório.
(b) Considerando ainda ex e o ponto de expansão a = 2, compare e2,5 com o resultado anterior,
também com 5 termos no somatório; interprete os resultados.
(c) Mostre que o seno de um ângulo pode ser calculado pela soma:
sen(x) = x− x
3
3!
+
x5
5!
− x
7
7!
+ . . .
=
∞∑
n=0
(−1)n x
2n+1
(2n+ 1)!
para todo x;
5. Integre as funções:
(a)
∫
sec t(sec t+ tan t)dt;
(b)
∫
x sen(x2)dx;
(c)
∫
sen(x) cos(x)dx;
(d)
∫ pi/6
0
sen2(t) cos(t)dt;
(e)
∫
5x− 4
2x2 + x− 1dx (OBS: use a técnica das frações parciais);
(e)
∫ 2
0
x2 + 12
x2 + 4
dx (OBS: use uma tabela de integrais);
6. Determine a área da região sob a curva y = x sen(x) entre as retas x = 0 e x = pi.
7. Suponha que o número de roupas de banho vendidas em uma loja de Salvador durante a semana
t de um certo ano pode ser modelada pela função:
B(t) = 40 + 36 cos
[ pi
52
(t− 4)
]
para
−20 6 t 6 32
(a) Em que semana as vendas são maiores?
(b) Com que taxa as vendas estão variando no fim de junho (para t = 26)?
(c) A experiência mostra que as vendas tendem a aumentar entre as semanas 0 e 8. Calcule a
média semanal de vendas durante esse período.
8. Se f(x, y) = x3 + x2y3 − 2y2, encontre fx(2, 1) e fy(2, 1) e interprete estes números como incli-
nações.
9. Encontre
∂z
∂x e
∂z
∂y se z é definido como uma função de x e y pela equação:
x3 + y3 + z3 + 6xyz = 1
96
A. Lista de exercícios n. 1: Revisão de matemática básica e cálculo diferencial e integral
10. Se f(x, y) = x3 + x2y3 − 2y2, encontre as derivadas segundas de f (ou seja: fxx, fxy, fyx, fyy).
11. (a) A lei dos gases para uma massa fixa m de um gás ideal na temperatura absoluta T , pressão
P e volume V é PV = mRT , onde R é a constante dos gases. Mostre que
∂P
∂V
∂V
∂T
∂T
∂P
= −1
e interprete fisicamente os termos da expressão acima e como ela pode ser utilizada na
prática para verificar a idealidade de um determinado gás.
(b) Mostre que, para um gás ideal
T
∂P
∂T
∂V
∂T
= mR
12. O índice de resfriamento do vento é modelado pela função
W = 13, 12 + 0, 6215T − 11, 37v0,16+
0, 3965Tv0,16
onde T é a temperatura (ºC) e v é a velocidade do vento (km/h). Quando T = -15 ºC e v = 30
km/h, quanto você espera que a temperatura aparente W irá cair se a temperatura verdadeira
cair de 1 ºC? E se a velocidade do vento aumentar de 1 km/h?
97
B. Lista de exercícios n. 2: Classificação,
campo de direções e EDOs de 1ª ordem
�
�
�
�
Universidade Federal de Uberlândia � Faculdade de Engenharia Química
Programa de Pós-Graduação em Engenharia Química � Métodos Matemáticos em Engenharia Química
Prof. Adilson J. de Assis
1. O valor de lim
x→∞
(1, 1)x
x1000
é:
(a) 0
(b) ∞
(c) 1
(d) 1, 1
(e) n.d.a.
2. Se f(x) =
√
x cos(1− x2), então df
dx
é:
(a) 2x3/2 sen
(
1− x2)+ 1
2
cos
(
1− x2)x−1/2
(b)
cos
(
x2 − 1)
2
√
x
− 2x 32
(c) 0
(d)
cos
(
x2 − 1)
2
− 2x 32 sin (x2 − 1) √x
(e) n.d.a.
3. O resultado de
∫
1
ex − 1 dx é:
(a)
−1
ex − 1
(b) log (ex − 1)− x
(c) x log (e)− log 1− x
(d) ln (ex − 1)− x
(e) n.d.a.
4. Classifique as seguintes EDs:
(a)
dy
dx + 2xy = sen(x)
(b) (1− x2) d2y
dx2
− 2x dydx + 6y = 0, com: − 1 < x < 1
(c)
d2y
dx2
+ a dydx + by = sen(ωx), com: ω = cte
(d) (x2y + x)dy + (xy2 − y)dx = 0
5. Mostre que y = 13x
3 + 32x
2 + C satisfaz a seguinte EDO: dydx = x
2 + 3x
98
B. Lista de exercícios n. 2: Classificação, campo de direções e EDOs de 1ª ordem
6. (a)Determine f de modo que y = Cx exp(−x) + f seja a solução geral de
x
dy
dx
+ (x− 1)y − x2 = 0
(a) Se y = 1 em x = 1, qual a solução particular desta EDO?
7. Considere o campo de direções obtido no Maxima com o comando:
plotdf(exp(-x)+y,[trajectory_at,0,0])
Para cada sentença abaixo, marque V se for Verdadeira e F se for Falsa.
( ) A curva (em vermelho) obtida é uma solução singular da EDO cujo campo de direções é
mostrado na figura;
( ) A curva (em vermelho) obtida é a solução do PVI composto da EDO, cujo campo de direções
é mostrado na figura, e a condição inicial xi = 0, yi = 0;
( ) Para qualquer condição inicial (xi, yi) > 0, tem-se que y > 0 se a constante de integração
= 1;
( ) A EDO correspondente ao campo das direções é:
dy
dx − x exp(−x) + yx
( ) A técnica do campo das direções é útil para se conhecer o comportamento das soluções de
uma EDO de 1ª ordem, mesmo antes de resolvê-la;
8. Suponha que em uma certa reação química autocatalítica o composto A reage para formar o com-
posto B. Além disto, suponha que a concentração inicial de A é CAo e que CB(t) é a concentração
de B no tempo t. Então, CAo − CB(t) é a concentração de A no tempo t. Como numa reação
autocatalítica a substância produzida estimula a reação, então a taxa da reação
CB
dt é proporci-
onal tanto a CB(t), quanto a CA(t) = CAo − CB(t). Portanto, dCB(t)
dt
= kCB(t)(CA0 − CB(t)).
Suponha que no início, t = 0, a concentração de B seja CBo.
(a) Classifique (quanto à ordem, à linearidade e ao tipo de coeficiente etc) a equação dada que
descreve a variação da concentração de B com o tempo;
(b) Resolva manualmente a EDO que descreve a variação da concentração de B com o tempo;
Resposta: CB(t) = 1/(C
−1
Ao + Ce
−kCAot), com:C = (CAo − CBo)/(CAoCBo)
9. A taxa com a qual uma epidemia se espalha em uma comunidade é conjuntamente proporcional ao
número de pessoas infectadas e ao número de pessoas suscetíveis que ainda não foram infectadas.
A seguinte EDO expressa o número de pessoas infectadas em função do tempo:
dQ
dt
= kQ(B−Q),
na qual k é uma constante de proporcionalidade (k > 0) e B é o número total de pessoas
suscetíveis.
(a) Mostre que a solução a seguir satisfaz a EDO: Q =
B
1 +Ae−Bkt
(A é uma constante);
(b) O que acontece com o número de pessoas infectadas no início da epidemia e para um
tempo muito longo, de acordo com este modelo? Resposta: Início: Q = B/(1+A);
tempo muito longo: Q = B;
10. Considere o campo (ou mapa) de direções mostrado na Figura B.0.1. Este é o campo de direções
de qual EDO? Justifique apropriadamente sua escolha.
(a) dy/dx− x = 0
(b) dy/dx− cos(x+ y) = 0
(c) dy/dx− y2 = x
(d) dy/dx− ex = 1
(e) dy/dx = cos(x)
99
B. Lista de exercícios n. 2: Classificação, campo de direções e EDOs de 1ª ordem
Figura B.0.1.: Campo de direções de uma determinada equação diferencial.
11. [solução manual]
1
Dois líquidos A e B estão em ebulição juntos num mesmo recipiente. Expe-
rimentalmente verifica-se que a razão das velocidades segundo as quais A e B estão evaporando
em qualquer instante de tempo é proporcional à razão da quantia de A (digamos, x) em relação
à quantia de B (digamos, y) ainda no estado líquido. Esta lei física é matematicamente repre-
sentada por (dy/dt)/(dx/dt) = k y/x, onde k é uma constante de proporcionalidade. Encontre
y em função de x, ou seja, como a quantia remanescente de B se relaciona com a quantia de A.
Resposta: y(x) = Cxk = Cek ln(x)
12. Resolva a EDO a seguir: (3x2 + 2y + 2 cosh(2x + 3y))dx + (2x + 2y + 3 cosh(2x + 3y))dy = 0
Resposta: x3 + 2xy + y2 + senh(2x+ 3y) = C
13. [solução manual] Considere a seguinte EDO: (x2− 1)dy
dx
= A− 2x y, onde A pode ser função
só de x [ou A(x)], função só de y [ou A(y)] ou função de ambas as variáveis [ou A(x, y)].
(a) Analise a forma que A deve necessariamente assumir para que a EDO seja exata.
(b) Se A = cos(x), qual função y satisfaz a EDO? Resposta: y = sen(x)+C
x2−1
14. Considere um tanque cilíndrico de área transversal constante, como o mostrado na Figura
B.0.2.Realizando um balanço de massa no sistema, tem-se:
(entra) - (sai) + (geração) = (acúmulo)
(entra) = q
(sai) = qo = Ao
√
2gh(t)
(geração) = 0
(acúmulo) =
dV
dt
então:
q −Ao
√
2gh(t) =
dV
dt
1
Significa que o Maxima poderá ser utilizado para auxiliar na solução, mas o desenvolvimento deverá ser manual.
100
B. Lista de exercícios n. 2: Classificação, campo de direções e EDOs de 1ª ordem
Como a área da seção transversal é constante, V = A.h(t) e, portanto:
A
dh(t)
dt
= q −Ao
√
2gh(t)
Figura B.0.2.: Tanque cilíndrico com entrada e saída de líquido.
(a) Classifique a EDO anterior e verifique se ela pode ser resolvida por separação de variáveis;
(b) Calcule a altura do líquido no estado estacionário, ou seja, no ponto de equilíbrio do sistema;
(c) Faça x =
√
h e substituindo na EDO do modelo, encontre uma EDO separável, resolvendo-a
de modo apropriado, encontrando uma expressão do tempo em função da altura de líquido;
Resposta: t =
(
qA
A2og
)
ln
(
q
q −Ao
√
2gh
)
− 2
(
A
Ao
)√
h
2g
(d) Quando não há vazamento, qo = 0, ou seja, Ao = 0. Nesta situação, encontre como o tempo
varia com a altura;
(e) Esboce num mesmo gráfico h× t o comportamento qualitativo para os dois casos anteriores,
destacando as semelhanças e diferenças entre eles;
15. [solução manual] Por volta das 22 h em uma noite animada de sexta-feira uma casa noturna
de dimensões 30 × 50 × 10 ft está lotada de pessoas se divertindo. Infelizmente, muitas destas
pessoas são fumantes, de tal forma que fumaça dos cigarros contendo 4% de monóxido de carbono
(CO) é introduzida no ambiente numa vazão de 0,15 ft
3
/min. Suponha que esta vazão não
varie significantemente durante a noite. Antes das 22 h não há traços de CO no ambiente;
e, felizmente, a casa noturna está equipada com bons ventiladores que permitem a formação
uniforme da mistura ar-fumaça no ambiente. Além disto, há um sistema de exaustão que retira
1,5 ft
3
/min da mistura contida no ambiente; isto equivale a uma vazão 10 vezes maior do que a
entrada de poluentes.
Você deseja dançar e se socializar, mas também quer preservar sua saúde. Uma exposição pro-
longada a concentrações de CO maiores ou iguais a 0,012% é considerada perigosa para a saúde,
de acordo com a OMS. Sabendo que o clube fecha as 4 h, é prudente ficar no ambiente até o fe-
chamento ou é conveniente sair um pouco mais cedo? Para ser mais preciso, você quer encontrar
o momento quando a concentração de CO atinge a concentração crítica de 0,012%.
Sendo C(t) a concentração de CO no ambiente (gramas de CO por pé cúbico de ar, ou g/ft3) em
qualquer tempo t, onde t = 0 representa 22 h, então Q(t), a quantia de poluente no ambiente
no tempo t, é descrita pela equação Q(t) = (volume da sala) ×C(t). Como a sala tem como
101
B. Lista de exercícios n. 2: Classificação, campo de direções e EDOs de 1ª ordem
dimensões 30 × 50 × 10 ft = 15 000 ft3, então: Q(t) = 15 000C(t). A taxa segundo a qual o CO
está entrando na sala é dada por:(
0, 15
ft
3
min
)(
0, 04
g
ft
3
)
= 0, 006
g
min
Similarmente, a taxa na qual o CO está saindo da sala (via sistema de ventilação) é
(
1, 5
ft
3
min
)
C(t).
Um balanço de massa no ambiente, para o CO, nos diz que o acúmulo de CO na sala é igual à
entrada de CO, menos a saída de CO, ou seja:
dQ(t)
dt
=
(
0, 15
ft
3
min
)(
0, 04
g
ft
3
)
−
(
1, 5
ft
3
min
)
C(t)
(a) Qual é a equação diferencial então que permite encontrar a concentração de CO em função
do tempo?
(b) Classifique esta equação diferencial quanto: à ordem, à homogeneidade, ao tipo de coeficiente
e quanto a linearidade;
(c) Trata-se ou não se um PVI? Se sim, explicite a equação diferenciale a(s) condição(ões)
inicial(ais);
(d) Existe um fator de integração para esta equação diferencial? Se sim, qual é? Resposta:
Sim, F = e0,0001t
(e) Qual é a solução geral do problema? Resposta: C(t) = 4× 10−3 + 4× 10−7ke−10−4t
(f) E a solução particular? Resposta: equação anterior, com k = −1× 10−4
(g) Após quanto tempo a concentração de CO atingirá o valor-limite de 0,012%? Então, qual
seria a hora-limite para ficar no ambiente sem prezuízo à saúde, de acordo com a OMS?
Resposta: após 5h5min depois das 22 h (22 h = início do cálculo ou t=0) ou às 3h05 da
madrugada
(h) Esquematize o mais detalhadamente possível C(t)× t até t →∞
102
C. Lista de exercícios n. 3: EDOs de 2ª ordem,
solução em série de potências de EDOs
�
�
�
�
Universidade Federal de Uberlândia � Faculdade de Engenharia Química
Programa de Pós-Graduação em Engenharia Química � Métodos Matemáticos em Engenharia Química
Prof. Adilson J. de Assis
1. Considere a EDO que descreve a perda de calor em flanges de tubos: r
d2T
dr2
+
dT
dr
− h
6k
r(T−T1) = 0
(a) Mostre que a substituição de variáveis y = T −T1 e x = r
√
h/6k leva à equação modificada
de Bessel x2
d2y
dx2
+ x
dy
dx
− x2y = 0
(b) Tente resolver a EDO do item (a) pelo Maxima com o comando ode2. O que acontece?
Agora tente com os comandos a seguir e interprete o resultado:
eq1:x^2*'diff(y,x,2)+x*'diff(y,x)-x^2*y=0;
load(contrib_ode);
sol2: contrib_ode(eq1,y,x);
(c) Compare a solução obtida com o Maxima com a solução padrão da equação modificada de
Bessel: y(x) = c1I0(x) + c2K0(x)
(d) Se um tubo de 1 in de diâmetro (externo) está conectado a flanges de 4 in de diâmetro e
1/2 in de espessura e os tubos carregam vapor a 250 ºC e sendo a condutividade térmica do
metal do flange igual a k = 220 Btu/h.ft2.ºF.ft−1, a vizinhança com T1 = 60ºC, o coeficiente
de transferência de calor h = 2 Btu/h.ft2.ºF, as condições de contorno apropriadas são: i)
na superfície do tubo: r = 1/2, T = 250; ii) na extremidade do flange, o calor que chega
por condução (
−pikr
12
dT
dr ) é igual ao calor que é dissipado por convecção (
2pirh∆r
144 (T −T1)) [no
estado estacionário] ou r = 2, −k12
dT
dr =
h
144(T −T1). Mostre que estas condições de contorno
são equivalentes a: i) x=0,0195; y = 190; ii) x = 0,078; dy/dx = -0,0195y.
(e) Usando as condições de contorno apropriadas e a solução padrão da equação modificada
de Bessel, mostre que T = 60 + 186, 5I0(0, 039r) + 0, 8636K0(0, 039r) e também calcule a
temperatura na extremidade do flange.
2. Em controle de processos, um sistema de 2ª ordem é aquele descrito por uma EDO de 2ª ordem.
Se o sistema for linear:
a2
d2y
dt2
+ a1
dy
dt
+ a0y = b.f(t)
(a) Divida todos os termos por a0 a fim de se obter:
τ2
d2y
dt2
+ 2ξτ
dy
dt
+ y = Kp.f(t)
onde: τ = a2a0 , 2ξτ =
a1
a0
,Kp =
b
a0
.
(b) Na forma de variável desvio, as condições iniciais são: y(0) = 0 e y′(0) = 0. Analise o tipo
de resposta qualitativa que será encontrada para um sistema de 2ª ordem para os casos: a)
ξ > 1 [chamado de sistema superamortecido]; b) ξ = 1 [chamado de sistema criticamente
amortecido]; ξ < 1 [chamado de sistema subamortecido]; esboce num mesmo gráfico y × t
os 3 comportamentos;
1
1
Curiosidade: Em controle de processos as raízes da equação característica são chamadas de pólos do sistema!
103
C. Lista de exercícios n. 3: EDOs de 2ª ordem, solução em série de potências de EDOs
3. No projeto e análise de reatores químicos é importante encontrar a relação para predizer o perfil
de composição em um reator tubular com recheio de catalisador empacotado processando uma
reação isotérmica, com cinética linear e com difusão axial. Este tipo de reator, além de ser
encontrado nas indústrias químicas, é amplamente utilizado no tratamento de efluentes líquidos
e gasosos. O tubo empacotado, que é um reator catalítico heterogêneo, é usado para converter
espécies B por meio da reação química:
B −→ produtos; RB = kCB
(
mols
tempo.volume do leito
)
em produtos sob condições (assumidas como sendo) isotérmicas. A difusão ao longo do compri-
mento axial é controlada por uma expressão do tipo Lei de Fick, de tal modo que, em paralelo
com o transporte por convecção devido a velocidade superficial v0, há também um fluxo do tipo
difusivo representado pela relação Fickiana, atuando ao longo da direção axial (coordenada z):
JE = −DE dCB
dz
(
mol
área.tempo
)
Para uma cinética linear, a aplicação da lei da conservação da massa (balanço de massa) aplicada
à espécie B, ao longo de um tubo com área de seção transversal A, resulta em:
−v0dCB
dz
− dJE
dz
−RB = 0
Legenda: DE = difusividade; k = constante da taxa da reação; v0= velocidade axial.
(a) Mostre que a substituição do fluxo e da expressão da taxa dados acima, na lei da conservação
da massa, resulta em:
DE
d2CB
dz2
− v0dCB
dz
− kCB = 0
(b) Classifique adequadamente, quanto ao tipo, à ordem, à linearidade e ao tipo de coeficiente,
a equação anterior.
(c) Mostre que a solução desta equação é dada por:
CB = exp(αz)[A exp(βz) +B exp(−βz)]
na qual:
α = v02DE , β =
√
1
4
(
v0
DE
)2
+
(
k
DE
)
4. (a) Resolva, usando o método dos coeficientes indeterminados, a seguinte EDO:
d2y
dx 2
− 8dy
dx
+ 16y − 6xe4x = 0
(b) Mostre que as funções que compõem a solução complementar são linearmente independentes.
5. Resolva, usando o método da variação dos parâmetros, o seguinte PVI:
2
d2y
dx 2
− 4dy
dx
+ 2y = 2xex , y(0) = 0, e y′(0) = 1.
6. Para cada sentença abaixo, marque V se for Verdadeira e F se for Falsa.
a) ( ) A EDO 2y′′ − 3y′ + y = 0 possui como raízes da equação característica r = 1, 1/2.
104
C. Lista de exercícios n. 3: EDOs de 2ª ordem, solução em série de potências de EDOs
b) ( ) A solução da EDO anterior é y = c1e
x + c2x
3/2e−x
c) ( ) A EDO y′′ + 2y′ + y = 0 é de 2ª ordem, não homogênea e de coeficientes constantes.
A função y1(x) = c1xe
x
é solução dessa EDO.
d) ( ) Se W (y1,y2) 6= 0(W é o wronskiano) então y1 não possui dependência linear com y2 e
y1 + y2 pode ser solução geral de uma EDO de 2ª ordem.
e) ( ) Se o wronskiano w de f e g é 3e4t e se f(t) = e2t, então g(t) = 3t+ e2t.
7. Encontre o valor de A para que y(x) = [c1cos(2x)+c2sen(2x)]e
−3x
seja solução de y′′+Ay′+13y =
0 .
8. Considere a EDO a seguir:
d2y
dx 2
− 2dy
dx
+ y = 0
Assumindo que a série de potências seja solução da referida EDO:
y =
∞∑
n=0
anx
n
(a) Mostre que a substituição de y, y' e y� na EDO conduz a:
∞∑
n=0
[an+2(n+ 2)(n+ 1)− 2(n+ 1)an+1 + an]xn = 0
(b) Mostre que a relação de recorrência é dada por (onde r são as raízes da equação característica
da EDO):
an = a0
rn
n!
(c) Mostre, finalmente, que a solução da EDO é então:
y(x) = a0
∞∑
n=0
xn
n!
+ b0x
∞∑
n=0
xn
n!
= a0 exp(x) + b0x exp(x)
(d) A solução anterior poderia ter sido alcançada através de outro método? Se sim, qual?
9. Considere a seguinte EDO de 2ª ordem e as respectivas condições inicias:
y′′ + 2xy′ + (1 + x2)y = 0, com , y(0) = 3 e y′(0) = −1
Para cada sentença abaixo, marque V se for Verdadeira e F se for Falsa.
a) ( ) Pode-se usar a série de potências y(x) =
∞∑
n=0
anx
n
, cuja expansão é em x0 = 0, para
resolver a EDO, pois o ponto de expansão considerado é um ponto ordinário.
b) ( ) Substituindo a série e suas respectivas derivadas na EDO, obtem-se:
∞∑
n=0
(n+ 2)(n+ 1)an+2x
n +
∞∑
n=1
2nanx
n +
∞∑
n=0
anx
n +
∞∑
n=2
an−2xn = 0
c) ( ) Da relação anterior obtem-se: 2a2 +a0 = 0, 6a3 + 3a1 = 0 e (n+ 2)(n+ 1)an+2 + (2n+
1)an + an−2 = 0 (para n > 2).
d) ( ) Da relação de recorrência válida para encontrar o termo genérico da solução em séries,
tem-se: a4 =
a0
8 ,a5 =
a1
8 ,a6 =
−a0
48 ,a7 =−a1
48 ,...
e) A solução em séries de potências para a EDO é: y(x) = 3(1− 1
2
x2 +
1
8
x4− 1
48
x6 +
1
384
x8−
...)− 1(x− 1
2
x3 +
1
8
x5 − 1
48
x7 +
1
384
x9 − ...)
105
C. Lista de exercícios n. 3: EDOs de 2ª ordem, solução em série de potências de EDOs
10. Considere a seguinte EDO: x2y′′ − 3xy′ + 4y = x2 lnx , x > 0
(a) Verifique que as funções dadas y1 e y2 satisfazem a equação homogênea associada: y1(x) =
x2, y2(x) = x
2 lnx;
(b) Mostre que as funções dadas y1 e y2 formam um conjunto fundamental de soluções da
equação homogênea associada;
(c) Mostre que Y (x) = 16x
2(lnx)3 pode ser uma solução particular da equação não homogênea
dada;
(d) A solução geral da EDO dada é:
11. Considere o comportamento do sistema mostrado neste vídeo: http://www.youtube.com/watch?
v=T7fRGXc9SBI.
(a) Que tipo de equação diferencial modela matematicamente o comportamento deste sistema?
(b) Faça uma breve pesquisa e escreva a EDO que descreve este comportamento, destacando as
hipóteses usadas no seu desenvolvimento;
(c) Em termos de EDO, qual a diferença entre o sistema mostrado no vídeo anterior e este outro
sistema? http://en.wikipedia.org/wiki/File:Oscillatory_motion_acceleration.ogv
12. Quais os objetivos dos comandos do Maxima a seguir?
eq1:x^2*'diff(y,x,2)+5*x*'diff(y,x)+4*y-x;
eq2:ode2(eq1,y,x);
eq3:ic2(eq2,x=0,y=0,'diff(y,x)=0.5);
eq3:ic2(eq2,x=1,y=0,'diff(y,x)=0.5);
eq4:rhs(eq3);
wxplot2d(eq4,[x,0.5,5],[ylabel,"y"]);
eq5:taylor(eq4,x,2,1);
eq6:taylor(eq4,x,2,2);
eq7:taylor(eq4,x,2,3);
eq8:taylor(eq4,x,2,4);
eq9:taylor(eq4,x,2,5);
wxplot2d([eq4,eq5,eq6,eq7,eq8,eq9],[x,0.5,5],[y,-0.5,1],[legend, false],[ylabel,"y"]);
13. Em que os sistemas de 2ª ordem descritos nos slides 4 e seguintes diferem dos exemplos dos slides
anteriores? http://goo.gl/fGvOuD [link para arquivo do Moodle, na forma curta do Google]
106
D. Lista de exercícios n. 4: solução de EDOs
por transformada de Laplace, sistemas de
EDOs
�
�
�
�
Universidade Federal de Uberlândia � Faculdade de Engenharia Química
Programa de Pós-Graduação em Engenharia Química � Métodos Matemáticos em Engenharia Química
Prof. Adilson J. de Assis
1. Resolva as seguintes EDOs por transformada de Laplace:
(a)
dy
dt
+ ky = 1, com y(0) = 1
(b)
d2y
dt2
+ 2
dy
dt
+ 2y = 0, sendo que em t = 0, y(t) = 1 e dydt = −1
(c) t
d2y
dt2
+
dy
dt
+ 4ty = 0, sendo que em t = 0, y(t) = 3 e dydt = 0. Resposta: y(t) = 3J0(2t)
2. Considere um reator de mistura perfeita (CSTR) isotérmico e com holdup constante (em outras
palavras, com volume constante) no qual ocorre uma reação de 1ª ordem, irreversível. A equação
que descreve como a concentração dentro do reator, CA, varia com o tempo é:
(acúmulo) = (entra) - (sai) + (geração) - (consumo) ≡ [massa/tempo]
d(V CA)
dt
= FCAo(t)− FCA + 0− V kCA
ou
dCA
dt
+
(
F
V
+ k
)
CA =
(
F
V
)
CAo(t)
para a qual CA(0) = 0.
(a) Sendo τ = V/F (tempo de residência no reator), mostre que:
G(s) ≡ CA(s)
CAo(s)
=
1/τ
s+ k + 1/τ
=
(
1
1+kτ
)
(
1
k+1/τ
)
s+ 1
=
Kp
τps+ 1
G(s) ≡ razão entre a saída e a entrada do processo, no domínio de Laplace, chamada de
função de transferência do processo;
Kp ≡ ganho do processo no estado estacionário;
τp ≡ constante de tempo do processo;
(b) Suponha que uma perturbação degrau de magnitude CAo seja aplicada ao processo:
CAo(t) =
{
CAo , t ≥ 0
0 , t < 0
Sabendo que CA(s) = G(s).CAo(s), mostre como a perturbação degrau afeta CA(t). Res-
posta: CA(t) = KpCAo(1− e−t/τp)
(c) Esboçe qualitativamente CA(t) versus t.
107
D. Lista de exercícios n. 4: solução de EDOs por transformada de Laplace, sistemas de EDOs
Figura D.0.1.: Dois tanques aquecidos.
3. Considere o sistema formado por dois tanques aquecidos, como mostrado na Figura D.0.1. Atra-
vés do balanço de energia, pode-se mostrar que as temperaturas T1 e T2 são dadas por:
ρCpV1
dT1
dt = ρCpF (T0 − T1) +Q1
ρCpV2
dT2
dt = ρCpF (T1 − T2)
Sendo: F = 90 ft3/min; Cp = 0,6 Btu/(lbm ºF); ρ = 40 lbm/ft
3
; V1 = 450 ft
3
; V2 = 90 ft
3
.
Encontre a expressão que fornece como T2(s) varia em função de T0(s) e Q1(s), na forma:
T2(s) = GL(s).T0(s) +GM (s).Q1(s)
Faça as considerações que forem necessárias.
4. A concentração de um reagente A varia dentro de uma partícula de catalisador cilíndrica
1
de
acordo com a EDO:
r
d2CA
dr2
+
dCA
dr
− r(k/DA)CA = 0
Sujeita a: em r = R, CA = CAs; em r = 0 [centro do catalisador],
dCA
dr = 0 [condição de
simetria para a concentração]. Encontre CA(r) utilizando transformada de Laplace, seguindo o
procedimento abaixo:
(a) Escreva a EDO na forma adimensional:
x
d2y
dx2
+
dy
dx
− xy = 0
com: y = CACAs e x = r
√
k/DA
(b) Aplique TL e encontre
Y (s) =
A√
s2 − 1
(c) Mostre que:
y(x) =
CA(r)
CAs
= AI0
(
r
√
k
DA
)
1
O catalisador usado na oxidação do SO2 a SO3, que é uma das etapas da produção de ácido sulfúrico, pode ser
aproximado por um cilindro, por exemplo!
108
D. Lista de exercícios n. 4: solução de EDOs por transformada de Laplace, sistemas de EDOs
Figura D.0.2.: Campo das direções de um circuito elétrico.
(d) Finalmente, mostre que:
CA(r)
CAs
=
I0
(
r
√
k
DA
)
I0
(
R
√
k
DA
)
(e) Usando: k/DA = 0, 5; 1; 1, 5, com R = 2, faça na mesma figura as três curvas de
CA(r)
CAs
versus r e interprete fisicamente os resultados;
5. Os autovalores e os autovetores da matriz
A =
[−4, 0 4, 0
−1, 6 1, 2
]
são (assinale a alternativa correta):
(a) λ1 = 2, λ2 = 0, 8;x
(1) = [1 0, 5]T , x(2) = [1 0, 8]T
(b) λ1 = −2, λ2 = −0, 8;x(1) = [2 1]T , x(2) = [1 0, 8]T
(c) λ1 = −4, λ2 = −1, 6;x(1) = [4 2]T , x(2) = [2 1, 6]T
(d) λ1 = 4, λ2 = 1, 6;x
(1) = [4 2]T , x(2) = [2 1, 6]T
(e) λ1 = −2, λ2 = −0, 8;x(1) = [2 1, 6]T , x(2) = [1 0, 8]T
6. Considere um circuito elétrico onde se deseja encontrar as correntes I1 e I2 dadas pelas equações
a seguir:
I ′1 = −4I1 + 4I2 + 12
I ′2 = −1, 6I1 + 1, 2I2 + 4, 8
Para cada sentença abaixo, marque V se for Verdadeira e F se for Falsa.
a) ( ) O equações formam um sistema de EDOs lineares, com coeficientes constantes e não
homogêneo;
109
D. Lista de exercícios n. 4: solução de EDOs por transformada de Laplace, sistemas de EDOs
Figura D.0.3.: Sistema massa-mola, com dois corpos e três molas.
b) ( ) Analisando a Figura D.0.2, percebe-se que independente da condição inicial do sistema,
para um tempo suficientemente grande, I1 = 3 e I2 = 0;
c) ( ) O ponto (3,0) no plano de fase mostrado na Figura D.0.2 é um nó;
d) ( ) Uma solução geral para o sistema é: J =
[
I1
I2
]
=
[
2c1 e
−2t +c2 e−0,8t +3
c1 e
−2t +0, 8c2 e−0,8t
]
e) ( ) Se ambas as correntes forem nulas no início, então: c1 = −4 e c2 = 5
7. Considere o problema de valor inicial (PVI) a seguir:
y′′ − y = t, y(0) = 1, y′(0) = 1
Assinale a alternativa correta:
(a) A transformada de Laplace é um método poderoso de resolução de EDOs, mas não é capaz
de resolver EDPs.
(b) Uma EDO no domínio original se transforma numa equação algébrica linear no domínio de
Laplace, a qual deve ser invertida para o domínio original.
(c) A transformada de Laplace da EDO anterior é: s2Y − sy(0)− y′(0)− Y = s
(d) A solução da equação algébrica no domínio de Laplace, do PVI anterior, é: Y = 1s−1 +
1
s2−1 − 1s
(e) A solução do PVI anterior é: y(t) = et + senh(t)− t
8. A Figura D.0.3 mostra um sistema massa-mola composto por 3 molas de constante k e dois
corpos, de igual massa, suspensos pela ação das molas. O deslocamento com o tempo do corpo
1 é y1 e do corpo 2, y2. Considere o código do Maxima a seguir.
kill(all);
eq1:'diff(y1(t),t,2)=-k*y1(t)+k*(y2(t)-y1(t))$
eq2:'diff(y2(t),t,2)=-k*(y2(t)-y1(t))-k*y2(t)$atvalue(y1(t),t=0,1)$
atvalue(y2(t),t=0,1)$
atvalue('diff(y1(t),t),t=0,sqrt(3*k))$
atvalue('diff(y2(t),t),t=0,-sqrt(3*k))$
assume(k>0);
sol1:desolve([eq1,eq2],[y1(t),y2(t)]);
sol2:ratsimp(sol1);
sol3:sol2[1];
110
D. Lista de exercícios n. 4: solução de EDOs por transformada de Laplace, sistemas de EDOs
sol4:sol2[2];
wxplot2d([ev(rhs(sol3),k=1),ev(rhs(sol4),k=1)],[t,0,15],
[style, [lines,3,1], [lines,2,2]],[legend, "eq6", "eq7"],
[xlabel, "t"], [ylabel, "y1,y2"]);
(a) A partir do exame do código do Maxima, reconstitua o PVI que está sendo resolvido;
(b) Explique, linha por linha, o que está sendo resolvido, interpretando convenientemente o(s)
resultado(s);
9. A transformada de Laplace inversa de
F (s) =
4s+ 10
s2 + 6s+ 8
é (assinale a alternativa correta):
(a) f(t) = e2t + e−4t
(b) f(t) = e−2t + 3e−4t
(c) f(t) = e2t + 3e4t
(d) f(t) = 3e−2t + e−4t
(e) f(t) = e−2t − 3e4t
10. Considere o seguinte PVI:
y′′ + 3y′ + 2y = sen(2t), com y(0) = 2 e y′(0) = −1
Para cada sentença abaixo, marque V se for Verdadeira e F se for Falsa.
a) ( ) O PVI anterior pode ser resolvido pelo método da transformada de Laplace.
b) ( ) Aplicando transformada de Laplace ao PVI, obtem-se: s2Y (s) − 2s + 1 + 3[sY (s) −
2] + 2Y (s) = 2
s2+4
c) ( ) Resolvendo a equação algébrica do item anterior, cuja incógnita é s, encontra-se:
s = 2s
3+5s2+8s+22
(s−4)(s2+3s+2)
d) ( ) Para encontrar a solução do PVI no domínio original, deve-se fazer a transformada de
Laplace inversa da solução no domínio de Laplace. Para tal, deve-se expandir a solução no
domínio de Laplace em frações parciais, neste caso em particular.
e) ( ) A solução deste PVI não contém o termo:
−1
20 cos(2t)
111
E. Lista de exercícios n. 5: EDPs
�
�
�
�
Universidade Federal de Uberlândia � Faculdade de Engenharia Química
Programa de Pós-Graduação em Engenharia Química � Métodos Matemáticos em Engenharia Química
Prof. Adilson J. de Assis
1. Considere o exemplo 4 do resumo de EDPs (determinação da difusividade térmica). Quanto
tempo uma almôndega congelada num freezer, de 5 cm de diâmetro, na qual se usou ovo na sua
formulação, necessita ficar em óleo de soja aquecido a 170 ºC, para garantir que uma possível
contaminação com salmonela seja eliminada? Faça o cálculo pela eq. 6 e pelo gráfico e compare
os resultados.
2. Considere uma aleta cilíndrica, de raio R e comprimento L na qual o calor se propaga por condu-
ção nas direções axial e radial. Uma das extremidades está imersa num líquido de temperatura
constante T1 e o meio que envolve a aleta está a uma temperatura constante T0. O modelo
matemático deste sistema é:
k
[
1
r
∂
∂r
(
r
∂T
∂r
)
+
∂2T
∂x2
]
= 0 (E.0.1)
Com as seguintes condições de contorno:
r = 0, ∂T∂r = 0
r = R, −k ∂T∂r = hG(T − T0)
x = 0, T = T1
x = L, ∂T∂x = 0
(E.0.2)
a) Classifique a EDP dada;
b) Interprete fisicamente a EDP e as condições de contorno dadas, incluindo nesta interpretação
um esquema físico do sistema;
3. Considere a equação diferencial parcial que descreve o comportamento difusivo em um meio
semi-infinito:
∂c
∂t
= D
∂2c
∂x2
(E.0.3)
Com as seguintes condições de contorno:
t = 0, c = 0, x ≥ 0
x = 0, D
(
∂c
∂x
)
x=0
= k (constante), t > 0
x → ∞, c = 0, t ≥ 0
(E.0.4)
a) Aplique transformada de Laplace na EDP, com respeito à variável tempo; i.e., obtenha a
seguinte EDO no domínio de Laplace:
sc(x, s)− c(x, 0) = Dd
2c(x, s)
dx2
(E.0.5)
b) Resolva a EDO no domínio de Laplace, obtendo:
c(x, s) = A1 exp
[(√
s
D
)
x
]
+A2 exp
[
−
(√
s
D
)
x
]
(E.0.6)
112
E. Lista de exercícios n. 5: EDPs
c) Aplique a TL às condições de contorno e após isto use tais condições de contorno já no
domínio de Laplace e encontre:
c(x, s) =
−k
sD
(
D
s
)1/2
exp
[
−
(√
s
D
)
x
]
(E.0.7)
d) Usando a inversa da transformada de Laplace, encontre:
c(x, t) =
−k
D1/2
[
2
(
t
pi
)1/2
exp
(−x2
4Dt
)
+
−x√
D
erfc
(
x
2
√
Dt
)]
(E.0.8)
A função erfc é a função erro complementar:
(%i1) wxplot2d(erfc(x),[x,-3,3]);
limit(erfc(x),x,inf);
limit(erfc(x),x,-inf);
erfc(0);
(%t1)
(%o1)
(%o2) 0
(%o3) 2
(%o4) 1
e) Utilizando o Maxima, faça um gráfico bidimensional de c(x) para vários t fixos; adote valores
razoáveis para as constantes envolvidas. Anexe o scrip do Maxima na sua resposta.
4. A temperatura numa barra de cobre que tem as laterais isoladas, 80 cm de comprimento, é
dada por u(x, t). O perfil inicial de temperatura é 100 sen(pi x/80) ºC. Dados físicos para o
cobre: densidade = 8,92 g/cm
3
; calor específico = 0,092 cal/(g.ºC); condutividade térmica =
0,95 cal/(cm.s.ºC).
(a) Encontre u(x, t);
(b) Quanto tempo demorará para que a temperatura máxima na barra caia para 50 ºC? Res-
posta: 388 s ( 6,5 min);
(c) Faça no Maxima e anexe o script usado: (1) um gráfico 3D de u(x, t); (2) as curvas de
nível u(x) para vários valores de t; (3) as curvas de nível u(t) para vários valores de x
(0 ≤ x ≤ 80 cm);
113
E. Lista de exercícios n. 5: EDPs
5. Considere o fluxo de água em um meio poroso, como areia, em um aquífero
1
. O fluxo de água é
bombeado por uma cabeça hidráulica, que é uma medida da energia potencial da água acima do
aquífero. Seja R: 0 < x < a, 0 < z < b uma seção vertical do aquífero. Em um meio homogêneo
uniforme, a cabeça hidráulica u(x, z) satisfaz a equação de Laplace em R:
uxx + uzz = 0
Se a água não puder fluir através dos lados e do fundo de R, então as condições de contorno serão
ux(0, z) = 0, ux(a, z) = 0, 0 ≤ z ≤ b
uz(x, 0) = 0, 0 ≤ x ≤ a
Finalmente, suponha que a condição de contorno em z = b é
u(x, b) = b+ αx, 0 ≤ x ≤ a
onde α é a inclinação do aquífero. Veja a Figura 5.
Figura E.0.1.: Esquema 2D de um aquífero.
(a) Encontre u(x, z) que satisfaz a equação de Laplace e as condições de contorno dadas;
(b) Sejam a = 1000, b = 500 e α = 0, 1. Desenhe curvas de nível da solução em R, ou seja,
desenhe algumas curvas de nível de u(x, z);
(c) A água flui dentro do aquífero ao longo de caminhos em R ortogonais às curvas de nível de
u(x, z) no sentido de u decrescente. Faça gráficos de alguns dos caminhos de fluxo.
Resposta do item (a):
u(x, z) = b+
αa
2
− 4αa
pi2
∞∑
n=1
cos[(2n− 1)pix/a] cosh[(2n− 1)piz/a]
(2n− 1)2 cosh[(2n− 1)pib/a] (E.0.9)
1
Um aquífero é uma formação ou grupo de formações geológicas que pode armazenar água subterrânea. São rochas
porosas e permeáveis, capazes de reter água e de cedê-la. Esses reservatórios móveis aos poucos abastecem rios e
poços artesianos e podem ser utilizadas pelo homem como fonte de água para consumo. É um reservatório natural
em que a água preenche por completo o volume de vazios do mesmo, até o solo ficar saturado. Geralmente é limitado
superior e inferiormente por materiais impermeáveis, em que a água é sujeita a uma pressão hidrostática superior à
pressão atmosférica. Na América do Sul, o aquífero mais famoso é o de Guarani, um imenso aquífero que abrange
partes dos territórios do Uruguai, Argentina, Paraguai e principalmente Brasil, ocupando 1 200 000 km
2
.
114
F. Revisão de matemática básica e de cálculo
diferencial e integral: uma abordagem
prática com o Maxima
115
Revisão de Matemática Básica e Cálculo Diferencial e Inte-
gral: uma abordagem prática com o Maxima
por Adilson J. de Assis
ajassis@ufu.br
• (x± y)2=x2± 2x y+ y2
(%i1) expand((x+y)^2)
(%o1) y2+2 x y+ x2
• x2− y2=(x− y) (x+ y)
(%i3) factor(x^2-y^2)
(%o3) −(y− x) (y+ x)
• Fatorial de n : n! = 1.2.3.� . n
(%i5) n!;
(%o11) n!
(%i12) 5!
(%o12) 120
• Para um triângulo retângulo:
* Teorema de Pitágoras: OB 2=OA 2+
AB 2
* sen2α+ cos2α=1(%i13) trigsimp(sin(x)^2+cos(x)^2);
(%o13) 1
• ap. aq= ap+q
(%i14) a^p*a^q
(%o14) aq+p
• ap/aq= ap−q
(%i15) a^p/a^q
(%o15) ap−q
• (ap)q= apq
(%i18) radcan((a^p)^q)
(%o18) apq
• a0=1, a� 0
(%i19) a^0
(%o19) 1
• 0a=0, a� 0
(%i20) 0^a
(%o20) 0
• 0
5
= 0, mas
5
0
é indefinido
(%i21) 5/0
expt: undefined: 0 to a negative
exponent.
-- an error. To debug this try:
debugmode(true);
• a−p=1/ap, a� 0
(%i22) a^(-p)
(%o22)
1
ap
• (a b)p= ap bp
(%i31) exp1:(a*b)^p;
(%o33) (a b)p
(%i34) radcan(%)
(%o34) ap bp
•
(
a
b
)−n
=
(
b
a
)n
=
bn
an
(%i39) (a/b)^(-n)
(%o39)
1( a
b
)n
(%i41) radcan(%)
(%o41)
bn
an
• a
n
m=
(
a
1
m
)n
=(an)
1
m
(%i45) (a^(1/m))^n
(%o45)
(
a
1
m
)n
(%i46) radcan(%)
(%o46) a
n
m
• an√ = a1/n
(%i47) a^(1/n)
(%o47) a
1
n
(%o48) a
1
n
1
• amn√ = am/n
• a/bn√ = an√ / bn√
• Cuidado!!!
* (a+ b)2� a2+ b2
* a2+ b2
√
� a+ b
*
3 a+ b
3 a+ c
�
b
c
*
d pi3
d x
� 3 pi2
• Função logarítmica� y= logax.
• logaMN = logaM + logaN
(%i61) log(%e) /* ATENÇÃO: log no
Maxima é ln */
(%o61) 1
(%i62) log(M*N)
(%o62) log (MN)
(%i64) %,logexpand=super
(%o64) log (N)+ log (M)
• logba= logca/logcb
(%i59) log(5)/log(10) /* logaritmo de
5 na base 10 */
(%o59)
log (5)
log (10)
(%i60) %, numer
(%o60) 0.69897000433601
(%i61)
• logaM/N = logaM − logaN
• logaM p= p logaM
• blogb x= x
(%i67) %e^log(x)
(%o67) x
• i= −1√
(%i68) sqrt(-1)
(%o68) i
• i2=−1
(%i70) %i^2
(%o70) −1
• −a√ = i a√ se a≥ 0
(%i73) sqrt(-a)-%i*sqrt(a)
(%o73) −a√ − i a√
(%i75) radcan(%)
(%o75) 0
• eiθ= cos θ+ i senθ
(%i76) %e^(%i*\theta)
(%o76) eiϑ
(%i77) rectform(%)
(%o77) i sin (ϑ)+ cos (ϑ)
• e−iθ= cos θ− i senθ
• sen θ= e
iθ− e−iθ
2 i
(%i81) sin(theta)
(%o81) sin (ϑ)
(%i82) exponentialize(%)
(%o82) − i (e
iϑ− e−iϑ)
2
• cos θ= e
iθ+ e−iθ
2
(%i83) cos(\theta)
(%o83) cos (ϑ)
(%i84) exponentialize(%)
(%o84)
eiϑ+e−iϑ
2
• Seno hiperbólico de x : senhx= e
x− e−x
2
• Cosseno hiperbólico de x : coshx= e
x+ e−x
2
• Equação quadrática ou do 2o grau:
a x2+ b x+ c=0;
raízes: x=
−b± b2− 4 a c
√
2 a
.
Se a, b e c são reais e se D= b2− 4 a c é o
discriminante, as raízes da equação são:
i). reais e desiguais se D> 0
ii). reais e iguais se D=0
iii). conjugadas e complexas se D< 0
• lim
x→a
b f(x)= b lim
x→a
f(x)
• lim
x→a
[f(x) + g(x)] = lim
x→a
f(x) + lim
x→a
g(x)
• lim
x→a
[f(x) g(x)] = lim
x→a
f(x) · lim
x→a
g(x)
• lim
x→a
f(x)
g(x)
=
limx→a f(x)
limx→a g(x)
, se lim
x→a
g(x)� 0
• lim
x→a
[f(x)]n= [ lim
x→a
f(x)]n
• lim
x→∞
ex=∞ e lim
x→−∞
ex=0
• lim
x→∞
ln (x)=∞ e lim
x→0+
ln (x) =−∞
2
• Descontinuidades:
* Removível: y= senx/x, quando x=0
* Infinita: y=1/x, quando x=0
* Salto: y= 10/(1+ e1/x), quando
x = 0+ y = 0+ ou x = 0 y = 0 ou
x=0− y= 10
• Se y = P (x), a derivada de y ou de P(x) em
relação à x é definida como
d y
d x
= lim
h→0
f (x+h)− f(x)
h
=
lim
∆x→0
f (x+∆ x)− f(x)
∆x
• A derivada é também designada por y ′, df/dx
ou f ′(x). O processo de derivar/obter uma
derivada é chamado diferenciação.
• d
d x
(c)= 0
(%i85) diff(c,x)
(%o85) 0
• d
d x
(c x) = c
(%i86) diff(c*x,x)
(%o86) c
• d
d x
(c xn)=n c xn−1
(%i87) diff(c*x^n,x)
(%o87) c n xn−1
• d
d x
(u± v±w±
 )= d u
d x
± d v
d x
± dw
d x
±
• d
d x
(c v)= c
d v
d x
(%i88) depends(v,x);
(%o88) [v(x)]
(%i88) diff(c*v,x)
(%o89) c
(
d
d x
v
)
• d
d x
(u v)= u
d v
d x
+ v
d u
d x
(%i90) depends([u,v],x)
(%o90) [u(x), v(x)]
(%i91) diff(u*v,x)
(%o91) u
(
d
d x
v
)
+
d
d x
u v
• d
d x
(
u
v
)
=
v
d u
d x
− u d v
d x
v2
(%i92) depends([u,v],x)
(%o92) [u(x), v(x)]
(%i94) diff(u/v,x)
(%o94)
d
d x
u
v
−
u
(
d
d x
v
)
v2
• d
d x
(un)=nun−1 d u
d x
(%i95) depends(u,x)
(%o95) [u(x)]
(%i96) diff(u^n,x)
(%o96) nun−1
(
d
d x
u
)
• Regra da cadeia: d y
d x
=
d y
d u
d u
d x
(%i104) depends(y,u);
(%o104) [y(u)]
(%i105) depends(u,x);
(%o105) [u(x)]
(%i106) diff(y,x)
(%o106)
d
d x
u
(
d
d u
y
)
• d u
d x
=
1
d x/d u
• d y
d x
=
d y/d u
d x/d u
• d
d x
senu= cosu
d u
d x
(%i107) depends(u,x)
(%o107) [u(x)]
(%i108) diff(sin(u),x)
(%o108) cos (u)
(
d
d x
u
)
• d
d x
cosu=−senu d u
d x
• d
d x
logau=
loga e
u
d u
d x
, a� 0, 1
• d
d x
lnu=
d
d x
loge u=
1
u
d u
d x
• d
d x
au= au ln a
d u
d x
(%i109) depends(u,x)
(%o109) [u(x)]
(%i110) diff(a^u,x)
(%o110) au log (a)
(
d
d x
u
)
3
• d
d x
eu= eu
d u
d x
• d
d x
uv =
d
d x
evlnu = evlnu
d
d x
[v ln u] =
v uv−1 d u
d x
+ uv lnu
d v
d x
• d
d x
senh u= cosh u
d u
dx
• d
d x
cosh u= senh u
d u
d x
• Derivadas superiores � segunda derivada:
d
d x
(
d y
d x
)
=
d2 y
d x2
= f ”(x)= y ”
(%i3) kill(all)
(%o0) done
(%i1) dependencies
(%o1) []
(%i2) ’diff(’diff(y,x),x) =
’diff(y,x,2)
(%o2)
d2
d x2
y=
d2
d x2
y
• Diferencial: d (u ± v ± w ± 
 ) = d u ± d v ±
dw±
• Diferencial: d (u v)= ud v+ v d u
• Seja f(x, y), uma função de duas variáveis x
e y. Então definimos a derivada parcial de
f(x, y) com relação a x , considerando y cons-
tante, por
∂ f
∂ x
= lim
∆x→0
f (x+∆ x, y)− f(x, y)
∆x
Igualmente a derivada parcial de f(x, y) com relação
a y , considerando x constante é definida por:
∂ f
∂ y
= lim
∆y→0
f (x, y+∆ y)− f(x, y)
∆ y
• A diferencial de f(x, y) (ou derivada total de
f) é definida como:
df=
(
∂ f
∂ x
)
y
d x+
(
∂ f
∂ y
)
x
d y
(%i3) depends(f,[x,y])
(%o3) [f(x, y)]
(%i4) diff(f)
(%o4)
d
d y
f del(y) +
d
d x
f del(x)
• Teorema de L’Hospital: Formas como 0/
0, ∞/∞, 0 × ∞ etc são chamadas de inde-
terminadas; se lim
x→a
f(x)
g(x)
=
0
0
ou
±∞
±∞ então
lim
x→a
f(x)
g(x)
= lim
x→a
f ′(x)
g ′(x)
. Exemplo: lim
x→0
sen(x)
x
=
lim
x→0
d sen(x)
d x
= lim
x→0
cos(x)
1
=1
• Se d y
d x
= f(x), então y é a função cuja derivada
é f(x) e é chamada de anti-derivada de f(x)
ou integal indefinida de f(x), designada por
y=
∫
f(x) d x
•
∫
a dx = a x + C [as constantes de integração
serão omitidas, mas estão implícitas]
(%i5) integrate(a,x)
(%o5) a x
•
∫
a f(x)dx= a
∫
f(x)dx
(%i6) ’integrate(a*f(x),x)
(%o6) a
∫
f(x) dx
•
∫
(u ± v ± w ± 
 ) dx =
∫
u dx ±
∫
v dx ±∫
w dx±
• Integração por partes:∫
u dv= u v−
∫
vdu
•
∫
un du=
un+1
n+1
, com n� −1
•
∫
d u
u
= lnu se u>0 ou ln (−u) se u<0 =ln |u|
(%i7) integrate(u^(-1),u)
(%o7) log (u)
•
∫
eudu= eu
•
∫
sen u du=−cos u
•
∫
cos udu= senu
• Transformações importantes:
*
∫
F (a x + b) dx =
1
a
∫
F (u) du onde
u= a x+ b
*
∫
F ( a x+ b
√
) dx =
2
a
∫
u F (u) du
onde u= a x+ b
√
*
∫
F ( a x+ bn
√
)dx=
n
a
∫
un−1F (u)du
onde u= a x+ bn
√
*
∫
F (eax) dx =
1
a
∫
F (u)
u
du onde
u= eax
*
∫
F (ln x) dx =
∫
F (u) eu du onde
u= lnx
4
• Métodos de integração1
* Fórmula direta: quando o integrando
pode ser transformado numa forma
facilmente integrável;
* Substituição trigonomética: aplica-
se em geral quando o integrando está na
forma de umradical;
* Substituição algébrica: aplica-se em
geral quando há funções da forma (a+
b x)1/n;
* Frações parciais: redução de funções
racionais na forma f(x)/g(x) para fra-
ções mais facilmente integráveis;
* Integração por partes
* Expansão em séries
•
∫
a
b
f(x)dx= g(x)|ab=g(b)− g(a)
(%i10) ’integrate(f(x),x,a,b)
(%o10)
∫
a
b
f(x) d x
•
∫
a
a
f(x)dx=0
•
∫
a
b
f(x)dx=−
∫
b
a
f(x)dx
•
∫
a
b
f(x) dx=
∫
a
c
f(x) dx+
∫
c
b
f(x) dx com
a<c< b
• Função Gama: Γ(n) =
∫
0
∞
tn−1 e−t dt com
n> 0
(%i11) gamma(n)
(%o11) Γ(n)
• Função Erro: erf(x)= 2
pi
√
∫
0
x
e−u2du
(%i12) erf(x)
(%o12) erf(x)
(%i13) erf(0.5)
(%o13) 0.52049987781304
• erf(x)= 2
pi
√
(
x− x
3
3.1!
+
x5
5.2!
− x
7
7.3!
+
)
• erf(−x)=−erf(x), erf(0)= 0, erf(∞)= 1
(%i14) erf(-x)
(%o14) −erf(x)
(%i15) erf(0)
1. Consultar a literatura de Cálculo Diferencial e Integral
para detalhes dos métodos.
(%o15) 0
(%i16) erf(inf)
(%o16) 1
• Série de Taylor [expansão de uma função
f(x) em torno de um ponto a]:
f(x) = f(a) + f ′(a) (x − a) +
f ′′(a) (x− a)2
2!
+ 
 +
f (n−1)(a) (x− a)n−1
(n− 1)! +
Rn
onde Rn é o resto após n termos. Se a =
0 a série é frequentemente chamada série de
Maclaurin
(%i20) taylor(f(x),x,a,5)
(%o20) f(a) +
(
d
d x
f(x)
∣∣∣∣
x=a
)
(x − a) +(
d2
d x2
f(x)
∣∣∣
x=a
)
(x− a)2
2
+(
d3
d x3
f(x)
∣∣∣
x=a
)
(x− a)3
6
+(
d4
d x4
f(x)
∣∣∣
x=a
)
(x− a)4
24
+(
d5
d x5
f(x)
∣∣∣
x=a
)
(x− a)5
120
+�
• ex=1+ x+ x
2
2!
+
x3
3!
+
 com −∞<x<∞
(%i21) taylor(%e^x,x,0,5)
(%o21) 1+ x+
x2
2
+
x3
6
+
x4
24
+
x5
120
+�
• ln (1 + x) = x − x
2
2
+
x3
3
− x
4
4
+ 
 com
−1<x< 1
• ln x =
(
x− 1
x
)
+
1
2
(
x− 1
x
)
2
+
1
3
(
x− 1
x
)
3
+
 com x≥ 1
2
• senx=x− x
3
3!
+
x5
5!
− x
7
7!
+
 com −∞<x<
∞
(%i23) taylor(sin(x),x,0,7)
(%o23) x− x
3
6
+
x5
120
− x
7
5040
+�
• cosx=1− x
2
2!
+
x4
4!
− x
6
6!
+
 com −∞<x<
∞
• Laplaciano: U = ∇2 U = ∇ · (∇ U) = ∂
2U
∂x2
+
∂2U
∂ y2
+
∂2U
∂ z2
(%i11) kill(all);
load(vect);
5
(%o0) done
(%o1)
/usr/share/maxima/5.32.1/share/vector/vect.mac
(%i6) depends(U,[x,y,z])
(%o6) [U(x, y, z)]
(%i22) laplacian(U)
(%o22) laplacian(U)
(%i23) express(%) /* laplaciano de U
*/
(%o23)
d2
d z2
U +
d2
d y2
U +
d2
d x2
U
(%i9) grad(U)$ express(%) /* gradiente
*/
(%o10)
[
d
d x
U ,
d
d y
U ,
d
d z
U
]
(%i13) div([U,U,U])$ express(%) /*
divergente */
(%o14)
d
d z
U +
d
d y
U +
d
d x
U
(%i15) curl([U,U,U])$ express(%) /* */
(%o16)
[
d
d y
U − d
d z
U ,
d
d z
U − d
d x
U ,
d
d x
U − d
d y
U
]
(%i19) [a,b,c]~[x,y,z] /* produto
cruzado */
(%o19) ([a, b, c], [x, y, z])
(%i20) express(%)
(%o20) [b z − c y, c x− a z, a y− b x]
• Equação diferencial de Bessel: x2 y ′′ + x y ′ +
(x2 − n2) y = 0 com n > 0. As soluções
desta equação são chamadas funções de Bessel
de ordem n.
• Funções de Bessel de primeira classe de ordem
n:
* Jn(x)=
∑
k=0
∞
(−1)k(x/2)n+2k
k!Γ (n+ k+1)
* J−n(x) =
∑
k=0
∞
(−1)k(x/2)2k−n
k!Γ (k+1−n)
* J−n(x) = (−1)nJn(x) n=0, 1, 2,
* J0(x)= 1− x
2
22
+
x4
22 · 42 −
x6
22 · 42 · 62 +
(%i25) bessel_j(0,0)
(%o25) 1
* J1(x) =
x
2
− x
3
22 · 4 +
x5
22 · 42 · 6 −
x7
22 · 42 · 62 · 8 +
(%i26) bessel_j(1,0)
(%o26) 0
* J0
′(x)=−J1(x)
(%i27) diff(bessel_j(0,x),x)
(%o27) −bessel˙j(1, x)
• Funções de Bessel de segunda classe de ordem
n: Yn(x)
(%i25) bessel_y(n,x)
(%o25) bessel˙y(n, x)
6
	Sumário
	Lista de Figuras
	Introdução: modelagem matemática, engenharia e equações diferenciais
	Revisão de Matemática Básica e Cálculo Diferencial e Integral
	Breve apresentação do software livre Maxima
	Iniciando o uso do Maxima
	Processando um comando
	Manipulação de variáveis
	Funções matemáticas elementares no Maxima
	Operações básicas com matrizes
	Solução de sistemas lineares
	Solução analítica de uma equação
	Solução numérica de uma equação
	Manipulação algébrica de equações e funções
	Simplificações diversas
	Gráficos bi e tridimensionais
	Limites e continuidade de funções
	Derivadas
	Derivadas simples e pontos críticos de uma função
	Derivadas parciais de uma função
	Integrais indefinidas, definidas e área abaixo de uma curva
	Programação no Maxima
	Instrução: if
	Instrução: for...end
	Instrução block
	Séries
	Série de Taylor
	Série de Fourier
	Regressão Linear
	Campo de direções no Maxima
	Exemplo 1
	Exemplo 2
	Exemplo 3
	Exemplo 4
	Função drawdf
	Solução analítica de equações diferenciais ordinárias (EDOs)
	Verificando se uma dada equação é solução de uma EDO
	Equações diferenciais de 1ª ordem
	Encontrando a solução geral de uma ED de 1ª ordem
	Encontrando a solução particular de uma ED de 1ª ordem
	EDO com solução implícita
	Solução de EDOs exatas
	Exemplo 1: variação da massa de um soluto com o tempo a volume constante
	Exemplo 2: lei do resfriamento de Newton
	Equações diferenciais lineares de 2ª ordem
	Solução geral
	Solução de PVIs
	Exemplo: difusão e reação de um componente gasoso em um líquido [exemplo 2.7]rice
	Solução de PVCs
	EDOs não homogêneas
	Exemplo: Dinâmica de um tanque agitado com aquecimento elétrico
	Wronkiano
	Pacote contrib_ode
	Solução usando série de potências para EDOs lineares de 2ª ordem
	Solução de EDOs por transformada de Laplace (EDOs simples e sistemas de EDOs)
	Transformada de Laplace
	Solução de EDOs por transformada de Laplace
	Sistemas de EDOs lineares de 1ª ordem no plano
	Autovalores e autovetores
	Retrato de fase de sistemas lineares
	Solução numérica de equações diferenciais ordinárias usando o método de Runge-Kutta
	Solução analítica de equações diferencias parciais (EDPs)
	Exemplos de EDPs na Engenharia Química
	Séries de Fourier
	Solução de EDPs por separação de variáveis
	Solução numérica de equações diferencias parciais (EDPs) usando o método das diferenças finitas
	EDP parabólica
	EDP hiperbólica
	EDP elíptica
	Referências Bibliográficas
	Lista de exercícios n. 1: Revisão de matemática básica e cálculo diferencial e integral
	Lista de exercícios n. 2: Classificação, campo de direções e EDOs de 1ª ordem
	Lista de exercícios n. 3: EDOs de 2ª ordem, solução em série de potências de EDOs
	Lista de exercícios n. 4: solução de EDOs por transformada de Laplace, sistemas de EDOs
	Lista de exercícios n. 5: EDPs
	Revisão de matemática básica e de cálculo diferencial e integral: uma abordagem prática com o Maxima