Aplicações Nas Ciencias

29
3 Modelagem Matemática
3.1 Modelo Matemático

Para se definir modelo matemático, tomaremos como base a definição de Bassanezi
(2002):
Modelo matemático é um conjunto de símbolos e relações matemáticas

que representam de alguma forma o objeto estudado. O modelo pode
ser considerado como uma síntese de reflexão sobre alguma parte da
realidade. Seu objetivo é explicar ou entender a situação estudada para,
eventualmente poder agir sobre ela e, mesmo as situações mais simples
fornecem motivações para uma iniciação científica (BASSANEZI, 2002,
p. 12).
Destacamos a definição do Professor Bassanezi, pois ele é um dos pioneiros no

estudo de Modelagem Matemática a qual tem se dedicado por muitos anos com publicações
e produções científicas, orientações e tem atuado como professor titular nos maiores centros
universitários brasileiros (CARVALHO, 2010).
Algumas etapas são necessárias para obter um modelo. Em sua obra Equações
Diferenciais Ordinárias: um curso introdutório, Bassanezi expõe sete etapas significativas:
(1) experimentação, (2) abstração, (3) formulação do modelo, (4) resolução, (5) validação,
(6) modificação e (7) aplicação.
Experimentação: para se compreender o problema, levanta-se dados a partir
de experimentos que por sua vez são bastante relevantes na estruturação, formulação e
possíveis modificações que podem vir a ser necessárias nos modelos e na sua consequente
validação por meio dos dados obtidos.
Abstração: nessa fase de estudo são levantadas hipóteses para serem testadas no
modelo e assim verificar as variáveis essenciais para sua evolução.
Formulação do Modelo: aqui a linguagem matemática entra em ação substi-
tuindo a linguagem usual anteriormente para a formulação das hipóteses. “A construção
do modelo segue de perto o uso de um dicionário que traduz as palavras chaves em alguma
estrutura matemática” (BASSANEZI, 2002, p. 13). Um exemplo inicial e de fácil verificação
é a variação de uma determinada população no tempo. Em linguagem matemática, temos
dP
a derivada de P (t) em relação ao tempo e que nesse caso é um tratamento de uma
dt
variação contínua. Também, podemos citar P (t2 ) − P (t1 ) que é uma variação discreta da
população em dois tempos distintos ou até mesmo a variação média dessa população nos
P (t2 ) − P (t1 )
mesmos dois tempos citados que é dada por .
t2 − t1
30 Capítulo 3. Modelagem Matemática
Se formularmos a hipótese de que a variação populacional é proporcional à população

dP
e se tratando, por exemplo, de uma variação contínua de P (t), então: = kP , em que k
dt
é uma constante de proporcionalidade.
Resolução: Como já citado anteriormente, a resolução de um modelo pode ser
dada analítica ou numericamente 1 dependendo do quão complexo é o fenômeno estudado
e modelado.
Validação: A comparação é feita por intermédio dos resultados obtidos com a
solução do modelo com os dados reais obtidos experimentalmente. Nesse ponto é necessário
verificar a proximidade dos dados com o que foi “projetado” inicialmente por meio das
hipóteses formuladas.
Modificação: Esta etapa depende da anterior, pois se a proximidade dos dados
não está de acordo com que foi “projetado” inicialmente, então novas variáveis devem
ser inserida no modelo ou até mesmo modificando sua lei de formação, iniciando, assim,
novamente o processo.
Aplicação: A aplicação do modelo permitirá um estudo da situação de acordo
com as variáveis selecionadas, permitindo, assim, o entendimento, a tomada de decisões,
previsões, como por exemplo, em função do tempo, fornecendo “ferramentas” para criar
estratégias de trabalho para determinar melhoras em diversas situações modeladas.
A partir dessa reflexão, acerca da Modelagem Matemática, podemos afirmar que
há uma tendência de crescimento da aplicabilidade em diversas áreas do conhecimento.
Bassanezi acrescenta:
a Matemática tem penetrado fortemente na Economia, Química, Biologia,

entre outras, na perspectiva da utilização de modelos, quase sempre
apoiados nos paradigmas que nortearam a Física – como as leis de
conservação e analogias consequentes. Outras áreas como Sociologia,
Psicologia, Medicina, Linguística, Música, e mesmo a História, começam
a acreditar na possibilidade de ter suas teorias modeladas por meio da
linguagem matemática (BASSANEZI, 1999, p. 13).
A partir das necessidades em um determinado momento, os profissionais dessas áreas

sentem a necessidade de uma explicação matemática de algum fenômeno de interesse ou da
satisfação, de algum modo, de uma necessidade humana. Isso nos leva a fazer suposições
a fim de construir uma estrutura básica do sistema a ser modelado. Se tais suposições
são consideradas suficientemente precisas, a consequência imediata disso são as equações
que regem o sistema do fenômeno estudado em questão. Um exemplo disso é que na
mecânica clássica a massa é considerada uma constante, porém na teoria da relatividade de
Einstein, a massa é considerada variável. Com isso em mente, Newton estudou corpos com
1
Reforçamos que nesse trabalho nos preocuparemos mais com as soluções analíticas, mas não descartamos
a grande importância das soluções numéricas estudadas na disciplina Cálculo Numérico.
3.1. Modelo Matemático 31
velocidades consideradas baixas, mas Einstein modelou o sistema considerando velocidades

próximas ou iguais à velocidade da luz gerando, assim, inconsistência acerca daquilo que é
estudado na teoria e o que é realmente observado.
Não podemos deixar de relatar que uma das partes consideradas mais difícil de
um modelo é o teste, pois é nesse momento que podem aparecer os possíveis erros e é
importante ter em mente as previsões de erros no modelo considerado. Se estivermos
trabalhando com dados empíricos, há uma grande necessidade de tomar cuidado porque o
fenômeno estudado pode ser muito específico para aquele sistema e é possível relatar a
realidade apenas daquele fenômeno e não apresentar uma generalização.
Se sabemos que é possível encontrar erros ao modelar um sistema, então podemos
ter uma base sólida para saber como lidar com eles. Caso estamos estudando um fenômeno
natural, por exemplo, o crescimento populacional de uma determinada bactéria, podemos
fazer previsões de possíveis erros ao conhecer o ambiente natural e sua variabilidade por
meio de estudos e pesquisas feitos anteriormente para esse específico sistema. Tendo em
mãos os dados dos testes realizados, podemos fazer ajustes no modelo para que ele possa
atender bem ao sistema e verificar, assim, se tais alterações atendem de forma eficaz os
dados obtidos anteriormente.
Com o que relatamos anteriormente, podemos notar que a interdisciplinaridade
da Matemática Aplicada, segundo Bassanezi (1999), nos permite utilizar das estruturas
matemáticas fora do seu próprio contexto para adentrar–se em um contexto alheio, isto é,
a Modelagem Matemática aplicada em outras ciências.
Além de Bassanezi, outro nome que se destaca nas pesquisas sobre Modelagem
Matemática é o de Dionísio Burak. Para esse autor,
A modelagem matemática constitui-se em um conjunto de procedimentos

cujo objetivo é construir um paralelo para tentar explicar, matematica-
mente, os fenômenos do qual o homem vive o seu cotidiano, ajudando-o
a fazer predições e a tomar decisões (BURAK, 1987, p. 21).
Isso quer dizer que um modelo matemático pode traduzir fenômenos que se asse-
melham a acontecimentos reais, sendo que o estudo de tais modelos permite que façamos
previsões futuras de situações que ocorrem no presente. Conforme cita Burak (1987) “a
adequação de um modelo pode ser julgada por seu êxito em ordenar os dados e, a partir
deles, fazer predições”, ou seja, quanto mais o modelo se aproximar da realidade mais
garantido será inferir sobre seus dados.
Burak (1987) apresenta, em sua dissertação, onde tem sido desenvolvido o uso da
Modelagem Matemática na educação: iniciação científica, aperfeiçoamento para o professor,
cursos regulares e pesquisa científica. A iniciação científica proporciona, com a modelagem
matemática, as aplicações da Matemática em outras áreas do conhecimento podendo
32 Capítulo 3. Modelagem Matemática
ser aprofundada a partir das áreas de pesquisa da disciplina em questão, por exemplo,
as Equações Diferenciais Ordinárias. O aperfeiçoamento para o professor proporciona a
especialização do docente para uma nova metodologia de ensino ou simplesmente renovar os
conhecimentos já adquiridos anteriormente a respeito da modelagem. Os cursos regulares,
por exemplo, de Cálculo Diferencial e Integral pode proporcionar um ambiente rico de
aplicações e modelos matemáticos em diversas áreas como Engenharia, Biologia, Física,
entre outras áreas do conhecimento, permitindo, assim, investigações por parte dos alunos
em sua própria área de estudo. E a pesquisa científica, de acordo com Burak (1987), é o
segmento que, por meio de projetos, faz largo uso dos modelos matemáticos em diversas
áreas do conhecimento e há um aumento significativo de programas de pós-graduação na
área da Modelagem Matemática.
Como forma de indagação, de acordo com Barbosa (2001), a Modelagem Matemática
propicia aos alunos oportunidades de questionamentos em diversas situações utilizando os
conceitos matemáticos sem a necessidade de procedimentos pré-fixados. O encaminhamento
fornecido aos estudantes é que apresentará o ambiente de aprendizagem por meio da
modelagem e as condições das ferramentas fornecidas é que possibilitará a estimulação e a
criatividade para solucionar as situações propostas.
Barbosa (2001) assume modelagem como “um ambiente de aprendizagem no qual
os alunos são convidados a indagar e/ou investigar, por meio da matemática, situações
oriundas de outras áreas da realidade”. Nesta pesquisa, em cada situação, estamos apre-
sentando as aplicações por meio de equações diferenciais para estudar problemas oriundos
de outras realidades além da matemática aplicada à própria matemática.
Com isso, a pesquisa será realizada a fim de investigar a utilização das Equações
Diferenciais Ordinárias para solucionar problemas modelados matematicamente em algumas
áreas do conhecimento e, para isso, utilizaremos as etapas descritas por Bassanezi (2002)
nas situações apresentadas.
33
4 Modelagem Matemática nas En-

genharias
Seria um tanto quanto injusto escolher uma definição única para Engenharia, mas,
para esse trabalho, uma “definição” apropriada é que Engenharia é uma ciência que utiliza
o conhecimento das Ciências Matemáticas e Naturais para desenvolver recursos e utilizar
de forma justa aquilo que é obtido na natureza em benefício da sociedade. Ou seja, uma
ciência que depende de outras ciências para cumprir com seu objetivo.
Apresentamos a descrição acima porque um curso de Engenharia tem uma carga
ampla das ciências exatas principalmente de Matemática. Nosso objetivo é mostrar Modelos
Matemáticos que são frequentemente utilizados para descrever situações dentro de cada
campo a fim de maximizar as ideias teóricas estudadas nos cursos de Cálculo e EDO
ligando-as contextualmente dentro das aplicações e desafios que surgirão dentro da carreira
profissional.
A construção e o estudo de um modelo são mais relevantes do que sua solução pois
envolvem os desafios proporcionados pela elaboração de hipóteses que levam até sua tese
final. Com isso, queremos destacar e até mesmo fazer uma ligação desses desafios ao que foi
citado no Capítulo 2 sobre as etapas necessárias para obter um Modelo Matemático, pois
um Engenheiro também experimenta, abstrai ideias, formula hipóteses, resolve problemas,
valida, modifica e, por fim, faz aplicações a partir de suas conclusões finais.
Os modelos que aqui serão apresentados são desafios e problemas além dos exercícios
propostos para ilustrar a teoria. Temos como pretensão fornecer aos futuros Engenheiros
subsídios de que os tópicos que são estudados em Cálculo Diferencial e Integral e Equações
Diferenciais Ordinárias são poderosas ferramentas para utilização em sua vida profissional
e não só acadêmica.
4.1 Aplicações das EDO na Engenharia Mecânica

Na Engenharia Mecânica um dos principais focos de estudo é analisar e estudar
os efeitos da força e do movimento nos projetos de sistemas em diversos ramos internos
dentro dessa própria Engenharia, como por exemplo, Estática, Dinâmica, Mecânica dos
Fluidos, Termodinâmica entre outras subáreas.
Um dos tópicos que tem uma carga horária grande de estudo em um curso de
Engenharia Mecânica é Mecânica Geral e de acordo com Kraige e Merian (1997), os alunos
encaram o curso como uma exigência difícil e com obstáculos acadêmicos desinteressantes
34 Capítulo 4. Modelagem Matemática nas Engenharias
porque focam-se muito no aprendizado por memorização, o que é muito diferente do

aprendizado a partir de fundamentos teóricos, pois com isso poderão atacar os problemas
da mecânica com raciocínios concisos a partir de formulação de hipóteses criando, assim,
uma forte ligação com a prática profissional. O ideal é motivar o aluno no aspecto de que
“a teoria pode apenas se aproximar do mundo real da mecânica, ao contrário da visão de
que o mundo real se aproxima da teoria.” (KRAIGE; MERIAN, 1997, citação de prefácio).
Isso vem ao encontro com o que foi citado anteriormente de que um modelo, matemático
ou não, pode ser considerado como uma reflexão sobre alguma parte da realidade.
Como primeiro exemplo, citaremos uma aplicação de EDO muito importante na
área da Mecânica que são as Vibrações de Sistemas. Exemplos conhecidos que envolvem
esse tipo de aplicação são: vibração de máquinas rotativas, por exemplo, ventiladores
industriais, máquinas operatrizes na indústria metalúrgica e até mesmo em vibrações com
grande escala a qual podemos citar os terremotos. Uma das soluções para amenizar os
problemas e danos decorrentes de tais vibrações é assentar os “corpos” sobre molas para
minimizar sua propagação.
De fato, nos concentraremos na vibração livre de partículas subdividida em duas
partes: (a) vibração livre não-amortecida e (b) vibração livre amortecida. Esses tipos de
vibrações são modelados com Equações Diferenciais Ordinárias de Segunda Ordem.
4.1.1 Vibração livre não-amortecida

Como forma introdutória e intuitiva, tomemos a Lei de Hooke que está relacionada
à elasticidade de determinados corpos, principalmente as molas. A força de um corpo
exercido na mola causa uma determinada deformação a partir de um comprimento de
equilíbrio. Desprezando o atrito, a força de tensão na mola é dada por:
F = −kx (4.1)
em que k é a constante elástica da mola e o sinal negativo significa a força contrária que a
mola faz para voltar ao seu comprimento de equilíbrio e x é o afastamento em relação ao
equilíbrio.
Analisando a Lei de Hooke, podemos transformá-la em uma Equação Diferencial
Ordinária. Sabendo que F = −kx e analisando a 2a Lei de Newton, sabendo que a aceleração
é a derivada segunda do espaço em relação ao tempo , então ma = −kx ⇒ mx = −kx.
Disso, podemos concluir que
mx + kx = 0 (4.2)
A expressão (4.2) é uma Equação Diferencial de Segunda Ordem e é conhecida

também como movimento harmônico simples e esse tipo de movimento é caracterizado por
oscilações do espaço em função do tempo (movimento da partícula sobre o eixo Ox). Como
4.1. Aplicações das EDO na Engenharia Mecânica 35
vimos no Capítulo 2, a solução dessa EDO Linear Homogênea de 2a Ordem é caracterizada

pelo Teorema 2.4.2. De fato, tomando (4.2), temos
mx + kx = 0
k
x + x = 0.
m
A equação característica é dada por
k
λ2 + =0
m

k k
cujas raízes são λ1 = i e λ2 = − i. Logo a solução geral da equação (4.2), de
m m
acordo com Teorema 2.4.2 do Capítulo 2, será
⎛ ⎞ ⎛ ⎞
k ⎠ k ⎠
x(t) = Acos ⎝ t + Bsen ⎝ t , A, B ∈ R (4.3)
m m
Podemos, na equação acima, fazer uma substituição conveniente. Ao analisarmos as

vibrações, recaímos no estudo das Frequências Circulares Naturais 1 . Da Física, sabemos
que

k k
ωn = ⇒ ωn2 =
m m
k
e substituindo em x + x = 0, então
m
x + ωn2 x = 0 (4.4)

k
em que ωn = é a velocidade angular constante cuja unidade é radianos por segundo.
m
Portanto, a equação (4.3) poderá ser escrita da seguinte forma
x = Acos (ωn t) + Bsen (ωn t) , A, B ∈ R (4.5)
Embora há vários tipos de molas, como molas não-lineares duras e moles, restrin-
giremos nosso estudo às molas lineares que é o caso já exemplificado anteriormente pela
Lei de Hooke. Se a força exercida pela massa m esticar ou comprimir esse tipo de mola,
ela exerce uma força restauradora −kx seja para a direita ou para a esquerda. Abaixo,
apresentaremos duas figuras que exemplificam a situação. Na primeira, a mola está em
equilíbrio e, na segunda, a mola sofre uma tensão (no caso, o esticamento da mesma).
1
Para um estudo mais detalhado das FCN’s recomendamos uma análise nas obras de Halliday, Walker
e Resnick (2010a) e Kraige e Merian (1997).
Figura 2 – Mola em equilíbrio.
Figura 3 – Mola após sofrer uma tensão (esticamento).
Para a análise do modelo de acordo com as etapas de modelagem descritas no

Capítulo 3, primeiramente tomamos a Lei de Hooke como hipótese para verificação das
variáveis necessárias para evoluirmos o modelo, por exemplo, a identificação da 2a Lei de
Newton, fase esta que chamamos de abstração.
A segunda etapa necessária foi a formulação do modelo. Identificadas as variáveis
necessárias, passamos a modelar a equação diferencial. Sabemos que a aceleração é a taxa
de variação instantânea da velocidade em função do tempo que, por sua vez, é a taxa de
variação instantânea do espaço em função do tempo e, portanto, temos que a aceleração é
a segunda derivada da função espaço. Fazendo as substituições necessárias, chegamos à
fórmula mx + kx = 0.
O próximo passo foi resolver o modelo por intermédio do conhecimento e processos
de resolução de equações diferenciais ordinárias. Esta etapa está caracterizada por resolução.
Utilizamos do Teorema 2.4.2 para obter a solução da equação diferencial.
A solução da equação diferencial (4.2) nos proporcionou avaliar a validade do
modelo de acordo com as hipóteses iniciais consideradas. A Lei de Hooke por si só não
considera o atrito no sistema, logo o dispositivo massa-mola sofrerá ação da força e ficará
oscilando até uma força externa parar o sistema e, como a solução da equação (4.2) é dada
por uma combinação linear de funções periódicas, então concluímos a etapa chamada de
validação.
Validado o modelo, a próxima etapa é a aplicação. Na subseção (4.1.3), apresentamos
a aplicação da vibração livre não-amortecida para verificação imediata do que o fenômeno
pode nos apresentar. Fornecemos tanto uma leitura algébrica quanto uma leitura gráfica
para uma visualização de diferentes ângulos do fenômeno.
4.1.2 Vibração livre amortecida

No estudo das vibrações livres não-amortecidas despreza-se o atrito na descrição do
fenômeno, porém “todo sistema mecânico possui algum grau inerente de atrito, que age
como um consumidor de energia mecânica” (KRAIGE; MERIAN, 1997, p. 422). Torna-se
complexo modelar matematicamente (com precisão) um sistema com forças de atrito
dissipativas.
Um exemplo que é encontrado em diversos sistemas para limitar as vibrações são
os amortecedores (ou pistões) com fluidos viscosos, por exemplo óleo ou ar. Nesses casos,
a equação (4.1) é acrescida de uma constante de proporcionalidade c conhecida como
coeficiente de amortecimento viscoso e como se trata de um novo movimento acrescido, ou
seja, do pistão dentro do cilindro, então existe uma velocidade (objeto percorrendo sobre o
eixo Ox) contrária à direção da massa e, portanto, a força sobre massa é −cx . A derivada
primeira utilizada é devido ao fato de que a velocidade em relação ao tempo é a derivada
do espaço x em relação ao tempo t.
Para esse estudo, a Segunda Lei de Newton nos fornece
mx = −cx − kx
mx + cx + kx = 0. (4.6)
Novamente, temos um fenômeno modelado matematicamente por uma Equação Diferencial

Ordinária de Segunda Ordem. Como fizemos anteriormente, introduziremos algumas
substituições convenientes. Novamente, temos
k
ωn2 =
m
e
c
c = 2mωn ζ ⇒ ζ =
2mωn
em que ζ é chamado de taxa de amortecimento que mede o grau de amortecimento.
Substituindo esses valores na equação (4.6), então
c k
mx + cx + kx = 0 ⇒ x + x + x=0
m m
x + 2ζωn x + ωn2 x = 0. (4.7)
A equação característica é dada por
λ2 + 2ζωn λ + ωn2 = 0
√ √
cujas raízes são λ1 = ωn −ζ + ζ 2 − 1 e λ2 = ωn −ζ − ζ 2 − 1 .
√
A análise dessa solução se dá pelo fator ζ 2 − 1 originando três categorias de
movimento amortecido:
1. Superamortecido. Quando ζ > 1 a equação característica fornece duas raízes reais λ1

e λ2 . Com efeito, a solução geral será
√ √
−ζ+ ζ 2 −1 ωn t −ζ− ζ 2 −1 ωn t
x(t) = Ae + Be . (4.8)
2. Criticamente amortecido. Quando ζ = 1 a equação característica fornece duas raízes

reais iguais λ1 = λ2 = −ωn t. Com efeito, a solução geral será
x(t) = e−ωn t (A + Bt) . (4.9)
3. Subamortecido. ζ < 1 a equação característica fornece duas raízes complexas, uma

√
conjugada da outra. Note que ζ 2 − 1 é negativo, então
√ 2 √ √
ζ − 1 = −1(1 − ζ 2 ) = i 1 − ζ 2 , em que i = −1.
Reescrevendo a equação (4.8) e fazendo os ajustes necessários conforme as manipula-
ções algébricas feitas acima, temos
√ √
−ζ+ ζ 2 −1 ωn t −ζ− ζ 2 −1 ωn t
x(t) = Ae + Be
√ √
−ζ+i 1−ζ 2 ωn t −ζ−i 1−ζ 2 ωn t
= Ae + Be
√ √
2 2
= Ae−ζωn t ei 1−ζ ωn t + Be−ζωn t e−i 1−ζ ωn t
e, por fim,
! √ √ "
−ζωn t 1−ζ 2 ωn t −i 1−ζ 2 ωn t
x(t) = e Ae i
+ Be , A, B ∈ R (4.10)
Com esses três casos citados acima, analisamos nosso modelo para três tipos de
soluções distintas que podem ocorrer em uma única Equação Diferencial Ordinária de
Segunda Ordem. Para efeito de curiosidade, podemos fazer uma nova substituição na
√
equação e usar a Fórmula de Euler eix = cos(x) + isen(x). Seja ωf = ωn 1 − ζ 2 . Então,

x(t) = e−ζωn t Aeiωf t + Be−iωf t
= e−ζωn t [A (cos(ωf t) + isen(ωf t)) + B (cos(ωf t) − isen(ωf t))]
= e−ζωn t [(A + B) cos(ωf t) + i (A − B) sen(ωf t)]
e, por fim,
x(t) = e−ζωn t [Ccos(ωf t) + Dsen(ωf t)] (4.11)
em que C = (A + B) e D = i (A − B), em que A, B ∈ R

Além de exemplificar o Modelo Matemático por meio de uma EDO, podemos
utilizar outros tópicos dentro da Matemática, no caso os números complexos, para so-
lucionar o problema. Muitas vezes ao estudá-los, os alunos questionam-se quanto à sua
aplicação prática, e, para o problema apresentado, foi essencial o conhecimento dos números
complexos para a sua solução.
O movimento harmônico é um dos fenômenos naturais utilizado nas Engenharias
que pode ser modelado matematicamente por meio das Equações Diferenciais Ordinárias.
De fato, sugerimos autores como Halliday, Walker e Resnick (2010a) e Kraige e Merian
(1997) aos leitores interessados no aprofundamento deste estudo.
Da mesma forma como apresentamos na vibração livre não-amortecida, o processo
passou pelas etapas descritas no Capítulo 3. Para o estudo das vibrações livres amortecidas,
nossa hipótese contém o acréscimo de um termo chamado de coeficiente de atrito viscoso,
uma vez que o fenômeno apresentado contém atrito. Como abstração, a nossa outra hipótese
continua sendo a 2a Lei de Newton, porém acrescida de um novo termo, ou seja, −cx
explicado na modelagem do fenômeno e, com isso, chegamos à fórmula mx + cx + kx = 0,
etapa denominada formulação do modelo, uma equação diferencial ordinária de segunda
ordem com coeficientes constantes, uma vez que a massa, o coeficiente de atrito e a
constante de elasticidade da mola não variam com o tempo.
Tendo a equação, passamos à resolução por intermédio das técnicas de resoluções
de EDO mostradas no Capítulo 2 por meio do Teorema 2.4.1 e subdividindo-o nas três
categorias de movimento apresentadas.
Com a função solução do sistema considerado, passamos a avaliá-lo de acordo
com as três categorias de movimento. Quando o valor de Δ é estritamente positivo, a
equação característica possui duas raízes reais e distintas gerando a solução da EDO como
x(t) = A eλ1 t + B eλ2 t , sendo A e B números reais. Como temos uma combinação linear
entre duas funções exponenciais, concluímos a não oscilação do sistema por não contemplar
funções periódicas e isso vem ao encontro com a hipótese considerada do sistema possuir
uma força externa, por exemplo, o atrito.
A avaliação é considerada, também, quando Δ é nulo. Neste caso, a equação
característica possui duas raízes reais iguais gerando a solução x(t) = A eλt + Bt eλt , sendo
A e B números reais. Novamente não temos um movimento oscilatório que é caracterizado
por oscilações. Aqui temos dois termos exponenciais sendo um deles multiplicado por um
termo linear.
E, por fim, a avaliação é dada quando Δ é estritamente negativo. Neste caso,

equação característica possui duas raízes complexas, uma conjugada da outra, gerando
a solução x(t) = eαt = [A cos(βt) + B sen(βt)], sendo A e B números reais. Observe que
agora temos funções periódicas na solução, porém multiplicadas por um termo exponencial.
O movimento, neste caso, será oscilatório porém de forma decrescente uma vez que a força
de atrito está agindo no sistema. Todos os termos exponenciais aqui apresentados são
decrescentes, pois a força de atrito é contrária à força de ação que age no fenômeno. A
etapa descrita acima é denominada validação.
Nas subseções (4.1.4), (4.1.5) e (4.1.6), apresentamos três aplicações diferentes
das vibrações livres amortecidas, uma para cada tipo de movimento. Tais aplicações
nos permitirá entender cada processo e verificar a escolha de variáveis para modelar os
sistemas de acordo com as necessidades. Esta última etapa é a aplicação. Podemos citar,
por exemplo, que para projetar uma suspensão de um carro não podemos ter molas
extremamente duras, porém não podemos ter molas extremamente moles e isto significa
que a escolha da constante elástica da mola em seu projeto deverá satisfazer condições
para que o amortecedor aja para fornecer conforto e segurança para os passageiros e as
formulações de hipóteses serão condizentes para cada fenômeno a ser projetado.
As aplicações que serão apresentadas a seguir fornecem uma leitura algébrica e
gráfica de modo que possamos interpretá-la de forma distinta. Tendo essas visões, podemos
tomar decisões para até mesmo modificar o modelo caso ele não esteja de acordo com os
dados coletados e as hipóteses formuladas.
4.1.3 Aplicações
Como aplicação inicial e mais generalizada do movimento livre não-amortecido, utili-
zaremos a Transformada de Laplace para resolver o seguinte PVI por meio da equação (4.4):
⎧
⎪
⎪ x (t) + ωn2 x = 0
⎪
⎨
⎪
x(0) = x
0 m
⎪
⎪
⎩ x (0) = V m/s
0
x (t) + ωn2 x(t) = 0
L {x } + ωn2 L {x} = 0

s2 L {x} − sx(0) − x (0) + ωn2 L {x} = 0
s2 L {x} − sx0 − V0 + ωn2 L {x} = 0

s2 + ωn2 L {x} = sx0 + V0
sx0 + V0
L {x} =
s2 + ωn2
Por fim, fazendo os ajustes necessários e aplicando a transformada inversa, obtemos
V0
x(t) = x0 cos(ωn t) + sen(ωn t) (4.12)
ωn
Exemplo 4.1.1. Suponhamos que um corpo com massa de 3 kg estica uma mola em
0,05 m. No instante t = 0 s, o corpo é solto de um ponto 0,10 cm abaixo da posição de
equilíbrio com uma velocidade de 0,20 m/s. Desprezando o atrito, determina a função x(t)
que descreve o movimento livre do sistema.
Solução. Para determinar a constante elástica da mola, da Lei de Hooke temos que
F = kx =⇒ mg = kx e com os dados do exercício temos 3 · 9, 8 = k · 0, 05 =⇒ k = 588
N/m. Daí, de (4.2) temos que
mx (t) + kx(t) = 0 =⇒ 3x (t) + 588x(t) = 0
=⇒ x (t) + 196x(t) = 0
k 588
Podemos notar que ωn2 = 196 rad/s, pois ωn2 = = = 196 rad/s. Com isso,
m 3
temos o seguinte PVI:
⎧
⎪
⎪ x (t) + 196x(t) = 0
⎪
⎨
⎪
x(0) = 0, 10 cm
⎪
⎪
⎩ x (0) = V (0) = 0, 20 m/s
De (4.12), temos
V0
x(t) = x0 cos(ωn t) + sen(ωn t)
ωn
0, 20
= 0, 10cos(14t) + sen(14t)
14
= 0, 10cos(14t) + 0, 0143sen(14t)
Portanto, a função que descreve o movimento livre do sistema é x(t) = 0, 10cos(14t)+

0, 0143sen(14t).
Como o sistema não possui atrito, então a mola oscilará por tempo indefinido até
que uma força externa pare o movimento. Graficamente podemos perceber isso.
Figura 4 – Gráfico: deslocamento versus tempo.
Podemos notar um movimento com uma característica peculiar: uma única frequên-
cia angular ωn . Esse tipo de sistema só ocorre em sistemas conservativos, isto é, em que
forças externas ou dissipativas não atuam ou, se houver tais forças, elas são irrelevantes
comparadas às aplicadas ao sistema. Podemos dizer também que a quantidade de energia
inicial é igual à final formando um sistema em que a energia mecânica é constante, ou seja,
Em = Ec + Ep = constante
em que Em é a energia mecânica, Ec é a energia cinética e Ep é a energia potencial.

Ressaltamos que estamos tratando de um sistema em que a constante da mola k é
um valor fixo, mas isso acontece em uma situação ideal. Porém, nas aplicações cotidianas,
ocorre um enfraquecimento da mola após um período de tempo, isto significa que a
constante de elasticidade da mola variará (decaindo com o tempo). Então, k da equação
(4.2) é substituída por K(t) = ke−αt , com k e α positivos resultando em
mx (t) + ke−αt x(t) = 0.
Para o movimento livre amortecido, temos os três casos apresentados neste capítulo.
Exemplo 4.1.2. Para o movimento subamortecido, quando ζ < 1, utilizaremos a Trans-

formada de Laplace para resolver o seguinte PVI por meio da equação (4.7):
⎧
⎪
⎪ x (t) + 2ζωn x (t) + ωn2 x(t) = 0
⎪
⎨
⎪
x(0) = x m
0
⎪
⎪
⎩ x (0) = V m/s
0
x (t) + 2ζωn x (t) + ωn2 x = 0
L {x } + 2ζωn L {x } + ωn2 L {x} = 0

s L {x} − sx(0) − x (0) + 2ζωn sL {x} − x(0) + ωn2 L {x} = 0
2
s2 L {x} − sx0 − V0 + 2ζωn sL {x} − 2ζωn x0 + ωn2 L {x} = 0
Reordenando os termos, temos

s2 + 2ζωn s + ωn2 L {x} = sx0 + 2ζωn x0 + V0
sx0 + 2ζωn x0 + V0
L {x} =
s2 + 2ζωn s + ωn2
No denominador podemos completar quadrado da seguinte forma:
s2 + 2ζωn s + ωn2 + (ωn2 ζ 2 − ωn2 ζ 2 ) = s2 + 2ζωn s + ωn2 ζ 2 + ωn2 − ωn2 ζ 2

= (s + ζωn )2 + ωn2 (1 − ζ 2 )

ωf2
= (s + ζωn )2 + ωf2
Então,
sx0 + 2ζωn x0 + V0
L {x} =
(s + ζωn )2 + ωf2
Fazendo os ajustes necessários para poder aplicar a Transformada inversa de

Laplace, temos
# $ # $
s + ζωn ζωn x0 + V0 ωf
L {x} = x0 2
+
2
(s + ωn ) + ωf ωf (s + ωn )2 + ωf2
# $
−ζωn t ζωn x0 + V0
x(t) = e x0 cos(ωf t) + sen(ωf t) (4.13)
ωf
Exemplo 4.1.3. Para o movimento com amortecimento crítico, basta fazer ζ = 1 em
sx0 + 2ζωn x0 + V0 sx0 + 2ωn x0 + V0

L {x} = 2 2
= 2
s + 2ζωn s + ωn s + 2ωn s + ωn2
# $
x0 1
= + ωn x0 + V0
s + ωn (s + ωn )2
Aplicando a Transformada inversa de Laplace, temos
x(t) = e−ωn t [x0 + (ωn x0 + V0 )t] (4.14)
Exemplo 4.1.4. Para o movimento superamortecido, quando ζ > 1, temos duas raízes
√ √
reais distintas λ1 = ωn −ζ + ζ 2 − 1 e λ2 = ωn −ζ − ζ 2 − 1 , sabendo que ωf =
√
ωn ζ 2 − 1. Da mesma forma, temos em
sx0 + 2ζωn x0 + V0 sx0 + 2ζωn x0 + V0

L {x} = =
2
s + 2ζωn s + ωn2 (s − λ1 )(s − λ2 )
Utilizando as frações parciais, podemos reescrever L {x} da seguinte forma
A B
L {x} = +
s − λ1 s − λ 2
Encontrando A e B a partir das condições iniciais já estabelecidas e aplicando a

Transformada inversa de Laplace, temos
# $ # $
x0 λ1 + 2x0 ζωn + V0 λ1 t x0 λ2 + 2x0 ζωn + V0 λ2 t
x(t) = e − e (4.15)
2ωf 2ωf
Nesta seção, apresentamos tópicos de onde as equações diferenciais ordinárias podem

ser aplicadas na Engenharia Mecânica sob os tópicos da vibração livre-não amortecida
e amortecida. Para a primeira, mostramos uma aplicação e a resolvemos por meio da
Transformada de Laplace, sendo sua solução apresentada algébrica e graficamente como
forma de analisar a resposta que o sistema nos retorna.
Na próxima seção, apresentaremos aplicações para as vibrações livres amortecida
que, diferentemente da não amortecida, não é uma situação idealizada.
4.1.4 Suspensões mecânicas automotivas

Em meados do século VIII, as suspensões veiculares já estavam presentes com o
uso de correntes em alguns sistemas de locomoção da época, porém causavam desconforto
aos passageiros pois os veículos balançavam muito e eram muito ruidosos. Com o passar
dos séculos e o avanço tecnológico, surgiram as molas de metal que, quando usadas como
suspensão, reduzia-se muito o balanço das charretes, mas a desvantagem era que o sistema
ficava bastante robusto, aumentando consideravelmente o peso do conjunto como um todo
(KENNETT, 2000).
Por volta do século XVIII, os processos de fabricação eram mais desenvolvidos,

sendo possível a realização de projetos automotivos mais resistentes e leves tendo como
consequência novos designs de suspensões. E nos dias atuais, com a tecnologia avançando
rapidamente, as suspensões tiveram uma atenção diferenciada tanto no quesito de carros
de corridas quanto em automóveis de passeio, passando a ser desenvolvidas suspensões
eletrônicas que corrigem imperfeições na pista mais rapidamente do que uma suspensão
comum, fornecendo, assim, mais estabilidade e desempenho.
Na atualidade, o tipo de suspensão mais utilizada no Brasil é a McPherson. São
utilizadas na parte dianteira dos veículos de pequeno e médio porte que possuem tração
dianteira. Faz parte do conjunto de suspensão independente que tem por finalidade isolar
o desnivelamento da pista para apenas uma roda, ou seja, caso ela passe por um desnível,
apenas esta será deslocada, sendo a vantagem de não transferir a oscilação para todo
veículo e também não transferir o ângulo de inclinação para o eixo tendo como consequência
a não inclinação, de mesma unidade angular, do chassi.
Para a suspensão traseira, uma das mais utilizadas é a do tipo barra de torção. De
material leve e compacto, é constituída de uma barra, circular ou retangular, montada
longitudinalmente interligando as torres de amortecimento e, como o próprio nome diz, a
torção é devido ao movimento da suspensão atuando, assim, como uma mola. Segundo
Fichera, Lacagnina e Petrone (2004), a barra de torção sofre deformação elástica por meio
do movimento das rodas e das cargas aplicadas, o que afeta as características geométricas
da suspensão e deve ser modeladas como um corpo flexível.
4.1.4.1 Aplicações
Exemplo 4.1.5. O Toyota IQ, lançado em 2008, de acordo com seu catálogo de 20122 ,
tem como especificação de massa 1.215 kg, sendo suas suspensões dianteira e traseira do
tipo McPherson e barra de torção, respectivamente. De acordo com Gillespie (1992), um
bom nível de conforto em relação à taxa de amortecimento ζ está entre 0,2 e 0,4. Sendo
ζ = 0, 4, podemos modelar o afastamento da mola em relação ao equilíbrio por meio da
equação (4.7). A força peso age sobre o carro, então
F = mg =⇒ F = 1.215 · 9, 8
=⇒ F = 11.907 N
Considerando que o afastamento em relação ao equilíbrio é x = 0, 04 m, então, pela

Lei de Hooke
F = kx =⇒ 11.907 = k · 0, 04
=⇒ k = 297.675 N/m
2
Disponível em <http://www.toyota.pt/Images/iq_catalogo_2012_tcm270-1198995.pdf>. Acesso em
12/03/2014.
A frequência circular, como já vimos, é dada por ωn = k/m, então

k 297.675
ωn = = = 15, 652 rad/s
m 1.215
Em (4.13), temos o termo ωf , daí
ωf = ωn2 (1 − ζ 2 ) = 245(1 − 0, 42 )
= 14, 35 rad/s
Considerando um afastamento inicial x0 = 0, 04 m com uma velocidade inicial

V0 = 0, 2 m/s e substituindo os valores na equação (4.7), temos o seguinte PVI:
⎧
⎪
⎪ x (t) + 12, 52x (t) + 245x(t) = 0
⎪
⎨
⎪
x(0) = 0, 04 m
⎪
⎪
⎩ x (0) = V = 0, 2 m/s
0
Novamente em (4.13), temos
# $
−ζωn t ζωn x0 + V0
x(t) = e x0 cos(ωf t) + sen(ωf t)
ωf
# $
−0,4·15,652t 0, 4 · 15, 652 · 0, 04 + 0, 2
=e 0, 04 cos(14, 35t) + sen(14, 35t)
14, 35
= e−6,2608t [0, 04 cos(14, 35t) + 0, 0314 sen(14, 35t)]
Logo, a função x(t) acima descreve o movimento oscilatório subamortecido do

amortecedor do Toyota IQ. Como a função x(t) está, também, em função de seno e cosseno,
notamos o movimento oscilatório do sistema. A exponencial decrescente faz com que a
amplitude do movimento harmônico diminua à medida que o tempo cresce. Para a situação
modelada acima, podemos observar que passado um pouco mais de um segundo a mola
retorna para sua posição de equilíbrio uma vez que o atrito, neste caso, é considerado.
O sistema é projetado para manter um equilíbrio necessário para não tirar o conforto e
a segurança dos passageiros deste automóvel. Se a suspensão ficar oscilando por muito
tempo, o desgaste será maior e é bem provável que a segurança do veículo diminua.
Observe que, considerando um afastamento inicial de 0, 04 m, o movimento oscila-

tório, por meio dos outros dados considerados, apresenta um aspecto gráfico de acordo
com a Figura 5:
Figura 5 – Movimento subamortecido considerando ζ = 0, 4.
Observe que à medida que o tempo aumenta, a suspensão tem um movimento

oscilatório periódico que decresce até a posição de equilíbrio.
Na próxima subseção, apresentaremos uma aplicação para o movimento criticamente
amortecido e para isto, utilizaremos o exemplo das válvulas tipo mola com retorno rápido.
Tais molas são usadas, por exemplo, em válvulas de segurança e atuadores pneumáticos
que requerem essa função elástica da mola.
4.1.5 Válvulas de segurança e alívio de pressão

Em relação a proteção de linhas condutoras de fluidos sobrepressão, as válvulas de
segurança e alívio de pressão são consideradas os dispositivos mais importante instalados.
Para condições de projetos, tais válvulas garantem que a pressão de operação não
exceda o valor limite projetado. Isto significa que, qualquer avaria na linha condutora, por
exemplo, o excesso de pressão causado por qualquer transiente hidráulico, a válvula entrará
em ação de modo automático, eliminando o excesso de fluido causador do aumento de
pressão na linha. A não instalação deste dispositivos pode acarretar em sérias consequências
e até mesmo graves acidentes.
Estamos interessados aqui em válvulas de segurança e alívio do tipo mola que
permitem um curso rápido de elevação e retorno para a posição inicial, que consideraremos
em nossa linguagem a posição de equilíbrio. Seu funcionamento ocorre, basicamente, da
seguinte forma: suponha uma linha hidráulica projetada para uma pressão de 10 bar e caso
ocorra uma avaria qualquer na linha e essa pressão aumentar para 11 bar, por exemplo,
a válvula de alívio de pressão, automaticamente, agirá de modo que o excesso de fluido
na linha faça pressão na mola e saia por uma descarga na própria válvula e com isso a
pressão volta a ser equilibrada. Porém, para a pressão da linha não ficar abaixo da de
trabalho, a mola retorna rapidamente à sua posição de equilíbrio, mantendo, assim, as
condições iniciais de trabalho.
O conjunto é composto por uma mola na vertical, sendo que em suas extremidades
há dois componentes chamados prato de mola superior e prato de mola inferior e o local
onde o fluido “empurra” a mola chama-se disco de vedação.
4.1.5.1 Aplicações
Exemplo 4.1.6. Considere uma válvula de alívio de pressão do tipo mola montada em
uma rede hidráulica, sendo diâmetro de entrada do fluido na válvula DE = 25 mm,
diâmetro de saída DS = 50 mm, pressão de trabalho P = 10 kgf
m2
. Suponhamos que houve
uma avaria na linha e a pressão aumentou fazendo com que a mola deslocasse no mesmo
valor que o diâmetro de saída, ou seja, x0 = 50 mm, para que aliviasse a pressão e
equilibrasse novamente o sistema. Como o processo requer um retorno rápido da mola,
consideraremos ζ = 1, ou seja, o movimento é criticamente amortecido. Com esses dados,
podemos modelar matematicamente a oscilação criticamente amortecida do sistema, isto
é, determinar a função deslocamento em cada instante de tempo para o retorno da mola
em sua posição zero.
Primeiramente, temos que determinar a constante elástica k da mola. Sabendo que
“pressão é força sobre área”, então
F πd2
P = =⇒ F = P A =⇒ F = P · =⇒ F = 0, 005 kgf = 0, 05 N
A 4
e, sendo x0 = 0, 050 m, então
F = kx =⇒ 0, 05 = k · 0, 05 =⇒ k = 1 N/m
Como a mola está na vertical, então ela está sujeita a ação da gravidade e à
massa, que consideraremos a massa do prato de mola inferior sem qualquer outro tipo de
interferência. Vamos determinar a massa desse dispositivo. Assim,
F = mg =⇒ 0, 05 = m · 9, 81 =⇒ m = 0, 0051 kg
com isso, temos que

k 1
ωn = =⇒ ωn = =⇒ ωn = 14 rad/s
m 0, 0051
No momento em que a mola começar a retornar à sua posição zero, ou seja, a

válvula fechada, consideraremos a velocidade inicial V0 = 0 m/s. Disso, temos o seguinte
PVI:
⎧
⎪
⎪ x (t) + 28x (t) + 196x(t) = 0
⎪
⎨
⎪
x(0) = 0, 050 m
⎪
⎪
⎩ x (0) = V (0) = 0 m/s
cuja solução por meio da Transformada de Laplace está generalizada pela equação (4.14).
Então,
x(t) = e−ωn t [x0 + (ωn x0 + V0 )t]

= e−14 t (0, 050 + 0, 7 t)
que tem como representação gráfica:
Figura 6 – Movimento criticamente amortecido considerando ζ = 1.

Observe que à medida que o tempo aumenta, o deslocamento tende a zero, ou seja,
o dispositivo por meio de mola está voltando para a sua posição zero, isto é, a válvula
de alívio de pressão fecha-se e o sistema retorna a trabalhar com sua pressão normal.
Dissemos, no início do exemplo, que quando o movimento é criticamente amortecido há
um retorno mais rápido para a posição zero. Isso, de fato, é uma característica desse tipo
de movimento. Iremos comparar essa mesma situação se ζ > 1, ou seja, o movimento
será superamortecido que faz com que o retorno seja mais lento que o anterior que é uma
característica real desse tipo de movimento.
Exemplo 4.1.7. Considere, agora, uma válvula de alívio de pressão com as mesmas
características do exemplo anterior, porém iremos fazer uma comparação entre ambos
exemplos considerando ζ = 1, 3. A ideia aqui é observar o movimento superamortecido e
comprovar que, de fato, o retorno da mola com essa informação é mais lento. Se tal sistema
fosse projetado dessa forma, seria possível que, por conta do retorno mais lento da mola,
a pressão na linha de fluido ficasse abaixo da pressão de operação por um determinado
período de tempo. Partindo desta ideia, temos o seguinte PVI:
⎧
⎪
⎪ x (t) + 36, 4x (t) + 196x(t) = 0
⎪
⎨
⎪
x(0) = 0, 050 m
⎪
⎪
⎩ x (0) = V (0) = 0 m/s
cuja solução por meio da Transformada de Laplace está generalizada pela equação (4.15).
Então,
# $ # $
x0 λ1 + 2x0 ζωn + V0 λ1 t x0 λ2 + 2x0 ζωn + V0 λ2 t
x(t) = e − e
2ωf 2ωf
sendo λ1 = −6, 5707, λ2 = −29, 8293 e ωf = 11, 6293 rad/s. Daí, substituindo os dados,
temos:
x(t) = 0, 06413 e−6,5707 t − 0, 01413 e−29,8293 t
Observe que a função x, agora, tem duas funções exponenciais decrescentes enquanto
a função x do exemplo anterior possui uma função exponencial decrescente multiplicando
uma função linear. Esse último fator, faz com que o decréscimo seja mais rápido. Por outro
lado, a soma das duas funções exponenciais caracterizará um retorno mais suave, mais
lento do que o anterior.
Tanto uma quanto outra tenderá a zero quando t → ∞, porém com velocidades
diferentes, ou seja, uma será mais rápida que a outra em determinados intervalos de tempo.
O gráfico do movimento superamortecido deste exemplo é dado por:
Figura 7 – Movimento superamortecido considerando ζ = 1, 3.
Note que, saindo da posição inicial, a mola retornará à posição zero para um tempo
suficientemente grande, porém, como já citado de um modo mais lento. Podemos comparar
isso observando ambas as situações no mesmo plano cartesiano.
Figura 8 – Comparação entre os movimentos criticamente amortecido e superamortecido.
A função em azul representa o movimento criticamente amortecido, já a vermelha

representa o movimento superamortecido. Para t = 0, a mola está em sua posição inicial

considerada e, após o aumento do tempo, a mola retorna para a posição onde a válvula
estará fechada em sua totalidade.
dx
A velocidade pode ser determinada fazendo v = dt
. Para o primeiro exemplo,
criticamente amortecido, temos:
dx
v= = −0, 7 e−14 t + 0, 7 e−14 t − 9, 8 t e−14 t
dt
=⇒ v(t) = −9, 8 t e−14 t
e para o segundo exemplo, superamortecido, temos:

dx
v= = −0, 4214 e−6,5707 t + 0, 4214 e−29,8293 t
dt

=⇒ v(t) = −0, 4214 e−6,5707 t − e−29,8293 t
que graficamente tem o seguinte aspecto:
Figura 9 – Comparação entre as velocidades de retorno à posição zero.
Observe que, em primeiro momento, a velocidade no movimento superamortecido é

maior e, após determinado instante de tempo, ela fica menor que a velocidade no movimento
4.2. Aplicações das EDO na Engenharia Elétrica 53
criticamente amortecido, em que podemos observar graficamente que a velocidade neste

último se aproxima de zero de forma mais rápida.
É evidente que, para condições de projeto, muitos outros fatores são considerados.
A mola tem um efeito bastante importante no sistema e, ao passar do tempo, a mesma
sofre desgastes consideráveis e precisa ser objeto de manutenção contínua e de modo
preventivo. Vale lembrar que os sistemas são projetados para atuarem em condições
normais de temperatura e pressão, porém tanto em sistemas com fluidos que trabalham
em temperaturas baixas ou muito altas é necessário um cuidado especial em relação a
avarias que fazem a pressão aumentar e tendo como consequência sérios danos ao projeto
em si.
Apresentadas algumas aplicações das vibrações livres amortecidas para a Engenharia
Mecânica, mostraremos, na próxima seção, algumas aplicações das equações diferenciais
ordinárias na Engenharia Elétrica.
4.2 Aplicações das EDO na Engenharia Elétrica

No século III a.C., os mesopotâmios conheciam a condutividade de materiais como
o cobre e o ferro e, segundo Battaglin e Barreto (2011), juntamente com os materiais
isolantes - betume e argila - construíram uma bateria nomeada Bateria de Bagdá que
gerava eletricidade. A bateria elétrica, a qual conhecemos hoje, foi inventada pelo físico
italiano Alessandro Volta (1745 - 1827), determinando, para a época, que os melhores
materiais para produzir eletricidade eram o zinco e a prata, contradizendo as afirmações
do físico e filósofo Luigi Galvani (1737 - 1798) que defendia que os metais só eram capazes
de gerar eletricidade em contato com tecido animal.
A civilização grega, por volta do século VI a.C., conhecia a magnetita assim como
os chineses e por meio de suas propriedades magnéticas construíam bússolas para auxiliar
na navegação pelo Mar Mediterrâneo. A importância do magnetismo, a partir de então,
foi explorada com mais rigor e até hoje é objeto de grandes estudos.
Em meados do século VI a.C., Thales de Mileto descobre uma resina vegetal
fóssil chamada de âmbar (em grego, elektron). Ao esfregá-la em pedaços de lã ou pele,
desempenhavam a propriedade de atrair objetos, por exemplo, a palha. Essa descoberta
deu origem à ciência conhecida hoje como eletricidade.
É também verdade que entre dois e três milênios antes de Cristo, a civilização
suméria conhecia a existência da eletricidade por meio de materiais condutores e uma das
aplicações era a deposição de prata em vasos de cobre (BATTAGLIN; BARRETO, 2011).
Já na idade média e mais precisamente na Inglaterra, Willian Gilbert (1544 - 1603)
por volta de 1600 editou sua obra De Magnete em que mostrou como construir ímas
permanentes por intermédio do tratamento térmico do ferro, descreveu e desenvolveu o

conceito de espectro de um campo magnético.
Adentrando os séculos XVIII, XIX e XX, temos novas descobertas. O físico Charles
Du Fay (1698 - 1790) distinguiu os dois tipos diferentes de carga elétrica, a positiva e a
negativa. Benjamin Franklin (1706 - 1790) estudou fenômenos provenientes da eletricidade
e atestou experimentalmente que a eletricidade na atmosfera é a causa de trovões e
relâmpagos e tal experimento é a famosa experiência da pipa (ou do papagaio).
O inglês Joseph Priestley (1733 - 1804) percebeu que a força de atração entre duas
cargas elétricas era proporcional ao produto das cargas e inversamente proporcional ao
quadrado da distância entre elas. Tais conclusões foram aceitas após Charles de Coulomb
(1736 - 1806) realizar um experimento com uma balança de torção.
O físico e matemático britânico James Clerk Maxwell (1831 - 1879) formalizou a
teoria moderna do eletromagnetismo. Os fenômenos abordados foram modelados matema-
ticamente e isso culminou nas Equações de Maxwell, um grupo de equações diferenciais
parciais que é a base do eletromagnetismo moderno. Seus estudos contribuíram de forma
ampla para a Segunda Revolução Industrial que do ponto de vista tecnológico culminou
no avanço e aprimoramento das tecnologias desenvolvidas na Primeira Revolução.
Muito do que conhecemos da evolução dos fundamentos da Engenharia Elétrica
também deve-se a cientistas como Johann Carl Friedrich Gauss (1777 - 1855) que descreveu
a Lei de Gauss como a relação entre o fluxo de um campo elétrico e as cargas elétricas
que geram esse campo e a Lei de Gauss para o Magnetismo que basicamente afirma que
os polos magnéticos são inseparáveis, ou seja, a linha de campo é contínua e fechada que
parte do polo norte e caminha para o polo sul por fora do íma e retorna ao polo norte por
dentro.
O físico inglês Michael Faraday (1791 - 1867) descreveu a variação de um campo
elétrico em relação ao tempo ou como esse campo induz um campo elétrico. Esse fato é
conhecido como Lei de Faraday. Outro importante contribuinte é o físico francês André-
Marie Ampère (1775 - 1836) que tem a Lei de Ampère em sua homenagem. Essa lei afirma
que os campos magnéticos podem ser gerados por meio de correntes elétricas e por campos
elétricos, este último com a correção de Maxwell, pois estes campos variam com o tempo.
Como já citado acerca da revolução tecnológica, o avanço nos estudos da eletricidade
culminaram em avanços em sistemas de controles automáticos. De acordo com Ogata
(2003), James Watt (1736 - 1819) foi o precursor nesses sistemas e projetou o regulador
centrífugo para controlar a velocidade de uma máquina a vapor.
Em se tratando de sistemas dinâmicos na Engenharia Elétrica, e principalmente as
conexões aos sistemas de controle, Ogata (2003) retrata sua importância nas seguintes
palavras:
Como os avanços no controle automático, na teoria e na prática, vêm

produzindo meios para otimizar o desempenho dos sistemas dinâmicos,
melhorar a produtividade, diminuir o trabalho árduo de várias rotinas de
operações manuais repetitivas, entre outros, a maioria dos engenheiros e
dos cientistas devem ter agora bons conhecimentos nessa área (OGATA,
2003, p. 1).
Isso reforça ainda mais o interesse e o avanço nas pesquisas nessa área do conheci-
mento. A tecnologia cresce no intuito de facilitar diversos processos, sejam eles de produção
ou não, automatizando o sistema de controle de uma ou mais variáveis de um sistema.
Os modelos matemáticos para um sistema dinâmico envolvem definir um conjunto
de equações que contribua para representar bem o funcionamento do sistema. Nosso
estudo está interessado em apresentar alguns modelos matemáticos aplicados à Engenharia
Elétrica descritos por equações diferenciais ordinárias e nesse sentido comparar os modelos,
por exemplo os circuitos LC e RLC, com as vibrações livres não-amortecidas e amortecidas
de sistemas encontrados na Engenharia Mecânica e com isso mostrar que muitas vezes,
modelos criados em uma área específica podem muito bem se enquadrar em uma outra
área do conhecimento que, a priori, são distintas. Essa conexão facilita o trabalho do
modelador, pois a dedicação maior de tempo estará mais focada no ajustes das variáveis
para cada área estudada.
4.2.1 Corrente elétrica

De acordo com a definição de Halliday, Walker e Resnick (2010b), corrente elétrica
é um movimento de partículas carregadas. Porém, os mesmos autores ressaltam que “nem
todas as partículas carregadas que se movem produzem corrente elétrica” (HALLIDAY;
WALKER; RESNICK, 2010b, p. 141). O exemplo citado na mesma obra é o fluxo de água
em uma mangueira. Os prótons das moléculas de água representam um movimento de
cargas positivas. Os elétrons das moléculas de água caracterizam a não existência de um
fluxo líquido de cargas por causa do movimento de suas cargas negativas compensando
o movimento das cargas positivas. A consequência imediata é que “a corrente elétrica
associada ao movimento da água no interior de uma mangueira é zero” (HALLIDAY;
WALKER; RESNICK, 2010b, p. 141).
Define-se corrente elétrica média de um circuito elétrico como

Δq q 2 − q1
i= = (4.16)
Δt t2 − t1
considerada também como a taxa de variação média entre a carga elétrica q e o tempo t
em dois pontos distintos do circuito elétrico, em que
• i é a intensidade da corrente;
• Δq é a taxa de variação da carga elétrica;
• Δt é a taxa de variação do tempo.
Se quisermos saber a intensidade da corrente elétrica em um instante de tempo t

basta, na função i(t), considerarmos a taxa de variação do tempo e da carga elétrica muito
pequena, ou seja, tendendo a zero, isto é
Δq q(t + Δt) − q(t) dq

i(t) = lim = lim =
Δt→0 Δt Δ→0 Δt dt
e, portanto,
dq
i= (4.17)
dt
que é a definição de corrente elétrica instantânea.
Podemos determinar a carga elétrica que passa pelo circuito no intervalo [0, t]. Por
separação de variáveis, resolvemos a seguinte equação diferencial
t
dq
i= =⇒ dq = idt =⇒ dq = i dt
dt 0
t
=⇒ q(t) = i · t | = i · (t − 0)
0
=⇒ q(t) = it (4.18)
Consideremos um circuito elétrico em que aplica-se a mesma diferença de potencial

nas extremidades e suponhamos materiais diferentes de condução elétrica, por exemplo,
cobre, vidro, alumínio, entre outros. As características desses materiais que caracterizam
essa diferença de potencial é chamada de resistência elétrica.
Selecionando dois pontos de um condutor e aplicando uma diferença de potencial V

e medindo a intensidade da corrente elétrica i resultante, temos que a resistência elétrica
R é dada por
V
R= (4.19)
i
isolando i, temos
V
i= (4.20)
R
porém, por definição de corrente elétrica
dq V
=
dt R
e separando as variáveis e considerante um intervalo de tempo [0, t], temos

V V t
dq = · dt =⇒ dq = dt
R R 0
V
=⇒ q(t) = ·t (4.21)
R
Em (4.20), quanto maior a resistência, menor será a intensidade da corrente. Já em

(4.21), fixados V e R, a quantidade de cargas elétricas em um instante de tempo t poderá
ser maior ou menor dependendo dos valores da diferença de potencial e da resistência e
não sendo mais analisados a partir do valor da intensidade da corrente i.
4.2.2 Força eletromotriz

A produção de energia elétrica estável depende, basicamente, de um dispositivo
tal que suas extremidades estabeleça uma diferença de potencial. Para isso, precisamos
realizar trabalho, realizado por dispositivos chamados de fonte de tensão que produz o que
chamamos de força eletromotriz E “que submete os portadores de carga a uma diferença
de potencial” (HALLIDAY; WALKER; RESNICK, 2010b, p. 167).
Há diversas fontes realmente úteis, por exemplo, a bateria que é utilizada em uma
variedade de equipamentos. As usinas hidrelétricas possuem geradores de eletricidade que
criam uma diferença de potencial nos consumidores finais.
Um circuito simples é mostrado na Figura 10 formado por uma fonte de tensão e
uma resistência R. Do terminal negativo para o positivo, representamos por meio de uma
seta a força eletromotriz E.
Figura 10 – Circuito elétrico simples.

Pela seção ss , passa uma carga dq em um intervalo dt. A movimentação da carga
dq depende da realização de trabalho dW pela fonte. Disso, definimos a força eletromotriz
E da fonte em relação ao trabalho realizado como
dW
E=
dq
isto é, a força eletromotriz é o trabalho realizado pela fonte por unidade de carga transferidas
do terminal de baixo potencial (terminal negativo) para o terminal de alto potencial
(terminal positivo). No SI, a unidade é volt (joule por Coulomb).
Apresentados os conceitos de corrente elétrica e força eletromotriz, a seguir introdu-
ziremos os circuitos RC, RL, LC e RLC. Os circuitos do tipo RC (resistor - capacitor) são
simples filtros eletrônicos que são utilizados como temporizadores de sinais. Um exemplo
de sua utilização são os controladores de temperatura de um sistema de ar condicionado.
Os circuitos do tipo RL (resistor - indutância) são utilizados, por exemplo, em válvulas
solenóides que, por ação de um campo magnético gerado por uma bobina (indutor),
abrem ou fecham para controlar a passagem de fluidos. Para os circuitos LC (indutância
- capacitância) e RLC (resistência - indutância - capacitância), suas aplicações estão
concentradas em filtros eletrônicos que atenuam características não desejadas de uma
frequência separando os sinais desejados.
4.2.3 Circuito RC
Nesta seção, vamos estudar circuitos elétricos em que as correntes variam com o
tempo. A Figura 11 mostra um circuito com um capacitor de capacitância C, resistência
R e uma fonte ideal com força eletromotriz E.
Figura 11 – Circuito RC.

A capacitância “é uma medida da quantidade de carga que precisa ser acumulada

nas placas (do capacitor) para produzir uma certa diferença de potencial entre elas”
(HALLIDAY; WALKER; RESNICK, 2010b, p. 112).
A carga q de um capacitor e sua diferença de potencial V são proporcionais, isto é,
q = CV
A partir do momento que o circuito é completado, surgem correntes elétricas

acumulando cargas q nas placas do capacitor estabelecendo, assim, uma diferença de
potencial V entre elas, ou seja,
q
V =
C
Considerando ainda a Figura 11, o carregamento do capacitor de capacitância C
é dado quando varia com o tempo a carga q, a diferença de potencial V e a corrente
elétrica i. A partir do terminal negativo, iremos aplicar a regra das malhas em que “a
soma algébrica das variações de potencial encontradas ao percorrer uma malha fechada é
sempre zero” (HALLIDAY; WALKER; RESNICK, 2010b, p. 170), temos que
q
E − iR − =0
C
mas, sabemos que
dq
i=
dt
então,
dq q dq q
E−R − = 0 =⇒ R + =E
dt C dt C
e colocando na forma padrão, dividindo ambos os membros por R, temos
dq 1 E
+ q= (4.22)
dt RC R
que é uma equação diferencial ordinária linear de primeira ordem que descreve a taxa de
variação instantânea da carga q em relação ao tempo t no capacitor. Podemos resolver
essa equação com o seguinte procedimento: seja μ(t) o fator integrante, então

1/RC dt
μ(t) = e = et/RC
Multiplicando a equação (4.22) por μ(t), temos

dq 1 t/RC E t/RC
et/RC + e q= e
dt RC R
O primeiro membro da equação é a derivada do produto da função h(t) = et/RC · q

em que q é uma função que depende da variável temporal t, daí
d t/RC E t/RC
e ·q = ·e
dt R
e integrando ambos os membros diferenciando t, temos

d t/RC E t/RC
e · q dt = e dt
dt R
et/RC · q = EC et/RC + K
e isolando a função q, segue que
q(t) = EC + Ke−t/RC (4.23)
Se fornecermos condições iniciais tais que t = 0 e q(t) = 0, então
0 = EC + K =⇒ K = −EC
logo,

q(t) = EC − EC e−t/RC =⇒ q(t) = EC 1 − e−t/RC (4.24)
Note que quando t → +∞, temos que q(t) = EC que é o valor final da carga no
capacitor.
Através da derivada de q(t), que sabemos que é a corrente, podemos determinar a
corrente de carregamento do capacitor. Com efeito,

dq 1 −t/RC
i(t) = = −EC − e
dt RC
E −t/RC
= e (4.25)
R
Para t = 0, o valor inicial da corrente é E/R e quando t → +∞, o valor da

corrente tende a zero. Isso nos leva a concluir que um capacitor em seu processo de
carregamento comporta-se inicialmente como um fio comum e após um longo período de
tempo, comporta-se como um fio interrompido.
Observe que em (4.24), podemos determinar a diferença de potencial do capacitor
em seu processo de carregamento. De fato,
q(t)
q(t) = EC 1 − e−t/RC =⇒ = E 1 − e−t/RC
C
mas,
q(t)
V (t) =
C
então,

V (t) = E 1 − e−t/RC (4.26)
Note que o comportamento quando t → +∞, a diferença de potencial do capacitor

tende para E que é o seu valor final, e em seu valor inicial, ou seja, t = 0, a diferença de
potencial é nula, isto é, o capacitor está descarregado em sua totalidade.
É possível descarregar o capacitor. Observe novamente Figura 11. Se deslocarmos
a chave S para a posição b, a força eletromotriz E não estará no circuito e o capacitor
será descarregado por meio da resistência R. Como E não faz mais parte do circuito, logo
E = 03 , então nossa equação diferencial (4.22) será
dq 1
+ q=0 (4.27)
dt RC
sendo seu processo de resolução análoga à EDO (4.22) cuja solução é dada por
q(t) = Ke−t/RC (4.28)
Para um tempo suficientemente grande, a carga do capacitor tende a zero.

Aqui também podemos determinar a corrente de descarga do capacitor. Com efeito,
dq(t) K −t/RC
i(t) = =− e (4.29)
dt RC
tendo, assim, um decrescimento exponencial e para t → +∞, a corrente no capacitor tende

a zero.
Observe o termo RC. Esse produto é chamado de constante de tempo capacitiva e
é representada pela letra grega “tau” (τ ), assim
τ = RC
e como o próprio o nome diz possui dimensão temporal que representa o tempo necessário
para que o capacitor atinja a carga ou a tensão um valor equivalente a 63% do seu valor
máximo que pode ser verificado quando t = τ na equação (4.24).
De acordo com as etapas para obtenção de um modelo descritas no Capítulo 3,
primeiramente, obtemos os dados do circuito, como a força eletromotriz E, a resistência R
e a capacitância C. Esta etapa é chamada de experimentação. Com os dados levantados,
aplicamos, como hipótese, a regra das malhas no circuito para ser testada no modelo,
sendo esta etapa chamada de abstração.
Aplicada a regra das malhas, usamos a linguagem matemática para formular nosso
modelo por meio de uma equação diferencial a partir dos dados levantados, por exemplo,
sabendo que que i = dq
dt
. Chamamos esta etapa de formulação do modelo. Tendo o modelo
formulado, passamos a resolvê-lo para determinar a função solução que determina a carga
3
Quando E = 0 a corrente no resistor cai para zero, mas não de forma instantânea. O decréscimo da
corrente pode ser calculada com essa conjectura.
de carregamento do capacitor em função do tempo que é a etapa de resolução. A aplicação

do circuito RC encontra-se na seção Seção 4.4.
Na próxima subseção, introduziremos os conceitos de indutores e indutância para
assim modelar matematicamente os circuitos elétricos do tipo RL (resistência - indutor).
4.2.4 Indutores e indutância

Com as propriedades de um capacitor, podemos produzir um campo elétrico como,
já citado, o capacitor cuja parte central são placas paralelas. Para produzir um campo
magnético, podemos utilizar um indutor, sendo um dos tipos mais simples o solenóide 4 .
Suponha um solenóide com um número de espiras N sendo utilizado como indutor
conduzindo uma corrente elétrica i que produz um fluxo magnético Φ na região central do
indutor. Define-se indutância do indutor como
NΦ
L= (4.30)
i
em que o produto N Φ é chamado de enlaçamento de fluxo magnético. E, portanto, a

indutância mede o enlaçamento desse fluxo no indutor por unidade de corrente sendo sua
unidade no SI T · m2 /A, isto é, tesla-metro quadrado por ampère que denominamos de
henry (H).
De acordo com a Lei de Faraday5 , surge no indutor uma força eletromotriz EL
induzida se a corrente que o atravessa varia, tendo como consequência a variação, nas
espiras, do fluxo magnético Φ.
Auto-indução é o nome que recebe esse processo e a força eletromotriz EL é
denominada força eletromotriz autoinduzida obedecendo, assim, a Lei de Faraday. Da
equação (4.30), temos que
N Φ = Li (4.31)
e da Lei de Faraday segue que
d(N Φ)
EL = − (4.32)
dt
Aplicando a derivada em relação a t na equação (4.31), temos
d(N Φ) di
=L
dt dt
4
Solenóide é uma bobina helicoidal formada por espiras circulares próximas que, ao passar uma corrente
elétrica i, forma um campo magnético.
5
Essa lei determina que o módulo da força E eletromotriz induzida por meio de uma bobina com N
espiras é igual à taxa de variação instantânea di fluxo magnético Φ que passa pela espira.
e substituindo em (4.32) segue que
di
EL = −L (4.33)
dt
que é a definição de força eletromotriz auto-induzida. O objeto por último definido é

importante para o estudo dos circuitos RL, pois a corrente que varia nesse tipo de circuito
faz com que exista, no indutor, uma força eletromotriz auto-induzida E. Com esses conceitos
presentes, passaremos a apresentar, na seguinte subseção, o circuito do tipo RL (resistência
- indutância).
4.2.5 Circuito RL
Como nos circuitos RC, a corrente apresenta um comportamento semelhando
quando o circuito contém uma força eletromotriz E, um indutor L e uma resistência
R. Pelo fato de o circuito agora possuir um indutor, então nele aparecerá uma força
eletromotriz autoinduzida EL opondo-se ao aumento da corrente de acordo com a Lei de
Lenz6 .
Figura 12 – Circuito RL.
Observe, na Figura 12, que agora temos duas forças eletromotrizes, sendo a da
fonte, cujo valor de E é constante e EL , no indutor, é variável com o tempo produzida pela
autoindução.
No circuito representado pela Figura 12, aplicaremos a regra das malhas no sentido
horário, ou seja, o sentido da corrente. Do ponto a para o ponto b, o potencial elétrico
varia de −iR. Com a variação da corrente i, significa a existência da força eletromotriz

autoinduzida EL = L dt di
. Pela Lei de Lenz, EL opõe-se ao aumento da corrente que é
para baixo, então do ponto b para o ponto c, o potencial varia de −L dt
di
. E, finalmente, do
ponto c para o ponto a, o potencial varia de +E que é a força eletromotriz da fonte, então,
6
Essa lei, que determina o sentido da corrente induzida em uma espira, afirma que “a corrente induzida
em uma espira tem um sentido tal que o campo magnético produzido pela corrente se opõe ao campo
magnético que induz a corrente” (HALLIDAY; WALKER; RESNICK, 2010b, p. 267).
pela regra das malhas, segue que

di di
−iR − L + E = 0 =⇒ L + iR = E (4.34)
dt dt
e, colocando na forma padrão, temos

di R E
+ i= (4.35)
dt L L
seguindo abaixo a demonstração da solução que é uma função i(t). Determinando o fator
integrante μ(t), temos
μ(t) = e R/L dt = e(R/L)t
Multiplicando a equação (4.35) por μ(t), obtemos

di R (R/L)t E
e(R/L)t + e i = e(R/L)t
dt L L
O primeiro membro da equação acima é a derivada do produto da função g(t) =

e(R/L)t · i, em que i é uma função que depende da variável temporal t, daí
d (R/L)t E (R/L)t
e ·i = e
dt L
Integrando ambos os membros da equação diferenciando t, temos

d (R/L)t E (R/L)t
e · i dt = e dt
dt L
E (R/L)t
e(R/L)t · i = e +K
R
e isolando a função i, temos

E
i(t) = + K e−(R/L)t (4.36)
R
Sendo a condição inicial t = 0 e i(t) = 0, podemos encontrar o valor da constante

K. De fato,
E E
0= + K =⇒ K = −
R R
e substituindo o valor de K na equação (4.36), obtemos uma solução particular tal que
E
i(t) = 1 − e−(R/L)t (4.37)
R
O termo R/L é chamado de constante de tempo indutiva e é representado pelo

símbolo τL (letra grega “tau” com índice L) e suas características são análogas a constante
de tempo capacitiva dos circuitos RC.
Considere a equação (4.34). As manipulações algébricas que faremos terão como

consequência a análise da energia armazenada em um campo magnético, principalmente
em forma de trabalho, e nas oscilações nos circuitos LC e RLC. Multiplicando a equação
(4.34) por i, temos
di
Li + i2 R = Ei (4.38)
dt
Sabemos que i = dq/dt o que implica Ei = E dq/dt que representa a taxa de

variação instantânea de fornecimento de energia ao circuito pela fonte. A fonte, desse
modo, realiza trabalho quando a quantidade de carga dq passa por ela num intervalo de
tempo dt.
No resistor, a energia fornecida ao circuito é dissipada como energia térmica e é
representada pela termo i2 R em (4.38). Porém, a energia que não é dissipada no resistor
é armazenada no campo magnético do indutor de acordo com a lei de conservação de
energia. Assim, o termo Li di/dt representa a energia potencial magnética UL armazenada
no campo magnético em função do tempo, ou seja, a taxa dUL /dt. Então,
dUL di
= Li =⇒ dUL = Li di
dt dt
e integrando ambos os membros, temos que

UL i 1 2
dUL = Li di =⇒ UL = Li
0 0 2
em que UL é a energia magnética armazenada pelo indutor L que passa uma corrente i e
que tem como analogia a energia armazenada por um capacitor de capacitância C que
passa uma quantidade de carga q, ou seja,
1 2
UC = q
2C
Com esses dados em mãos, apresentaremos os circuitos LC e RLC comparando

suas oscilações com as da vibração livre não-amortecida e amortecida, respectivamente.
Na etapa experimentação, obtivemos os dados deste circuito, como a força eletro-
motriz E, o indutor L e a resistência R. Porém, para os circuitos RL, aparecerá uma força
eletromotriz autoinduzida EL que está de acordo com a Lei de Lenz. Tendo esses dados,
aplicamos, como hipótese, a regra das malhas no circuito RL para testarmos no modelo,
etapa esta denominada de abstração.
E como, pela Lei de Lenz EL = L dt di
, então a formulação do modelo por meio
de uma equação diferencial foi direta, sem a necessidade de buscar outras informações
dos dados obtidos inicialmente. É importante ressaltar que, para o circuito RL, estamos
interessados em encontrar uma função que determina a variação da corrente elétrica i em
função do tempo, diferentemente do circuito RC que estamos interessados em encontrar a

função da carga de carregamento do capacitor.
A próxima etapa foi a resolução do modelo por meio das técnicas de resoluções de
equações diferenciais ordinárias e sendo a aplicação apresentada na Seção 4.4.
Apresentamos, nessa subseção, o modelo matemático que descreve a corrente elétrica
i em função do tempo para os circuitos elétricos do tipo RL. A seguir, apresentaremos o
modelo matemático que descreve as oscilações da carga elétrica q nos circuitos do tipo LC
e aproveitaremos a oportunidade para comparar essas oscilações com o modelo matemático
presente nas vibrações livre não-amortecida apresentada no Capítulo 4.
4.2.6 Circuito LC
4.2.6.1 Uma comparação com a vibração livre não-amortecida
As oscilações de um circuito LC podem ser comparadas às oscilações do movimento

livre não-amortecido apresentada na seção Seção 4.1, uma vez que não há a resistência R
no circuito. A figura abaixo mostra uma representação de um circuito LC.
Figura 13 – Circuito LC.
Mostramos que o movimento livre não-amortecido é descrito pela seguinte equação

diferencial
d2 x
+ kx = 0
dt2
porém, para modelar matematicamente a oscilação do circuito LC, usaremos como recurso
as consequências provenientes da Lei de Hooke. Dessa lei, sabemos que
F = −kx =⇒ ma = −kx
dv
=⇒ m = −kx
dt
dv dx dx
=⇒ m · + kx =0
dt dt dt
dv dx
=⇒ mv + kx =0
dt dt
⎛ ⎞
d ⎜
⎜1 1 ⎟
⎟
=⇒ ⎜ mv 2 + kx2 ⎟ = 0
⎝
dt 2 ⎠
2
U
O termo entre parênteses (U ) é a soma das energias cinética e potencial do corpo e

da mola, respectivamente, ou seja,
U = UK + UM
mas como não há atrito, a energia total U não varia ao passar do tempo, logo
dU
=0
dt
que é a relação a que chegamos.

De modo análogo, a energia total, qualquer que seja o instante de tempo, em um
circuito LC é
1 2 1 2
U = UL + UC =⇒ U = Li + q
2 2C
e como afirmamos que não há resistência no circuito, então não há transformação de

energia em energia térmica, ou seja, a energia total é constante e
dU
=0
dt
daí,
dU di q dq
= Li + =0
dt dt C dt
mas,
dq di d2 q
i= =⇒ = 2
dt dt dt
então,
d2 q q d2 q 1
Li 2
+ i = 0 =⇒ L 2
+ q=0
dt C dt C
que em sua forma padrão é dada por

d2 q 1
+ q=0 (4.39)
dt2 LC
que é a equação diferencial que modela as oscilações em relação de um circuito LC.

Ao compararmos com a equação
d2 x
m + kx = 0
dt
a indutância L corresponde à massa m e 1/C corresponde à constante elástica k da mola.

Vimos no estudo do movimento livre não-amortecido que a frequência circular é
dada por

k
ωn =
m
e fazendo as substituições de acordo com oscilação do circuito LC, temos que

1
ωn =
LC
logo, a solução geral da equação (4.39) de acordo com o Teorema 2.4.2 é dada por
⎛ ⎞ ⎛ ⎞
1 ⎠ 1 ⎠
q(t) = A cos ⎝ t + B sen ⎝ t , A, B ∈ R (4.40)
LC LC
que é a carga em função do tempo.

Para determinar a corrente deste circuito, basta derivar q em relação a t, ou seja,
⎡ ⎛ ⎞ ⎛ ⎞⎤
dq(t) 1 ⎣ 1 ⎠ 1 ⎠⎦
i(t) = = B cos ⎝ t − A sen ⎝ t
dt LC LC LC
= ωn [B cos (ωn t) − A sen (ωn t)]
Em se tratando da comparação, podemos resumi-la por meio da seguinte tabela:

Tabela 1 – Comparação entre elementos do movimento li-

vre não-amortecido e do circuito LC.
Movimento livre não-amortecido Circuito LC
k (constante elástica da mola) 1/C (inverso da capacitância)

m (massa do corpo) L (indutância)
O inverso da capacitância 1/C também é conhecido como elastância que é análogo

à propriedade mecânica da constante elástica da mola k.
Para este tipo de circuito, obtivemos informações a partir dos dados do modelo
matemático do circuito RC ao multiplicarmos a equação (4.34) por i obtendo a equação
(4.38). A partir daí, analisamos cada termo da equação sendo eles Ei = dq dt
, i2 R e Li dt
di
.
dUL
Este último representa a taxa dt (energia potencial magnética) donde obtemos uma nova
informação: UL = 12 Li2 ; energia magnética armazenada pelo indutor. Analogamente, a
1
energia armazenada pelo capacitor é UC = 2C q 2 que é uma outra nova informação. Disso,
descrevemos as etapas de experimentação e abstração.
Para a formulação do modelo, iniciamos, primeiramente, com uma comparação
com o movimento livre não-amortecido e as consequências da Lei de Hooke em termos
de energia total do sistema e que, de modo análogo, podemos estudar a energia total no
circuito LC e, tendo essas informações, transformamo-as em linguagem matemática por
meio de uma equação diferencial ordinária.
Para a resolução, utilizamos o Teorema 2.4.2 demonstrado para apresentar a
solução geral da equação que é a oscilação da carga em função do tempo. A aplicação está
apresentada na Seção 4.4.
Observa-se a analogia entre o modelo matemático para descrever as oscilações da
carga elétrica e do modelo matemático que descreve as oscilações de um sistema massa-
mola. É notável essa comparação, pois para sistemas em áreas distintas podemos usar a
mesma linguagem matemática adequando as variáveis para cada tipo de fenômeno.
Essa analogia também será feita na próxima seção para os circuitos do tipo RLC. A
comparação será feita com o modelo matemática que descreve as oscilações nas vibrações
livre amortecidas.
4.2.7 Circuito RLC

4.2.7.1 Uma comparação com a vibração livre amortecida
É chamado RLC um circuito que contém uma resistência R, uma indutância L e

uma capacitância C. Diferentemente do circuito LC, a energia total U não é mais constante
devido à existência da resistência R e diminui ao passar do tempo. Com isso, parte da
energia será dissipada como energia térmica na resistência R. A perda de energia fará com
que as oscilações diminuam de amplitude de forma contínua e isso caracteriza o sistema

com amortecimento. Abaixo segue uma figura que representa um circuito RLC.
Figura 14 – Circuito RLC.
A potência elétrica em um circuito dá-se à existência da diferença de potencial,

corrente e resistência. Em se tratando de energia, a potência está relacionada com a
velocidade com que a energia é utilizada, sendo que esta é a capacidade de realizar
trabalho. O secador de cabelo é um exemplo de conversão da energia elétrica em energia
térmica e a potência associada nessa conversão é a taxa de variação instantânea dU/dt.
Nosso interesse é analisar a queda da energia total e essa taxa é dada por
dU
= −i2 R
dt
mas, sabemos que
dU di q dq
= Li +
dt dt C dt
então
di q dq
Li + = −i2 R
dt C dt
Sendo i = dq/dt e di/dt = d2 q/dt2 e dividindo ambos os membros por i, segue que
d2 q dq 1
L 2
+R + q=0 (4.41)
dt dt C
que é a equação diferencial de segunda ordem com coeficientes constantes que modelam as
oscilações amortecidas em um circuito RLC. Para expor sua forma padrão, basta dividir
ambos os membros da equação por L. Logo,
d2 q R dq 1
+ + q=0 (4.42)
dt2 L dt LC
Para a solução geral da equação diferencial, determinemos as raízes da equação

característica. De fato,
R 1
λ2 + λ+ =0
L LC

2 2
R R 1 R R 1
cujas raízes são λ1 = − + − e λ2 = − − − .
2L 2L LC 2L 2L LC
Uma expressão mais compacta é λ1 = −α + α2 − ωn2 e λ1 = −α − α2 − ωn2 em
que
R 1
α= e ωn = √
2L LC
fator de amortecimento e frequência de ressonância (rad/s), respectivamente.

Neste caso, também faremos uma análise da solução por intermédio do fator
α2− ωn2 originando as três categorias de amortecimento: superamortecido, criticamente
amortecido e subamortecido.
1. Superamortecido: quando α2 > ωn2 . Esse fato gera duas raízes reais distintas λ1 e λ2 .
Logo, a solução geral da equação diferencial (4.42) é dada por
√ √
−α+ α2 −ωn2 t −α− α2 −ωn
2 t
q(t) = A e +Be
2. Criticamente amortecido: quando α2 = ωn2 . Neste caso, a equação característica gera

duas raízes reais iguais. Logo, a solução geral da equação (4.42) é dada por
q(t) = e−αt (A + Bt)
3. Subamortecido: quando α2 < ωn2 a equação característica fornece duas raízes comple-
xas, uma conjugada da outra. O termo α2 − ωn2 é negativo, então
α2 − ωn2 = − (ωn2 − α2 ) = i ωn2 − α2
daí,
λ1 = −α + i ωn2 − α2 e λ2 = −α − i ωn2 − α2
Logo, a solução geral da equação (4.42) é dada por
! "
−αt
q(t) = e A cos ωn2 − α2 t + B sen ωn2 − α2 t
e sendo γ = ωn2 − α2 , então
q(t) = e−αt [A cos (γt) + B sen (γt)]

Em se tratando da comparação, podemos resumi-la por meio da seguinte tabela:
Tabela 2 – Comparação entre elementos do movimento livre amortecido e

do circuito RLC.
Movimento livre amortecido Circuito RLC
k (constante elástica da mola) 1/C (inverso da capacitância)

m (massa do corpo) L (indutância)
c = 2mωn ζ (coeficiente de amortecimento viscoso) R (Resistência)
No caso do circuito RLC, o inverso da capacitância 1/C é, também, a elastância

que tem propriedade análoga da constante elástica da mola k.
Da mesma forma que o circuito LC, obtivemos as informações iniciais da energia
potencial descrita a partir da equação (4.38). A hipótese inicial é analisar a queda da energia
total dada por dUdt
= −i2 R. Outra informação que levantamos foi que dU dt
di
= Li dt + Cq dq
dt
.
Aqui temos as etapas experimentação e abstração.
Com as duas informações acima, adentramos na formulação do modelo matemático
por meio da linguagem matemática (EDO). Daí, temos o modelo que descreve as oscilações
amortecidas de um circuito RLC. A etapa da resolução acompanha as etapas do Teorema
2.4.1 e subdivide-se nas três categorias de amortecimento: superamortecido, criticamente
amortecido e subamortecido. A aplicação do modelo é apresentada na Seção 4.4.
Nesta seção, apresentamos modelos matemáticos descritos por equações diferenciais
ordinárias aplicadas a alguns tipos circuitos elétricos nos quais buscamos informações
iniciais do sistema até chegar à formulação do modelo e sua solução. Na próxima seção,
apresentaremos brevemente o estudo das ressonâncias que está presente tanto na Engenha-
ria Mecânica quanto na Engenharia Elétrica. Para este trabalho, iremos trazer o exemplo
da ressonância nas pontes que podem entrar em colapso por causa deste fenômeno. Tal
estudo é muito importante e gera grande preocupações para engenheiros e físicos.
4.3 Ressonância: uma preocupação para grandes os-

cilações em sistemas sem amortecimento
Grandes oscilações em um sistema constituído por um peso e uma mola, provavel-
mente ultrapassaria o limite de sua elasticidade.
A problemática consiste em situações em que o movimento forçado está livre de
amortecimento e com uma força periódica externa agindo no sistema oscilante, força esta
próxima ou igual à frequência das vibrações presentes no sistema que podem lhe causar
sérios danos e até mesma graves tragédias. Esse fenômeno é conhecido como ressonância.
4.3. Ressonância: uma preocupação para grandes oscilações em sistemas sem amortecimento 73
Quando abordamos as vibrações livres não-amortecidas, notamos que a equação

diferencial que descreve o movimento é homogênea, ou seja, livre de uma atuação de forças
externas, nesse caso, forças oscilatórias.
Suponha um sistema descrito por meio das seguintes condições: x (t) + ωn2 x(t) =
γ0 sen(ωn t), x(0) = 0 e x (0) = 0, em que γ0 é uma constante e g(t) = γ0 sen(ωn t) é a
função que descreve o movimento oscilatório por consequência da força externa.
Por meio do método dos coeficientes a determinar 7 , solucionaremos o problema.
Primeiramente, determinaremos a solução da equação homogênea associada (ou comple-
mentar) x + ωn2 = 0, que é dada por
xc (t) = c1 cos(ωn t) + c2 sen(ωn t)
Para uma solução particular, vamos supor que a mesma é dada por
xp (t) = At cos(ωn t) + Bt sen(ωn t)
sendo
xp (t) = A [cos(ωn t) − tωn sen(ωn t)] + [sen(ωn t) + tωn cos(ωn t)]
xp (t) = −2Aωn sen(ωn t) + 2Bωn cos(ωn t) − Aωn2 t cos(ωn t) − Bωn2 t sen(ωn t)
Substituindo em nossa equação diferencial a fim de encontrarmos os valores de A e

B, temos, após realizados os cálculos
−2Aωn sen(ωn t) + 2Bωn cos(ωn t) = γ0 sen(ωn t)
γ0
em que A = − e B = 0. Portanto, a solução geral é dada por
2ωn
x(t) = xc (t) + xp (t)
γ0
= c1 cos(ωn t) + c2 sen(ωn t) − tcos(ωn t) (4.43)
2ωn
Das condições iniciais fornecidas pelo PVI, mostra-se que os valores das constantes
γ0
da solução complementar são c1 = 0 e c2 = . Logo, a solução do problema de valor
2ωn
inicial é dada por
7
Ver Boyce e Diprima (2006) , capítulo 3, seção 3.6.
γ0
x(t) = [sen(ωn t) + tcos(ωn t)] (4.44)
2ωn
Uma representação gráfica é fornecida pela figura abaixo.
Figura 15 – Gráfico do deslocamento em função do tempo - ressonância.
Graficamente, podemos perceber que, à medida que o tempo aumenta, a oscilação

do sistema aumenta. Dependendo do projeto, isso pode ser um risco muito grande, pois
grandes oscilações após um período de tempo pode ocasionar danos sérios em todo o
sistema.
Duas situações que não deveriam acontecer são abordadas por Zill e Cullen (2001).
A primeira foram dois acidentes aviários que ocorreram entre o final do ano de 1959 e
começo de 1960 no Texas e Indiana, respectivamente. Segundo o relato, em ambos uma
asa separou-se do avião.
Após a perícia técnica, verificaram que, ao ultrapassar a velocidade de 640 km/h,
uma hélice e um motor iniciaram uma trepidação causando uma força externa que não
poderia ser absorvida pela caixa do motor e logo foi transferida para a asa que, por
condições necessárias de projeto já possui oscilação e, ao passar um período de tempo,
essa força coincidiu com a frequência máxima de oscilação da asa causando a ressonância,
que não deveria ocorrer no sistema, o que foi suficiente para a quebra da asa.
O segundo relato foi o conhecido colapso da ponte Tacoma Narrows, no estado de
Washington, gerado por fortes ventos. Há controversas quanto à causa que levou a ponte
4.4. Aplicações 75
cair, porém, segundo Zill e Cullen (2001), durante meio século acreditava-se que a causa
da queda era a ressonância. Os ventos de aproximadamente 70 km/h geraram movimentos
de torção e a ponte veio abaixo.
Por exemplo, a equação x (t) + ωn2 x(t) = γ0 sen(ωn t) do nosso exemplo é linear,
porém pesquisas posteriores apontam que efeitos não lineares foram fatores principais que
causaram as oscilações na estrutura da ponte Tacoma Narrows.
Esse caso desencadeou, posteriormente, em uma série de cuidados em projetos
de estruturas grandes, por exemplo, no caso das pontes. Apesar de que a ressonância
não é considerada a causa principal do colapso na ponte do estado de Washington, os
Engenheiros Civis procuram assegurar que as oscilações não gerem vibrações inaceitáveis
que venham a formar o fenômeno da ressonância.
4.4 Aplicações
Exemplo 4.4.1 (Circuito RC ). Considere um circuito RC com resistência R = 6M Ω e
um capacitor com capacitância C = 2, 5 μF ligados em série juntamente com uma bateria
com força eletromotriz E = 12V . Por meio desses dados, podemos determinar a carga em
função do tempo de carregamento do capacitor. De fato, pela equação (4.22), temos
dq 1 E dq 1
+ q= =⇒ + q = 2 · 10−6
dt RC R dt 15
e considerando as condições iniciais t = 0 e q(t) = 0, temos o seguinte PVI:
⎧
⎪
⎪ dq 1
⎪
⎪
⎨ dt + q = 2 · 10−6
15
⎪
⎪
⎪
⎪
⎩ q(0) = 0
cuja solução, de acordo com a equação (4.24), é dada por

q(t) = EC 1 − e−t/RC =⇒ q(t) = 3 · 10−5 1 − e−t/15
E, portanto, a função q é a carga de carregamento do capacitor, sendo a constante

de tempo τ = RC = 6 · 106 · 2, 5 · 10−6 = 15 s.
Graficamente, temos a função q representada da seguinte forma:
Figura 16 – Gráfico da carga de carregamento do capacitor.
Quando t → ∞, obtemos o valor final da carga no capacitor e esse valor é q = 30 μC.

É possível notar esse fato graficamente. A função exponencial é assintótica, sendo a assíntota
a reta q = 3 · 10−5 .
Podemos obter a corrente i de carregamento do capacitor. Da equação (4.25), que
é a taxa de variação instantânea da carga no capacitor, temos que
i(t) = 2 · 10−6 e−t/15
cuja representação gráfica é dada por:
Figura 17 – Gráfico da corrente de carregamento do capacitor.

Observe que para t = 0, a corrente inicial é i = 2 · 10−6 A e, ao aumentar o tempo,

a corrente elétrica diminui exponencialmente sendo que, para um tempo suficientemente
grande, a mesma tende a zerar.
É possível determinar a diferença de potencial entre as placas do capacitor. Por
meio da equação (4.26), podemos fazer tal análise. A diferença de potencial entre as placas
é dada por

V (t) = 12 1 − e−t/15
Note que para t = 0, V = 0, instante em que o capacitor está totalmente descar-

regado e quando t → ∞, a diferença de potencial tende para E e a carga do capacitor
tende para seu valor final. O gráfico abaixo representa o comportamento da diferença de
potencial V .
Figura 18 – Gráfico da diferença de potencial entre as placas do capacitor.
Note que a diferença de potencial no capacitor, à medida que o tempo aumenta,

tende para E = 12 V que é o seu valor e esse crescimento é dado exponencialmente
conforme a função V .
Suponhamos, agora, “retirar” a força eletromotriz E do circuito, isto significa ter
E = 0. Com esse dado considerado, estaremos descarregando o capacitor. Daí, nossa nova
equação diferencial é dada por
dq 1
+ q=0
dt 15
A equação (4.26) nos fornece a solução geral, porém vamos considerar como
condições iniciais t = 0 e q(0) = q0 , logo a solução do PVI é dada por
q(t) = q0 e−t/RC =⇒ q(t) = q0 e−t/15
onde q0 = CV0 . Por exemplo, considerando V0 = 4 V e a capacitância C = 2, 5μF do

exemplo considerado, então a descarga do capacitor é dada por
q(t) = 10 e−t/15
e sabendo que a corrente é a taxa de variação instantânea da carga, determinamos também

a corrente de descarga, ou seja,
2
i(t) = − e−t/15
3
Em ambos os casos, temos um termo exponencial decrescente que, de fato, condiz
com a nossa conjectura inicial de descarga no capacitor. É plausível pensar que para um
período de tempo suficientemente longo, o capacitor descarregará por completo.
O gráfico abaixo mostra a descarga do capacitor em função do tempo.
Figura 19 – Gráfico da descarga do capacitor em função do tempo.
Para t = 0, a carga no capacitor é q = 10 C e, como a função é exponencial

decrescente, à medida que o tempo aumenta, a carga no capacitor diminui e tende a zero
para t suficientemente grande.
O gráfico a seguir mostra a corrente i de descarga do capacitor em função do tempo.

Note que, para a corrente elétrica, i aumentará em função do tempo. Basta observar que a
derivada de i em relação a t é sempre positiva.
Figura 20 – Gráfico da corrente elétrica de descarga do capacitor em função do tempo.
A análise desse gráfico é análoga em relação a carga q do capacitor. Para t = 0,

i= − 23A. O sinal negativo quer nos dizer o sentido contrário da corrente no circuito. E,
novamente, para um tempo t suficientemente grande, a corrente elétrica no circuito tende
para zero e o capacitor estará totalmente descarregado.
Exemplo 4.4.2 (Circuito RL). Considere um circuito RL com resistência R = 9, 0 Ω,

indutor com indutância L = 2, 0 mH e uma fonte com força eletromotriz E = 18, 0 V .
Com esses dados, podemos determinar a corrente elétrica i no instante de tempo t. Pela
equação (4.35), temos que
di R E di
+ i = =⇒ + 4, 5 · 103 i = 9 · 103
dt L L dt
e considerando as seguintes condições iniciais t = 0 e i(t) = 0, temos o seguinte PVI:
⎧
⎪
⎪ di 3 3
⎨ dt + 4, 5 · 10 i = 9 · 10
⎪
⎪
⎪
⎪
⎪
⎪
⎩ i(0) = 0
cuja solução, de acordo com a equação (4.37), é dada por

E 3

i(t) = 1 − e−(R/L) t =⇒ i(t) = 2 1 − e−4,5·10 t
R
E, portante, a função i é a corrente elétrica no circuito, sendo a constante de tempo

L
indutiva τL = R = 0, 22 ms.
Graficamente, temos a função i representada da seguinte forma:
Figura 21 – Gráfico da corrente elétrica no circuito RL.
Observe que quando t = 0, ou seja, a chave fechada, a corrente no circuito é nula

neste instante. À medida que o tempo aumenta, a corrente no circuito tende para o valor
final i = 2, 0 A. A função exponencial é assintótica e a reta horizontal i = 2, 0 é a assíntota.
di
Podemos determinar a diferença de potencial no indutor dado que VL = L dt e
L = 2, 0 mH e também a diferença de potencial no resistor sendo que VR = iR e R = 9, 0 Ω.
Primeiramente, determinemos no indutor. Com efeito,
3

VL = (2, 0) · 9, 0 · 103 e−4,5·10 t
3
= 18 · 103 e−4,5·10 t
e a diferença de potencial no resistor é dada por

3

VR = 2 1 − e−4,5·10 t
· (9, 0)
3

= 18 1 − e−4,5·10 t
Segue abaixo a representação gráfica da diferença de potencial no indutor L:
Figura 22 – Gráfico da diferença de potencial no indutor L.
E abaixo segue a representação gráfica da diferença de potencial no resistor R:
Figura 23 – Gráfico da diferença de potencial no indutor R.
Exemplo 4.4.3 (Circuito LC ). Suponha que um capacitor com capacitância C = 2, 0 μF

é carregado por meio de uma fonte com força eletromotriz E = 36 V . Quando a fonte que
carrega o capacitor é desligada, um indutor com indutância L = 12 mH é ligado entre

os terminais do capacitor. Com esses dados, podemos determinar as oscilações da carga
elétrica q no circuito LC. Da equação (4.39) temos que
d2 q 1 d2 q
+ q = 0 =⇒ + 24 · 106 q = 0
dt2 LC dt2
e considerando as condições iniciais t = 0, q(t) = 30 μC e q (t) = i(t) = 0, temos o seguinte

PVI:
⎧ 2
⎪
⎪ dq
⎪
⎪
⎪ 2
+ 24 · 106 q = 0
⎪
⎪ dt
⎪
⎪
⎪
⎨
⎪
⎪ q(0) = 30 μC
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩ q (0) = i(0) = 0 A
Aplicando a Transformada de Laplace na equação diferencial, temos:

d2 q
L + 24 · 106 L {q} = 0
dt2
s2 L {q} − sq(0) − i(0) + 24 · 106 L {q} = 0

L {q} s2 + 24 · 106 = 30 · 10−3 s
30 · 10−3 s
L {q} = √ 2
s2 + 2 · 103 6
e fazendo os ajustes necessários, aplicamos a Transformada inversa de Laplace obtendo

√
q(t) = 30 · 10−3 cos 2 · 103 6 t
√
= 0, 03 cos 2 · 103 6 t
que é função que modela as oscilações do circuito elétrico do tipo LC. Comparando com o
movimento livre não-amortecido, podemos notar que a carga elétrica q ficará oscilando,
por causa da função cosseno presente na solução do PVI, até que uma força externa aja
no sistema e pare tal oscilação.
dq
Podemos determinar a corrente i neste circuito. Sabemos que i = dt
, então
dq √ √
= i(t) = −60 6 sen 2 · 103 6 t
dt
Segue abaixo a representação gráfica das oscilações da carga no circuito LC:
Figura 24 – Gráfico da oscilação da carga no circuito LC.
E abaixo segue a representação gráfica das oscilações da corrente no circuito LC:
Figura 25 – Gráfico da oscilação da corrente no circuito LC.
Exemplo 4.4.4 (Circuito RLC ). Considere um circuito RLC em série com um resistor
de resistência R = 2, 0 Ω, indutor de indutância L = 10 mH e capacitor de capacitância
C = 1, 5 μF . Este tipo de circuito apresenta três categorias de amortecimento e com esses
1
dados vamos analisar em qual categoria ele encontra-se. Vimos que α = R L
e ωn = √LC
e se α2 > ωn2 é superamortecido, se α2 = ωn2 é criticamente amortecido e se α2 < ωn2 é
subamortecido. De fato,
R 2 1 1 1
α= = = e ωn = √ =√
2L 20 10 LC 15
logo, α2 < ωn2 e, portanto, a oscilação é subamortecida. Considerando as condições inicias

t = 0, q(t) = 30 μC e q (t) = i(t) = 0, temos o seguinte PVI:
⎧ 2
⎪
⎪ d q 1 dq 1
⎪
⎪ + + q=0
⎪
⎪ dt2 5 dt 15
⎪
⎪
⎪
⎪
⎨
⎪
⎪ q(0) = 30 μC
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩ q (0) = i(0) = 0 A
que tem como solução, de acordo com o Teorema 2.4.1, é dada por
⎡ ⎛ ⎞ ⎛ ⎞⎤
−0,1t ⎣0, 03 cos ⎝
17 ⎠ 300 17 ⎠⎦
q(t) = e t + 0, 003 sen ⎝ t
300 17 300
e abaixo segue a representação gráfica da função q:
Figura 26 – Gráfico da oscilação subamortecida da carga no circuito RLC.
Em comparação com o movimento livre amortecido, a resistência R do circuito é

equivalente ao coeficiente de amortecimento viscoso c. Diferentemente do circuito LC, as
oscilações serão amortecidas até num determinado tempo t e, assim, a carga, conforme
mostra o gráfico, se anulará.
dq
Podemos determinar a corrente i deste circuito. Sabemos que i = dt
, então
⎡⎛ ⎞ ⎛ ⎞⎤
dq 300 17 ⎠ 17 ⎠⎦
= i(t) = −e−0,1t ⎣⎝0, 0003 + 0, 03 sen ⎝ t
dt 17 300 300
Figura 27 – Gráfico da oscilação subamortecida da corrente no circuito RLC.
e é possível, também, observar o caráter oscilatório da corrente elétrica, uma vez que a
oscilação se dá em termos da função senoidal e o amortecimento se dá pelo fator exponencial
decrescente contido na função i.
Exemplo 4.4.5 (Circuito RLC ). Considere, agora, um circuito RLC em série com
resistência R = 6 Ω, indutor com indutância L = 10 mH e capacitor com capacitância
C = 1, 5 μF . Sabemos que
R 6 3 1 1
α= = = e ωn = √ =√
2L 20 10 LC 15
logo, α2 > ωn2 e, portanto, a oscilação é superamortecida.

Considerando as condições inicias t = 0, q(t) = 30 μC e q (t) = i(t) = 0, temos o

seguinte PVI:
⎧ 2
⎪
⎪ d q 3 dq 1
⎪
⎪ + + q=0
⎪
⎪ dt2 5 dt 15
⎪
⎪
⎪
⎪
⎨
⎪
⎪ q(0) = 30 μC
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎩ q (0) = i(0) = 0 A
Para a solução deste PVI, iremos utilizar a equação (4.15) em que foi aplicada
a Transformada de Laplace para o movimento superamortecido para a aplicação na
Engenharia Mecânica e que é equivalente para esse caso bastando ajustar os termos
referentes àquele sistema para com o circuito RLC. A equação (4.15) é:
# $ # $
x0 λ1 + x0 2ζωn + V0 λ1 t x0 λ2 + x0 2ζωn + V0 λ2 t
x(t) = e − e
2ωf 2ωf
e a equação equivalente para o circuito RLC é dada por:
⎡ ⎤ ⎡ ⎤
⎢ q 0 λ 1 + q0
R R
⎢ + i0 ⎥
⎥
⎢ q0 λ 2 + q 0
⎢ + i0 ⎥
⎥ λ t
⎢
q(t) = ⎢ L ⎥ λ1 t ⎢
−⎢ L ⎥e 2
2 ⎥e 2 ⎥
⎣ 2 R ⎦ ⎣ 2 R ⎦
√ · 2 2
−1 √ · 2 2
−1
LC 4ωn L LC 4ωn L

3 1 7 3 1 7
sendo λ1 = − + e λ2 = − − as raízes da equação característica
10 10 3 10 10 3
dada por:
3 1
λ2 + λ+ =0
5 15
Substituindo os valores na função q, temos que a solução do PVI é dada por (com
valores arredondados):
q(t) = 0, 0445 e−0,1473 t − 0, 01446 e−0,4528 t
Observe que a função q não está em termos de funções trigonométricas, logo a

oscilação não formará períodos. Como está em termos de funções exponenciais decrescentes,
então à medida que o tempo aumenta, a carga no circuito tende a zero.
O fato acima citado pode ser observado na seguinte representação gráfica da função
q.
Figura 28 – Gráfico da oscilação superamortecida da carga no circuito RLC.
Podemos determinar, também, a corrente elétrica i neste circuito. Com efeito,
dq
i(t) = = −0, 00655 e−0,1473 t + 0, 00654 e−0,4528 t
dt
Figura 29 – Gráfico da oscilação superamortecida da corrente no circuito RLC.

Novamente, podemos notar a diferença da categoria de amortecimento entre o

subamortecido e o superamortecido. A corrente, conforme o gráfico acima, apresenta a
curva da corrente elétrica que tem como característica um retorno para a origem de uma
forma mais lenta e isso se dá pelo fato de conter duas exponenciais decrescentes na função
i. Para o movimento criticamente amortecido, esse retorno será muito mais rápido que é
análogo às propriedades dos osciladores bloco-mola.
Exemplo 4.4.6 (Circuito RLC ). Considere um circuito RLC em série de resistor com
√
resistência R = 43 15 Ω, indutor com indutância L = 10 mH e capacitor com capacitância
C = 1, 5 μF . Sabemos que
√ √ √
R 4 15 1 15 1 1 15
α= = · = e ωn = √ =√ =
2L 3 10 15 LC 15 15
logo, α2 = ωn2 e, portanto, a oscilação é criticamente amortecida. Considerando as condições

inicias t = 0, q(t) = 30 μC e q (t) = i(t) = 0, temos o seguinte PVI:
⎧ 2 √
⎪
⎪ d q 2 15 dq 1
⎪
⎪
⎨ dt2 + 15 dt + 15 q = 0
⎪
⎪ q(0) = 30 μC
⎪
⎪
⎩
q (0) = i(0) = 0 A
Da mesma forma que o exemplo anterior, iremos utilizar a equação em que foi
aplicada a Transformada de Laplace, equação (4.14), para o movimento criticamente
amortecido que também é equivalente para este caso bastando fazer os ajustes necessários
dos dados para o circuito RLC. A equação (4.14) é:
x(t) = e−ωn t [x0 + (ωn x0 + V0 ) t]
e a equação equivalente para o circuito RLC é dada por:
# $
√ 1
−1/ LC t
q(t) = e q0 + √ q 0 + i0 t
LC
Substituindo os dados na função q, temos que a solução do PVI é dada por:
# $
√ 0, 03
q(t) = e−1/ 15 t
0, 03 + √ t
15
Note que esta função também não se dá em termos de funções trigonométricas,

porém ela é diferente do caso superamortecido. Agora, além da função exponencial
decrescente, temos, também, um termo linear que, como nos casos dos osciladores bloco-
mola, influenciará em um retorno mais rápido ao ponto inicial que o caso superamortecido.
Figura 30 – Gráfico da oscilação criticamente amortecida da carga no circuito RLC.
Graficamente, a função q é representada da seguinte forma:
Podemos determinar, também, a corrente elétrica i neste circuito. Com efeito,
dq 0, 03 √
i(t) = = − √ t e−1/ 15 t
dt 15
Figura 31 – Gráfico da oscilação criticamente amortecida da corrente no circuito RLC.

Como já dito, a oscilação criticamente amortecida por conter um termo linear faz
com que o retorno para a corrente e a carga inicial ocorra mais rapidamente. Comparando
com o movimento livre amortecido em sistemas com molas a reação é análoga, o retorno
para a posição de equilíbrio é mais rápido do que o movimento superamortecido.
Com isso, apresentamos a aplicação das equações diferenciais ordinárias na Enge-
nharia Elétrica fazendo uma análise de acordo com as etapas para criação de um modelo
matemático descritas no Capítulo 3.
91
5 Modelagem Matemática na Fí-

sica
5.1 Físicos e matemáticos: a fundamentação do Cál-

culo Diferencial e Integral
Quando citamos as palavras Cálculo Diferencial e Integral é possível que, primeira-
mente, dois nomes surgem em nossa mente: Newton e Leibniz.
A Matemática, no século XVII, se desenvolveu muito devido às pesquisas em novas
áreas abrindo-se, assim, um vasto campo de estudo (EVES, 2004). O Cálculo, invenção
de Newton e Leibniz, tornou-se objeto de estudo e com sua divulgação transformou-se
numa poderosa ferramenta que implica em uma série de aplicações e solucionou diversos
problemas, por exemplo, tangentes a uma curva e o estudo de máximos e mínimos sendo
estes últimos objetos de publicações de Leibniz em 1684 por meio do artigo Um novo
método para máximos e mínimos e também para tangentes.
Como já citado, Newton e Leibniz são considerados os inventores do Cálculo, mas,
historicamente, este assunto já estava implícito em diversos estudos muitos séculos antes.
Vale ressaltar que a história do desenvolvimento deste assunto segue em ordem contrária
ao que vemos na maioria dos textos didáticos sobre o assunto (EVES, 2004). Normalmente,
apresenta-se o conceito de limite, após isso, a derivada de fuma função definida como
um limite e então a integração. O método de exaustão proposto por Eudoxo é um dos
problemas históricos do Cálculo que envolvia obter áreas, volumes e comprimentos de
arcos e também a divisão de grandezas.
Arquimedes foi um dos cientistas antigos a ser um precursor do Cálculo Integral.
Seu método, chamado de Método de Equilíbrio para determinar áreas e volumes, consistia,
segundo (EVES, 2004), em seccionar uma região em um número muito grande de tiras
paralelas e imaginá-las penduradas na extremidade de uma alavanca, estabelecendo o
equilíbrio com uma figura de área ou volume de centróides conhecidos1 .
Passado aos estudos de Arquimedes, o processo de desenvolvimento da integração
ficou paralisado, só recebendo novas atenções entre os séculos XV, XVI e XVII. Johann
Kepler utilizou este processo “para determinar as áreas envolvidas em sua segunda lei do
movimento planetário e os volumes de que se ocupou em seu tratado sobre a capacidade
dos barris de vinho” (EVES, 2004, p. 424).
1
Para um detalhamento deste método para obter o volume da esfera encontra-se em Eves (2004, p.
422).
92 Capítulo 5. Modelagem Matemática na Física
Estes são só alguns exemplos. Existem outros nomes e trabalhos desenvolvidos que
recorrem a esse método mesmo sem o rigor que conhecemos hoje.
Em se tratando do Cálculo Diferencial e Integral, não podemos esquecer de dois
importantes matemáticos predecessores de Newton: John Wallis (1616 - 1703) e Isaac
Barrow (1630 - 1677). O primeiro é considerado um dos mais importantes matemáticos
do seu tempo produzindo em diversos segmentos inclusive métodos para ensinar surdos
e mudos (EVES, 2004). Wallis utilizou sistematicamente as séries em que contribuiu
bastante e, assim, Newton pode estudá-la mais a fundo posteriormente. Sua dedicação e
contribuição ficou por conta da teoria da integral. Em sua obra Arithmetica infinitorum,
de acordo com Eves (2004), surge a afirmação da seguinte fórmula de integração
1 1
xm dx =
0 m+1
em que m pode ser inteiro diferente de −1 ou um número em forma de fração.

O segundo, Isaac Barrow, Professor de Newton, teve suas contribuições mais
ligadas ao processo de diferenciação e esse fato encontra-se em sua obra Lectiones opticae
et geometricae em que tal processo se dá pelo uso do triângulo diferencial ou triângulo
de Barrow. Nesta mesma obra citada, é provado o que conhecemos hoje como Teorema
Fundamental do Cálculo, pois Barrow percebeu e é isso que o teorema garante que derivar
e integrar são operações inversas uma da outra.
Muitos outros matemáticos desta época e anteriores a estes contribuíram para o
cálculo, porém faltava uma formalização e simbologia adequada para trabalhar com essa
poderosa ferramenta. Nesta questão, os primeiros a darem suporte neste sentido foram
Newton e Leibniz de forma independente um do outro que contribuíram para sistematizar
simbologia e regras de diferenciação e integração. Newton denominou o processo de
diferenciar como fluxões sendo que a notação adotada para derivada de uma função x é ẋ.
Leibniz adotou uma simbologia diferente. Para ele, diferenciar a função y em relação a
dy
variável independente x era denotado como dx , notação esta muito mais aceita na época.
O símbolo de integração que conhecemos hoje é devido a Leibniz e tinha como obje-
tivo a somatória de partes indivisíveis, ou seja, dx ou dy. E, portanto, estava estabelecido
que integrar uma função y com relação a variável independente x se dava por

y dx
Em sua obra Um novo método para máximos e mínimos e também para tangentes,
apresenta como encontrar tangentes por meio do cálculo diferencial e encontrar quadraturas
por meio do cálculo integral. Na mesma obra demonstra propriedades operatórias das
derivadas, por exemplo:
d du dv
[u · v] = v · +u·
dx dx dx
5.1. Físicos e matemáticos: a fundamentação do Cálculo Diferencial e Integral 93
du dv
! " v· −u·
d u dx dx
=
dx v v2
d n du
[u ] = nun−1
dx dx
Foi Leibniz que introduziu termos como função, coordenadas, oscilações, entre
outros termos e juntamente com os irmãos Bernoulli deu o nome que conhecemos hoje
de Cálculo Integral. As notações modernas das operações aritméticas também é devido
a Leibniz. É considerado o pai da Lógica Moderna e o precursor da matemática mais
avançada do seu tempo.
5.1.1 A era Bernoulli
A partir dos estudos de Leibniz sobre o Cálculo Diferencial e Integral, os irmãos

Jacques Bernoulli (1654 - 1705) e Jean Bernoulli (1667 - 1748) foram grandes entusiastas
em aprender sobre o tema para compartilhá-lo, uma vez que o cálculo de fluxões de Newton
era mais conhecido na Inglaterra (BOYER, 2010).
Da histórica família de cientistas, Jacques Bernoulli é considerado o primeiro a se
destacar na Matemática, viajando para outros países para interagir com outros matemáticos
e, assim, segundo Boyer (2010), esteve a par de diversos problemas importantes, por
exemplo, as equações da catenária, tratriz e isócrona e, por meio desses problemas,
percebeu o poder do Cálculo.
A equação de Bernoulli, em homenagem a Jacques, foi uma de suas contribuições
no estudo de equações diferenciais ordinárias. Como já mostrado, a substituição z = y 1−n
transforma a equação x + p(t)x = q(t)xn em uma equação diferencial linear do tipo
z + (1 − n)p(t)z = (1 − n)q(t). Juntamente com Leibniz, Jacques procurava uma solução
para a curva braquistócrona 2 e, por fim, provou que tal curva era uma cicloide.
Apresentamos, aqui, uma breve consideração em que matemáticos e físicos deram
contribuições para o Cálculo Diferencial e Integral. As investidas para criar essa ferramenta
deu-se muito pela necessidade de solucionar problemas científicos da época de cada um
destes cientistas. Para uma abordagem mais completa da história do cálculo, recomendamos
a leitura de Eves (2004) e Boyer (2010).
2
Braquistócrona é o nome que se dá à trajetória de uma partícula sujeita a um campo gravitacional
constante, sem atrito e com velocidade inicial nula que se desloca entre dois pontos no mais curto
período de tempo.
5.2 Corpo em queda livre com influência do atrito

A equação que descreve um corpo em queda livre num meio com atrito é dada por
dv
m = mg − cv 2 (5.1)
dt
em que m é a massa do corpo em kg, g é a aceleração devido à gravidade que é aproxi-

madamente 9,81 m/s2 , v é a velocidade em m/s, t é o tempo em segundos (s) e c é uma
constante de proporcionalidade conhecida como coeficiente de arrasto cuja unidade é kg/s.
Observemos o seguinte fato na equação (5.1).
dv dx dv dv dv
= · ⇒ =v .
dt
dt dx dt dx
v
Substituindo esse último termo em (5.1), então

dv dv
mv = mg − cv 2 ⇒ mv − mg = −cv 2
dx dx
e dividindo toda equação por mv e ajustando-a, temos

dv c
+ v = gv −1 (5.2)
dx m
A equação (5.2) é uma Equação Diferencial Ordinária mais conhecida como Equação
de Bernoulli. A formulação geral desse tipo de equação é dada por
dy
+ p(x)y = f (x)y n
dx
Para solucionar a equação (5.2), é necessária uma substituição a fim de recairmos
em uma EDO Linear de Primeira Ordem e, assim, encontrar o fator integrante.
Comparando essa equação com a Equação de Bernoulli temos n = −1, então
1
w = v 1−n ⇒ w = v 2 ⇒ v = w 2 . Derivando ambos os membros em relação a x, temos
dv 1 dw
= 1 ·
dx 2w 2 dx
Substituindo em (5.2), temos
1 dw c 1 1
1 · + w 2 = gw− 2
2w 2 dx m
1
Multiplicando a equação por 2w 2 , então
dw c
+ 2 w = 2g (5.3)
dx m
5.2. Corpo em queda livre com influência do atrito 95
Com isso, temos uma equação do tipo y (x) + p(x)y(x) = q(x) e podemos resolvê-la
c
encontrando o fator integrante3 μ(x). Sendo p(x) = 2 , então
m

c c
μ(x) = exp 2 dx = exp 2 x
m m
Multiplicando a equação (5.3) pelo fator integrante μ(x), temos

c dw c c c
exp 2 x · + exp 2 x · 2 w = exp 2 x · 2g
m dx m m m
Observemos que o primeiro membro da equação acima é a derivada do produto

d c
w · exp 2 x
dx m
Então,
d c c
w · exp 2 x = exp 2 x · 2g
dx m m
Integrando ambos os membros diferenciando x, temos

d c c
w · exp 2 x dx = exp 2 x · 2g dx
dx m m

c c
w · exp 2 x = 2g exp 2 x dx
m m
Como queremos a função w, solução do problema, e resolvendo a integral por substituição,
então

c mg c mg c
w · exp 2 x = · exp 2 x + K ⇒ w = + Kexp −2 x
m c m c m
Mas, no início do problema, fizemos a substituição w = v 2 , então

2 mg c
v = + Kexp −2 x
c m

mg c
v(x) = + Kexp −2 x , K∈R (5.4)
c m
Logo, a equação (5.4) é a função que representa a velocidade do corpo em queda

livre, considerando v(x) > 0.
Ao manipularmos a equação inicial, recaímos na forma da Equação de Bernoulli.
Esse é um dos muitos Modelos Matemáticos existentes que são descritos por EDO cuja
solução é resolver uma equação do tipo y + p(x)y = q(x)y n .
3
Considere, para o fator integrante, exp(x) = ex .
5.3 Aplicação
Exemplo 5.3.1. Considere um paraquedista em queda livre em que sua massa é m = 70 kg

e o coeficiente de arrasto é c = 12, 5 kg/s. Sendo a aceleração da gravidade g = 9, 81 m/s2 ,
vamos determinar a função velocidade do paraquedista em queda livre sendo que x(0) = 0 m.
Com essas informações, temos o seguinte PVI:
⎧
⎪
⎪ dv −1
⎪
⎨ dx + 0, 1786 v = 9, 81 v
⎪
⎪
⎪
⎪
⎪
⎩ v(0) = 0 m/s
Da equação (5.4), sabemos que:
v(x) = 54, 94 + Ke−0,3571x
e sendo v(0) = 0 m/s, então a solução do PVI é dada por:
v(x) = 54, 94 − 54, 94 e−0,3571x
Figura 32 – Gráfico da velocidade em função do deslocamento.

5.4. Tratamento numérico para EDO e questão do corpo em queda livre 97
Observe que, para x suficientemente grande, a velocidade do paraquedista tende

√
para o valor 54, 94 ≈ 7, 4122 m/s. Isto significa dizer que, para as condições fornecidas
na problemática, a velocidade final do paraquedista próximo à superfície da terra é
aproximadamente 7, 4122 m/s.
De acordo com as etapas para formulação de um modelo apresentadas no Capítulo

3 baseadas em Bassanezi (2002), como hipótese para testar o modelo utilizamos a equação
que descreve a velocidade de um corpo em queda livre com influência do atrito, porém
uma velocidade que depende do tempo. Essa etapa é a abstração. Como formulação do
modelo, a partir da equação (5.1), a manipulamos para obter a descrição da velocidade
do corpo dependendo do espaço percorrido e, com isso, chegamos numa EDO do tipo
Bernoulli e, com essa ferramenta, a resolução foi obtida por transformá-la numa equação
diferencial linear de primeira ordem. E como última etapa, apresentamos uma aplicação.
5.4 Tratamento numérico para EDO e questão do

corpo em queda livre
Diversos problemas encontrados em diversas áreas do conhecimento são formulados
através de equações diferenciais e alguns exemplos são mostrados neste trabalho. Como
já vimos, podemos resolver esses problemas analítica ou numericamente, porém, segundo
Franco (2006), a grande maioria das equações que descrevem situações reais não são
solucionadas por métodos analíticos o que nos remete a empregar soluções numéricas para
tais equações.
Historicamente, Isaac Newton, Gottfried Wilhelm Leibniz e Leonhard Euler foram
os pioneiros em se tratando de solução numérica para equações diferenciais. A partir dos
trabalhos, principalmente, de Euler os estudos desses métodos foram impulsionados nos
quais conhecemos hoje. Um dos trabalhos que Euler publicou consistia em um método
interativo que permitia aproximar um problema de valor inicial em um determinado ponto
conhecido hoje como Método de Euler.
Matemáticos como Cauchy e Lipschitz impulsionaram melhorias e um tratamento
rigoroso nesse estudo e a consequência direta dessa importância dedicação foram as
dedicações nas melhorias das soluções numéricas de equações diferenciais ordinárias por
Karl Heun, Carl Runge e Martin Wilhelm Kutta que apresentaram métodos bastante
precisos e de fácil utilização.
A balística é uma ciência que estuda o movimento de projéteis e na época começou
a exigir resultados através de tratamento numérico. Os problemas que surgiam eram muito
complexos e exigiam uma grande quantidade de cálculos que eram necessários através da
ajuda dessas ferramentas.
Assim como o Cálculo Diferencial e Integral e os números reais mereceram uma

atenção especial através da Análise Matemática, os métodos numéricos foi e vem recebendo
atenção especial através da Análise Numérica. Com o avanço da tecnologia diversos
matemáticos tem trabalhado de forma contínua para a melhoria dessa poderosa ferramenta
que fornece resultados para sistemas onde a solução analítica dos problemas ainda não
foram descobertas pelos estudiosos.
A solução numérica dos problemas nos garante informações precisas sobre o modelo
matemático estudado. Neste capítulo apresentaremos de forma sucinta quatro métodos
para resolver equações diferenciais de primeira ordem: método de Euler, Euler melhorado,
Heun e Runge-Kutta de quarta ordem.
5.4.1 Método de Euler

Segundo Zill e Cullen (2001) esse método, conhecido também como método das
tangentes, é considerado um dos mais simples métodos para aproximar soluções de Equações
Diferenciais Ordinárias. A expressão geral do Método de Euler consiste em:
yn+1 = yn + hf (xn , yn ) (5.5)
onde yn+1 é a aproximação para a solução exata no passo n + 1, n é o número de divisões

no intervalo [a,b] e h é chamado passo de integração. Para encontrar o valor de h, seja o
intervalo definido [a,b] no domínio da função, então:
b−a
h=
n
5.4.2 Método de Euler melhorado

Consiste em uma técnica mais precisa do que a anterior para encontrar soluções
numéricas. A expressão geral do método de Euler melhorado é dada por:

1
yn+1 = yn + h f (xn , yn ) + f (xn+1 , ȳn+1 ) , (5.6)
2
em que
ȳn+1 = yn + hf (xn , yn )
5.4.3 Método de Heun

O Método de Heun é uma das técnicas mais conhecidas para os Métodos de
Runge-Kutta de ordem três e três estágios. A expressão geral desse método é dada por:
1
yn+1 = yn + h(k1 + 3k3 ), (5.7)
4
5.5. Aplicação: o movimento de queda livre 99
em que
k1 = f (xn , yn )

1 1
k2 = f xn + h, yn + hk1

3 3
2 2
k3 = f xn + h, yn + hk2
3 3
Nesse caso note que o termo k2 não aparece explicitamente no método, porém ele
deve ser calculado a cada passo na resolução do problema.
5.4.4 Método de Runge-Kutta de quarta ordem

O Método de Runge-Kutta de quarta ordem é uma das técnicas mais precisas
para se encontrar soluções numéricas de equações diferenciais ordinárias. De acordo com
Franco (2006), a obtenção desse método consiste de um sistema de onze equações e treze
incógnitas, porém a expressão abaixo é uma das infinitas soluções do sistema:
1
yn+1 = yn + h(k1 + 2k2 + 2k3 + k4 ) (5.8)
6
sendo:
k1 = f (xn , yn )

1 1
k2 = f xn + h, yn + hk1

2 2
1 1
k3 = f xn + h, yn + hk2
2 2
k4 = f (xn + h, yn + hk3 )
onde k1 , k2 , k3 e k4 são os estágios obtidos através do desenvolvimento da Série de Taylor

em torno do ponto (xn , yn ).
5.5 Aplicação: o movimento de queda livre

A aplicação que apresentaremos será baseada na modelação matemática de um corpo
em queda livre. O problema consiste em apresentar uma função que melhor representa a
velocidade de um corpo em queda livre e assim estimar sua velocidade próxima à superfície
da Terra. De acordo com os objetivos deste trabalho, tal modelo é representado por uma
equação diferencial ordinária. De acordo com a Segunda Lei de Newton, temos
F = ma (5.9)
em que F é a força resultante que age no corpo em newtons (N ), m é a massa do corpo

em quilogramas (Kg) e a é a aceleração dada em m/s2 .
Podemos reescrever a equação (5.9) da seguinte forma
F
a= (5.10)
m
Através da equação (5.10), modelaremos nosso sistema para determinar a já citada

velocidade terminal do corpo próxima à superfície da Terra. Em termos de taxa de variação
instantânea, sabemos que a aceleração de um objeto é dada pela derivada da função
velocidade em relação ao tempo, logo a equação (5.10) é dada por
dv F
= (5.11)
dt m
em que v é a velocidade (m/s) e t é o tempo (s). Caso a força resultante que atua sobre o
corpo for positiva, ele irá ter uma aceleração e caso for negativa, ele terá uma desaceleração.
Para um corpo em queda livre sob a gravidade da Terra, duas forças opostas
resultam em F : a força gravitacional para baixo FD e a força de resistência do ar, oposta
a gravitacional FU . Logo a equação da resultante é
F = FD + FU (5.12)
sendo os termos que compõe a equação (5.12)
FD = mg (5.13)
FU = −cv (5.14)
em que g é a constante gravitacional e seu valor é aproximadamente 9,81 (m/s2 ) e c é

a constante de proporcionalidade conhecida como coeficiente de arrasto (kg/s) devido
à resistência do ar. Combinando as equações de (5.11) até (5.14), obtemos a seguinte
equação
dv mg − cv
= (5.15)
dt m
e simplificando a equação, temos
dv c
=g− v (5.16)
dt m
A equação (5.16) está em termos da taxa de variação instantânea dv/dt, com isso
ela é chamada de equação diferencial ordinária na variável que estamos interessados em
estudar, isto é, a velocidade do corpo.
Obteremos a solução dessa equação através do PVI (problema de valor inicial)

abaixo
⎧
⎪
⎪ dv c
⎨ dt = g − m v
⎪
⎪
⎪
⎪
⎪
⎪
⎩ v(0) = 0
cuja solução é dada por

gm
v(t) = 1 − e−(c/m)t (5.17)
c
A partir da (5.17) utilizaremos os diferentes métodos numéricos para solução de

equações diferenciais ordinárias e compararemos os resultados obtidos com os fornecidos
pela solução analítica. A partir de então, analisaremos os erros absolutos dos métodos
descritos neste capítulo. O software utilizado nas simulações numéricas foi o Matlab.
R O
problema em estudo considera para os cálculos c = 12, 5 kg/s e m = 68, 1 kg nos primeiros
quinze segundos.
5.5.1 Tabela para os métodos numéricos

Abaixo, apresentaremos a comparação entre a solução exata e o método de Euler e
o erro absoluto entre ambos.
Tabela 3 – Velocidades (m/s) - Método de Euler e solução

exata.
Tempo (s) Solução exata Método de Euler Erro Absoluto
1 8,95318 9,8 0,84682

2 16,40498 17,80117 1,39619
3 22,60717 24,33371 1,72654
4 27,76929 29,66717 1,89787
5 32,06577 34,02165 1,95589
6 35,64175 37,57685 1,93510
7 38,61807 40,47949 1,86142
8 41,09528 42,84933 1,75405
9 43,15708 44,78418 1,62709
10 44,87314 46,36389 1,49075
11 46,30142 47,65363 1,35221
12 47,49019 48,70663 1,21644
13 48,47961 49,56635 1,08674
14 49,30312 50,26827 0,96516
15 49,98852 50,84135 0,85283
Podemos observar que o Método de Euler não é muito preciso, mas consegue
representar de forma satisfatória a velocidade do corpo em queda livre. Note que se
aumentarmos o tempo de análise tal método se aproximará da solução exata, pois o erro
absoluto diminui cada vez mais para tempos suficientemente longos.
A próxima tabela apresenta uma melhoria no método de Euler. Isso poderá ser
reparado ao analisarmos o erro absoluto entre a solução exata e tal método.
Tabela 4 – Velocidades (m/s) - Método de Euler Melho-

rado e solução exata.
1 8,95318 8,90059 0,05259

2 16,40498 16,31738 0,08760
3 22,60717 22,49773 0,10943
4 27,76929 27,64778 0,12151
5 32,06577 31,93927 0,12649
6 35,64175 35,51534 0,12641
7 38,61807 38,49525 0,12282
8 41,09528 40,97838 0,11690
9 43,15708 43,04756 0,10952
10 44,87314 44,77179 0,10135
11 46,30142 46,20858 0,09284
12 47,49019 47,40584 0,08435
13 48,47961 48,40351 0,07609
14 49,30312 49,23487 0,06825
15 49,98852 49,92763 0,06089
Observe que ao analisarmos as colunas da solução exata e do método apresentado,

a velocidade do paraquedista em queda livre encontra-se na mesma casa dos inteiros
diferentemente da Tabela 3 em que a diferença se dava logo na casa dos inteiros aumentando,
assim, o erro absoluto.
À medida que o tempo aumenta, observa-se que o erro absoluto diminui gradativa-
mente. Para tempo suficientemente grande, tal erro tenderá a zero. Podemos notar que o
Método de Euler Melhorado se comporta muito melhor que o Método de Euler.
A seguir, apresentaremos uma comparação entre a solução exata e o Método de
Heun por meio da Tabela 5. Os métodos de Runge-Kutta são considerados mais precisos e
isso poderemos notar na coluna do erro absoluto considerado.
O Método de Heun é chamado, também, de Runge-Kutta de terceira ordem porque
contém três estágios.
Tabela 5 – Velocidades (m/s) - Método de Heun e solução

exata.
1 8,95318 8,95562 0,00243

2 16,40498 16,40903 0,00405
3 22,60717 22,61223 0,00506
4 27,76929 27,7749 0,00562
5 32,06577 32,07161 0,00584
6 35,64175 35,64759 0,00584
7 38,61807 38,62374 0,00567
8 41,09528 41,10067 0,00539
9 43,15708 43,16213 0,00505
10 44,87314 44,87780 0,00467
11 46,30142 46,30569 0,00427
12 47,49019 47,49407 0,00388
13 48,47961 48,48311 0,00349
14 49,30312 49,30625 0,00314
15 49,98852 49,99132 0,00279
Observe que a diferença, em sua grande maioria, se dá na segunda casa decimal

o que faz com que o erro absoluto seja menor que os métodos anteriores. Nota-se que à
medida que aumentamos a ordem do método aumenta-se a precisão do mesmo tornando-o
bastante eficaz em estudos em que a equação diferencial não possui solução analítica.
É importante ressaltar que quanto maior a ordem do método mais numerosos são
os cálculos de f nos k estágios considerados, então, como o Método de Runge-Kutta de
quarta ordem é bastante preciso, não é necessário ir além desta ordem. A Tabela 6 a seguir
apresenta a comparação entre a solução exata e o Método de Runge-Kutta de quarta
ordem, considerado um dos métodos mais precisos e mais simples e possui precisão dos
Métodos de Taylor.
Tabela 6 – Velocidades (m/s) - Método de Runge-Kutta

de 4a ordem e solução exata.
1 8,95318 8,95309 0,00008

2 16,40498 16,40483 0,00015
3 22,60717 22,60698 0,00019
4 27,76929 27,76908 0,00021
5 32,06577 32,06555 0,00022
6 35,64175 35,64153 0,00022
7 38,61807 38,61786 0,00021
8 41,09528 41,09508 0,00019
9 43,15708 43,15689 0,00019
10 44,87314 44,87297 0,00017
11 46,30142 46,30569 0,00016
12 47,49019 47,49005 0,00014
13 48,47961 48,47948 0,00013
14 49,30312 49,303 0,00012
15 49,98852 49,98842 0,00010
Nota-se que a diferença na velocidade, em sua grande maioria, encontra-se na

quarta casa decimal o que mostra a grande precisão do método. Neste caso, este método
converge para a solução exata mais rapidamente que os apresentados anteriormente.
As tabelas acima foram apresentadas na ordem de precisão dos métodos. Notemos
que o método de Euler é o menos preciso e em compensação o método de Runge-Kutta de
quarta ordem é o mais preciso sendo seus erros absolutos dados na quarta ou quinta casa
decimal. Todos os métodos apresentados aproximam-se satisfatoriamente da velocidade
calculada pela solução exata.
A velocidade terminal do corpo é dada quando fazemos t → ∞ na equação (5.17).
Logo, a equação se reduz em
mg
v(t) = (5.18)
c
e de acordo com as condições iniciais, a velocidade terminal será
68, 1 · 9, 81
v(t) = = 53, 44 m/s
12, 5
Abaixo, apresentaremos um gráfico comparativo entre a solução exata, o método

de Euler e o de Runge-Kutta de quarta ordem.
Figura 33 – Comparação das velocidades do corpo em queda livre
De fato, podemos notar por meio da Figura 33 que o método de Runge-Kutta de

quarta ordem está muito bem ajustado com a solução exata e em conjunto com a Tabela
6 por intermédio do cálculo do erro absoluto nota-se que a diferença é muito pequena.
107
6 Modelagem Matemática na Biolo-

gia
6.1 Dinâmica de crescimento de um tumor

Tumor é um aumento de tamanho em algum tecido do corpo causada por um
desequilíbrio no sistema de divisão celular causando uma multiplicação excessiva das
células. O mecanismo de divisão tornam-se inoperantes e as células passam a proliferar-se
anormalmente sem que as necessidades do organismo controlem tal proliferação. O tumor
pode ser considerado benigno (não cancerígeno) ou maligno (cancerígeno).
Segundo o Instituto Nacional do Câncer - INCA (2008), o câncer é responsável por
cerca de 13% de todas as causas de óbito no mundo e isso representa que cerca de sete
milhões de pessoas vão a óbito por causa dessa doença, segundo dados retirados da União
Internacional Contra o Câncer (em inglês International Union Against Cancer - UICC).
No Brasil, as principais causas de morte por câncer são: pulmão, próstata e estômago nos
homens e mama, pulmão e instestino nas mulheres.
De acordo com os dados do INCA (2008), o percentual tão alto deve-se à maior
exposição desses indivíduos aos fatores de risco cancerígeno, por exemplo, padrão de vida
em relação a nutrição, trabalho, hábitos nocivos à saúde tal que o consumo de tabaco e
álcool, entre outros hábitos.
Alguns pesquisadores modelaram o processo de crescimento de tumores, por exem-
plo, Murray Eden do Massachusetts Institute of Technology, Josef Smolle juntamente
com Haro Stettner do Departamento de Dermatologia da Universidade de Graz (Graz
and Institute of Mathematics) e Benjamin Gompertz, matemático inglês do século XIX
membro da Sociedade Real de Londres. Há outros modelos também como os contínuos que
envolvem equações diferenciais envolvendo condições iniciais e de contorno, modelos com
formação de padrão fractal, por exemplo, crescimento de tumores malignos no cérebro e
assim por diante.
6.1.1 Modelo de Gompertz para o crescimento de tumores

O modelo de Gompertz, segundo Domingues (2011), é dado por

dN K
= rN ln (6.1)
dt N
em que
• N (t) é a população de células tumorais no instante t;

108 Capítulo 6. Modelagem Matemática na Biologia
• t é o instante considerado para cada quantidade de população de células;
• r é a constante positiva de crescimento interno da célula;
• K é o tamanho máximo que o tumor pode atingir com os nutrientes disponíveis, ou

seja, nossa capacidade de suporte.
Para solucionar (6.1), temos que fazer uma mudança de variável. De fato,

dN N
= −rN ln
dt K
seja

N
u = ln ⇒ N = Keu
K
dN du
⇒ = Keu substituindo em (3.1)
dt dt
du du
Keu = −ruKeu ⇒ = −ru
dt dt

du du
= −rdt ⇒ = − rdt
u u

N
ln(u) = −rt + C ⇒ ln = e−rt eC
K
−rteC
N = Kee e tendo N (0) = N0 , então
rt N0
ln( )
N (t) = Ke−e K (6.2)
que é a solução do seguinte PVI:

⎧
⎪
⎪ dN K
⎪
⎪ = rN ln
⎨ dt N
⎪
⎪
⎪
⎪
⎩ N (0) = N
0
De acordo com os estudos de Friberg e Mattson (1997), a carga letal de células

tumorais está entre 1012 − 1013 células porque as populações de células não normais e não
6.1. Dinâmica de crescimento de um tumor 109
metastáticas não podem ser excedidas por esse limite suporte. Então, iremos considerar,
para os cálculos, que nossa capacidade de suporte será K = 1013 células.
Os parâmetros da solução da equação de Gompertz serão baseados no estudo de
Domingues (2011), onde r = 0, 0060, K = 1013 e N0 = 109 . Logo, a equação (6.2) com
esses valores declarados torna-se:
−0,0024 ln(10) t
N (t) = 1013 e−e
Figura 34 – Gráfico da evolução temporal de crescimento populacional de células tumorais.
É importante ressaltar que o crescimento do tumor, ao passar do tempo, está sendo

considerado sem algum tipo de tratamento contra o mesmo. Domingues (2011) realizou
um estudo computacional onde é inserido um fator de tratamento baseado em endostatina
e percebe-se que tal fator, teoricamente, impede o crescimento acelerado das células, mas
para um determinado tempo t suficientemente grande, mesmo com o tratamento, alcançará
a capacidade suporte.
6.1.2 Modelo Logístico para o crescimento de tumores

O modelo de crescimento logístico, apresentado por Pierre François Verhulst em
1838, nos moldes de nosso estudo é dado por:

dN N
=r 1− N (6.3)
dt K
que complementa a equação do crescimento populacional exponencial proposto por Malthus

que produz taxas infinitas de populações com o crescimento do tempo que pode vir a
descrever bem inicialmente, mas para tempos suficientemente grandes foge da realidade
das populações reais. O crescimento de tumor pelo estudo de Malthus seria
dN
= rN (6.4)
dt
que tem como solução, por separação de variáveis, e sendo as condições iniciais N (0) = n0 ,
a função
N (t) = n0 ert (6.5)
De acordo com os dados propostos por Domingues (2011), a equação (6.5) fica da
seguinte forma:
N (t) = 1013 e0,006 t
mas, observe que para t → ∞, N (t) → ∞ que para o crescimento tumoral não é real
para tempos indefinidamente grandes. Porém, para o modelo logístico, isso não acontece e
veremos esse fato ao resolvermos a equação. Sendo N (0) = n0 , temos o seguinte PVI:
⎧
⎪
⎪ dN N
⎪
⎪
⎨ dt
=r 1− N
K
⎪
⎪
⎪
⎪
⎩ N (0) = N
0
De fato,

dN K −N K
=r N =⇒ dN = r dt
dt K N (K − N )
Por frações parciais, temos que

K 1 1
= +
N (K − N ) N K −N
daí, integrando ambos os membros da equação, temos que

1 1
+ dN = r dt
N K −N

1 1
dN + dN = rt + C
N K −N
ln N − ln(K − N ) = rt + C

N
ln = rt + C
K −N
6.1. Dinâmica de crescimento de um tumor 111
N
= ert · ek =⇒ N = (K − N )ert · eC
K −N
=⇒ N = K ert · eC − N ert · eC
=⇒ N + N ert · eC = K ert · eC
K ert · eC
=⇒ N (t) =
1 + ert · eC
e, das condições iniciais e fazendo os ajustes necessários, a solução do PVI é dada por
KN0
N (t) = (6.6)
N0 + (K − N0 ) e−rt
Com os dados já citados anteriormente, a equação (6.6) fica da seguinte forma:
1022
N (t) =
109 + (1013 − 109 ) e−0,006 t
que tem como representação gráfica
Figura 35 – Gráfico da evolução temporal de crescimento populacional de células tumorais

- Modelo logístico.
Observe que, diferente do modelo de crescimento exponencial, quanto t → ∞,

N (t) → 1013 , ou seja, para um tempo suficientemente grande, a população de células
tumorais tende para a capacidade de suporte K. Novamente, ressaltamos que o modelo
não contém um fator de tratamento que, como já dito, impede o crescimento acelerado
das células.
Note que ambos modelos, para t → ∞, a população de células tende à capacidade
suporte. Porém, graficamente, observamos que há diferenças em como as células tumorais
estão evoluindo com o tempo. No modelo logístico, a população de células cresce bem mais
lentamente do que o modelo de Gompertz. Ambos, porém, para um tempo suficientemente
grande, tendem à capacidade de carga. O gráfico abaixo apresenta ambos modelos plotados
no mesmo plano.
Figura 36 – Gráfico da evolução temporal: Modelo de Gompertz e Modelo logístico.
Observe que, logo de início, o modelo de Gompertz tem crescimento rápido, já

o modelo logístico forma uma curva chamada sigmóide, onde para tempos iniciais o
crescimento é mais lento. Outro ponto importante é que o modelo logístico, ao analisar
graficamente, demora, praticamente, duas vezes e meia a mais que o modelo de Gompertz
para chegar à capacidade suporte.
Em comparação com o estudo de Domingues (2011) para o modelo que está inserido
um fator de tratamento, o modelo logístico está mais próximo desse fato, mesmo não
considerando tal fator. De acordo com a análise do mesmo autor, a inserção do fator
6.2. Modelo matemático para absorção de drogas (medicamentos) 113
de tratamento “representaria um ganho de tempo e de qualidade de vida ao paciente”

(DOMINGUES, 2011, p. 111).
Para ambos os modelos apresentados neste capítulo, não partimos de hipóteses
iniciais e nem levantamento de dados conforme as etapas descritas por Bassanezi (2002)
apresentadas no Capítulo 3, estas denominadas experimentação e abstração. Apresentamos
a problemática já modelada.
Com os modelos descritos por equações diferenciais ordinárias, passamos a resolvê-lo
analiticamente com as técnicas para resolver equações diferenciais apresentadas no Capítulo
2, etapa esta denominada resolução. E como validação, comparamos ambos modelos com o
estudado por Domingues (2011) e, por fim, a aplicação, nos permitiu fazer reflexões acerca
da situação tanto algebricamente quanto graficamente.
Para uma situação como essa, a aplicação gera dados importantes para a tomada
de decisões, principalmente em melhores maneiras de fornecer um tratamento adequado
para o paciente ou até mesmo executando mudanças de variáveis no modelo para melhor
atendê-lo.
Na próxima seção, apresentaremos um outro modelo interessante que explica
matematicamente o funcionamento da absorção de drogas no organismo. Com este estudo
é possível entender o motivo para o qual tomamos uma medicação de tempos em tempos,
por exemplo, de oito em oito horas e assim por diante. Veremos que a taxa de queda do
efeito da medicação se dá por uma função exponencial decrescente.
6.2 Modelo matemático para absorção de drogas (me-

dicamentos)
Da mesma forma que na dinâmica de crescimento de um tumor, nossa proposta
para uma abordagem no ensino superior é mostrar e tentar diminuir as inquietações acerca
da motivação de estudar o Cálculo Diferencial e Integral em um curso de Biologia. Nesse
contexto, apresentaremos um outro modelo muito importante: a absorção de drogas no
sangue.
A área que estuda a absorção e o tempo de absorção de drogas no organismo é a
farmacocinética que tem como objetivo estudar o caminho percorrido pelo medicamento
no organismo desde sua administração até sua eliminação. Isso já nos remete a refletir no
decaimento do efeito da droga no organismo.
O processo é simples: através da via de administração o corpo (1) absorve a droga,
após isso (2) distribui pelo organismo, (3) ocorre a chamada biotransformação e por último
(4) a excreção. Estes fatores estão intimamente associados à dose fornecida que determinam
a concentração do composto nos seus locais de ação.
A farmacocinética é muito importante no estudo para determinar a posologia

adequada, pois para cada ser e faixa etária é necessário doses distintas e até mesmo
reajustes quando necessário. Segundo a Professora Ivy Alcoforado Felisberto1 esse estudo
nos fornece uma melhor compreensão da ação do medicamento, como por exemplo, a
intensidade do efeito, duração do efeito no organismo e sua toxicidade.
O funcionamento do modelo é simples. A primeira dose aplicada no instante t0
tem como concentração de droga no sangue zero e a consequência é a distribuição no
organismo. Disso ocorre a biotransformação ou metabolização aumentando a concentração
da substância no organismo. Porém, há um determinado momento em que a concentração
para de aumentar havendo, assim, um declínio do efeito. Passado esse determinando
instante de tempo, precisamos, no instante t1 tomar uma nova dose da droga e o processo
virá um ciclo. Diante desse fato é que os médicos indicam de quanto em quanto tempo
precisamos tomar certo medicamento.
Caso não tivermos essa regulamentação precisa para ingerir a medicação podemos
correr diversos perigos se não respeitarmos os instantes de tempo corretos, como por
exemplo, insuficiência renal e hepática, respostas clínicas não adequadas, entre outras
situações mais graves.
Portanto, o interesse do modelo farmacocinético é nos inferir das seguintes situação:
• estudar o tempo entre a entrada e “saída” de um medicamento no organismo;
• verificar o pico de concentração máxima da droga no sangue e a partir disso analisar

o decaimento do efeito no organismo;
• comparar as respostas dos diferentes tipos de aplicações da droga (oral, sublingual,

intravenosa, intramuscular, entre outras);
• estudar os possíveis acontecimentos caso ocorra o esquecimento de tomar a medicação

ou tomá-la antes do tempo correto;
• entre outras situações que possam ocorrer.
Apesar de existirem diversos modelos, cada um com sua área de interesse, nosso
modelo é simples baseando-se apenas no decaimento do efeito da droga no sangue em
função do tempo. Supondo que a taxa de variação da concentração é proporcional à
concentração existente na corrente sanguínea em cada instante t e sabendo que a dosagem
(da concentração) inicial seja C(0) = C0 , então,
dC
= −γC (6.7)
dt
em que,
1
Reflexão obtida em sua nota de aula denominada Farmacocinética.
• C(t) é a concentração de droga (medicamento) no sangue;
• C0 é a dosagem inicial ministrada absorvida pelo sangue;
• t é o tempo de ação da concentração;
• −γ é a constante de proporcionalidade negativa (devido ao decaimento da concen-

tração da droga no sangue).
Da mesma forma que o modelo anterior, a solução se dá por separação simples de

variáveis. Com efeito,

dC dC t
= −γC =⇒ = −γdt
dt C t0
=⇒ ln(C) = −γ(t − t0 ) + k
=⇒ C(t) = e−γ(t−t0 ) .ek mas, como C(0) = C0
=⇒ C(t) = C0 e−γ(t−t0 ) (6.8)
que é a solução do seguinte PVI:

⎧
⎪
⎪ dC
⎨ dt = −γC
⎪
⎪
⎪
⎪
⎪
⎪
⎩ C(0) = C
0
Para encontrar o valor de γ é necessário saber o tempo de meia-vida t1/2 biológico

que é o tempo em que a concentração do fármaco cai pela metade do seu valor inicial.
Então, em (6.8) C(t) = 12 C0 obtemos
ln(1/2)
t1/2 = .
γ
Como exemplo, o Fenobarbital, princípio ativo do Gardenal,

R tem como tempo de
meia-vida entre 50 − 140 minutos em uma pessoa adulta. De acordo com os estudos farma-
cológicos, entre 4 e 6 meias-vidas o medicamento quase atinge sua concentração máxima
plasmática e quanto mais curta for a meia-vida, mais rápido alcança-se a concentração
máxima.
É evidente que a cada dose aplicada do medicamento ele possui um período de
tempo de ação no organismo e sendo eliminado durante tal período. Vamos analisar agora
a ação do Fenobarbital com uma concentração inicial C0 = 0, 03mg/ml em um adulto de

acordo com a bula do medicamento e γ = −0, 007296 por intermédio de um tempo de
meia-vida de 95 minutos nos fornecendo a seguinte equação
C(t) = 0, 03 e−0,007296(t−t0 ) (6.9)
onde t0 é o tempo inicial do intervalo na primeira dose da medicação ou na dose após a

queda do efeito da mesma. Por exemplo, para o intervalo de tempo em minutos [0, 480[,
t0 = 0 e para o intervalo [480, 960[, t0 = 480 e assim por diante.
O gráfico da função C é representado pela figura abaixo:
Figura 37 – Gráfico: concentração versus tempo.
Podemos notar que a cada 480 minutos (8 horas) uma nova dose de Fenobarbital é
aplicada, ou seja, o tempo inicial de aplicação é a cada 480 minutos e a partir de então a
concentração da medicação decai em função do tempo.
Observe que para t suficientemente grande, a concentração da medicação tende a
concentração inicial, então o paciente, no tempo designado, tem que tomar uma outra
dose, criando, assim, um ciclo até o fim do tratamento estipulado pelo médico responsável.
Em se tratando das etapas descritas por Bassanezi (2002) apresentadas no Capí-
tulo 3, também não partirmos de levantamento de dados e hipóteses iniciais e sim da
problemática já modelada, ou seja, partirmos com a formulação do modelo pronta e, assim,
passamos a resolvê-lo partindo da etada resolução.
Com os resultados obtidos, passamos a verificar a validação a partir da aplicação por

comparar a solução do modelo com a bula do medicamente que indica que a concentração
plasmática máxima ocorre, em adultos, dentro de aproximadamente oito horas.
119
7 Modelagem Matemática na Quí-

mica
7.1 Função de onda

As leis da mecânica clássica estabelecidas por Isaac Newton tinham como principal
objeto de estudo o movimento dos objetos e planetas em que suas velocidades eram
consideradas “baixas”. No final do século XIX muitos experimentos mostraram que, para
partículas tão pequenas quanto o elétron, as leis da mecânica clássica não se aplicavam. A
consequência disso foi o surgimento de conceitos e equações de uma nova mecânica para a
descrição do movimento dessas partículas chamada mecânica quântica.
No contexto proposto nesta pesquisa que é apresentar modelos matemáticos descri-
tos por equações diferenciais ordinárias, a função de onda é uma aplicação Ψ (x) (estamos
considerando a função independente do tempo, ou seja, dependendo apenas da posição x
da partícula) solução da EDO de segunda ordem conhecida como Equação de Schrödinger.
7.1.1 Mecânica clássica versus Mecânica quântica

As três principais formulações da mecânica clássica são: Mecânica Newtoniana,
Mecânica Lagrangeana e a Mecânica Hamiltoniana. Esses estudos estão concentrados na
análise do movimento, variações da energia e as forças que atuam num determinado corpo.
A mecânica newtoniana pode ser analisada por meio de três etapas que descrevem suas
bases:
1. o vetor velocidade é a taxa de variação instantânea do vetor posição em relação ao

tempo, isto é,
d r
v = .
dt
2. O vetor momento linear é o produto da grandeza escalar massa pela grandeza vetorial
velocidade, isto é,
p = m v .
3. Força é a taxa de variação instantânea do momento linear em relação ao tempo, isto

é,
d p
F =
dt
120 Capítulo 7. Modelagem Matemática na Química
No item 3, é possível observar que, sabendo que a aceleração é a taxa de variação

instantânea da velocidade em relação ao tempo, então
d p d(m v )
F = =
dt dt
d v
=m = m a
dt
ou seja, a grandeza vetorial força é o produto da massa pela grandeza vetorial aceleração
que é a segunda lei de Newton.
É importante aqui ressaltar o campo de validade da Mecânica Newtoniana. As
três leis de Newton são válidas em sistemas com velocidades baixas em comparação com
a velocidade da luz e para corpos com massa grande em comparação com partículas
elementares e suas massas (NETO, 2004).
Neto (2004), em sua obra Mecânica Newtoniana, Lagrangeana e Hamiltoniana,
apresenta uma tabela referente aos campos de atuação das teorias físicas. Quando a
velocidade do corpo estudado é muito menor do que a da luz, encontramo-nos em um
domínio newtoniano para massas muito maiores que a massa do elétron e no domínio da
teoria quântica não relativística quando a massa do corpo é de mesma grandeza que a
massa do elétron. Porém, quando a velocidade do corpo é de mesma grandeza que a da luz,
encontramo-nos na teoria relativística para massas muito maiores que a massa do elétron
e na teoria quântica relativística para massas de mesma grandeza que a massa do elétron.
A Mecânica Newtoniana descreve muito bem movimentos de corpos com velocidades
bem inferiores a da velocidade da luz, porém, por meio de observações experimentais, as
leis de Newton forneciam resultados incoerentes quando se tratava de velocidades grandes.
Por volta de 1905, Einstein apresentou a Teoria da Relatividade propondo equações para
substituir as conhecidas equações da Mecânica Newtoniana que descreviam os movimentos
de partículas em velocidades grandes e que estavam em conformidade com os dados
experimentais realizados. A Mecânica Newtoniana é um caso particular da Mecânica
Relativística quando tratamos de velocidades de partículas bem abaixo da velocidade da
luz.
A Mecânica Lagrangeana é uma outra formulação da mecânica clássica. Nesta
formulação, as leis de Newton não são mais o ponto inicial, mas sim o Princípio de Hamilton
também conhecido como princípio da mínima ação. Este princípio, nesta mecânica, é
o equivalente à segunda lei de Newton e está interessado em estudar a evolução de um
sistema do instante t1 para o instante t2 .
Segundo Neto (2004), o sistema é caracterizado por uma função escalar L que
depende de n coordenadas e n velocidades generalizadas que, de forma compacta, é
7.1. Função de onda 121
representada da seguinte forma:
L = L(qi , q˙i , t) (7.1)
denominada como a lagrangeana do sistema em que qi podem ser
grandezas sem qualquer relacionamento com os sistemas de coordenadas

usuais, eles são chamados de coordenadas generalizadas. Inclusive, essas
coordenadas não precisam ser necessariamente ângulos e comprimentos.
Podem ser quantidades gerais realmente, como, por exemplo, amplitudes
numa expansão em série de Fourier (NETO, 2004, p. 256).
É chamado de espaço das configurações este espaço que é formado pelas coordenadas
generalizadas em que as derivadas de primeira ordem de tais coordenadas são denominadas
velocidades generalizadas.
Nesta formulação a ação é definida pela seguinte integral em um intervalo fechado
[t1 , t2 ], ou seja
t2
S= L(q, q̇, t) dt (7.2)
t1
A função lagrangeana a qual dissemos ser uma função escalar e, considerando uma
única partícula no sistema, ou seja, i = 1, é dada pela diferença entre a energia cinética T
e a potencial V , isto é,
1
L = T − V = m(ẋ)2 − V (x)
2
que, por conseguinte, gera, de forma generalizada, as equações de Lagrange para um
sistema com i partículas, ou seja,
d ∂L ∂L d2 xi dV (xi )
= =⇒ m 2 = − (7.3)
dt ∂ ẋi ∂xi dt dxi
em que são as equações equivalentes às de Newton para um sistema com i partículas.

A outra formulação da Mecânica Clássica é a Mecânica Hamiltoniana. Neste caso,
haverá a substituição da lagrangeana L por uma função H denominada função hamiltoniana
e, da mesma forma, terá como ponto de partida o princípio da mínima ação.
Tanto a Mecânica Lagrangeana como a Hamiltoniana trabalham com os rudimentos
do Cálculo variacional ou também conhecido como Cálculo das variações. Esta teoria
preocupa-se em determinar extremos, ou seja, máximos ou mínimos de funcionais. Um
funcional é uma aplicação de um espaço de funções ao invés de um espaço vetorial n
dimensional. Assim, o cálculo variacional é a seguinte aplicação:
E : F −→ R
em que F é um espaço de funções. A grosso modo, o funcional associa uma função de uma
determinada classe a um número real.
Um exemplo clássico de um funcional é a distância entre dois pontos num plano ao
longo de uma curva S escrita como y = y(x). A ideia é determinar o menor comprimento de
curva que passa por tais pontos, x1 e x2 . Considerando ds é um comprimento infinitesimal
de S, então, pelo teorema de Pitágoras
(ds)2 = (dx)2 + (dy)2 =⇒ ds = (dx)2 + (dy)2

-
. 2
. dy
/
=⇒ ds = 1+ dx
dx
e, portanto, o comprimento da curva considerada entre dois pontos é um funcional da

função y = y(x) dada por x2
S[y] = 1 + (y )2 dx
x1
De uma forma geral, um funcional pode ser escrito da seguinte forma:

x2
J[y] = f y(x), y (x), x dx
x1
e que no caso do comprimento da curva a função f é f = 1 + (y )2 . Logo, a problemática

consiste em encontrar uma função tal que J[f ] assuma um valor extremo, ou seja, um
máximo ou mínimo. Numa situação simples, os problemas podem ser solucionados ao
resolver a Equação de Euler-Lagrange a seguir:

∂f d ∂f
− =0 (7.4)
∂y dx ∂y
Voltando à formulação de Hamilton para a Mecânica Clássica, a substituição da

lagrangeana L se dá pela função hamiltoniana definida por H = H(q, p, t) em que q são
as coordenadas generalizadas, p é chamado de momento conjugado ou momento canônico
dado por pi = ∂∂L
q̇i
para um sistema com i partículas e t é a variável temporal. A função H
pode ser escrita da seguinte forma:

n
H(q, p, t) = pi q̇i − L(q, q̇, t) (7.5)
i=1
A mecânica de Hamilton é descrita por um conjunto de equações diferenciais de

primeira ordem dadas por:
∂H
q̇i = (7.6)
∂pi
∂H
ṗi = − (7.7)
∂qi
em que estas equações são chamadas de Equações de Hamilton. Para obtê-las é necessário
fazer o uso da Equação de Euler-Lagrange (NETO, 2004).
A Mecânica Hamiltoniana tem como vantagem a “parte formal da teoria, principal-
mente no que se refere à passagem para a mecânica quântica” (NETO, 2004, p. 331). Em
nossa apresentação das equações diferenciais ordinárias aplicadas à Química, analisamos
a Equação de Schrödinger para determinar a função de onda, equação esta dada por
ĤΨ = EΨ em que Ĥ é chamado de operador Hamiltoniano. De fato, a mecânica formulada
por Hamilton tornou-se uma ponte para a Mecânica Quântica.
Inicialmente, pensava-se que o movimento de partículas microscópicas poderia ser
expresso a partir da formulação da Mecânica Clássica em relação às leis de Newton. Erwin
Schrödinger, por intermédio de observações experimentais, evidenciou que as partículas
em sistemas microscópicos deveriam obedecer às leis do movimento ondulatório. E a
partir desse fato foi necessário formular uma mecânica que descrevesse o movimento das
partículas na escala atômica, ou seja, a teoria proposta por Schrödinger
é uma generalização, que inclui a teoria de Newton como um caso

particular no limite macroscópico, assim como a teoria da relatividade
de Einstein é uma generalização, que inclui a teoria de Newton como um
caso particular no limite de baixas velocidades (LIMA, 2009, p. 3).
A teoria de Newton é considerada um caso particular porque “experimentos cui-

dadosos [. . .] mostraram que a mecânica clássica falha ao analisar as transferências de
quantidades muito pequenas de energia e os movimentos de corpos com massa muito
pequena” (ATKINS; PAULA, 2008, p. 219). A Mecânica Clássica é determinística, ou seja,
está preocupada, em um dos seus estudos, em determinar a posição de uma partícula em
função do tempo, porém, a Mecânica Quântica é probabilística, ou seja, abandona-se a
noção de trajetória e preocupa-se estatisticamente com a probabilidade de a partícula
encontrar-se numa região do espaço. Então, a teoria quântica surgiu para descrever, tam-
bém, sistemas físicos em que as dimensões são iguais ou abaixo da escala atômica. Alguns
cientistas que trabalharam na base dessa nova mecânica são Erwin Schrödinger, Werner
Heisenberg, Louis de Broglie, Max Planck, Albert Einstein, entre outros.
Nesta seção, apresentamos uma breve explanação da formulação da Mecânica
Clássica dada por Newton, Lagrange e Hamilton e da formulação da Mecânica Quântica.
Para um aprofundamento do estudo da Mecânica Clássica, recomendamos a leitura de
Neto (2004). Nesta obra, o autor apresenta a formulação dessa mecânica do ponto de vista
comparativo e aplicativo em relação aos três cientistas mencionados.
7.1.2 A equação de Schrödinger

A equação que apresentaremos a seguir foi modelada pelo físico Erwin Schrödinger
e seu objetivo era determinar a função de onda de qualquer sistema. Estudaremos o caso
em que a equação não depende da variável temporal t, pois a equação tornaria-se uma
EDP (equação diferencial parcial) que está fora dos limites desta pesquisa, porém, como
curiosidade, apresentaremos sua forma tridimensional sem a resolução do problema em
questão. Sejam m e E, massa e energia de uma partícula, respectivamente, então a Equação
de Schrödinger é dada por
2 d2 Ψ
− + V (x)Ψ = EΨ (7.8)
2m dx2
em que V (x) é a energia potencial no ponto x, E é a soma das energias potencial e cinética
e é uma constante definida da seguinte forma
h
= = 1, 05457 × 10−34 Js
2π
tal que h é a constante de Planck cujo valor é aproximadamente 6, 626069 × 10−34 Js.
A equação (7.8) descreve sistemas unidimensionais. Para sistemas tridimensionais,
temos
2 2
− ∇ Ψ + V Ψ = EΨ
2m
em que ∇2 (lê-se “nabla dois”) é dado por
∂2 ∂2 ∂2
∇2 = + + .
∂x2 ∂y 2 ∂z 2
E de uma forma mais geral, podemos escrever a equação de Schrödinger como
ĤΨ = EΨ
em que Ĥ é o operador Hamiltoniano dado por
2 2
Ĥ = − ∇ + V.
2m
Como veremos a seguir, a solução da equação diferencial para sistemas unidimensi-

onais é complexa, ou seja, há um fator i imaginário que nos remete a concluir que não
é possível visualizar experimentalmente a posição da partícula pois seu movimento é
totalmente imprevisível e isso teve como consequência um princípio muito importante na
mecânica quântica que é o princípio da incerteza de Heisenberg: “é impossível especificar,
simultaneamente e com a precisão que se quiser, o momento e a posição de uma partícula”
(ATKINS; PAULA, 2008, p. 243).
7.1.3 Partícula na caixa unidimensional

Este problema consiste em analisar uma partícula com massa m confinada em uma
caixa entre duas paredes no intervalo [0, L]. Segundo Atkins e Paula (2008), “esse modelo
é uma idealização da energia potencial de uma molécula na fase gasosa que é livre para se
mover num recipiente unidimensional” (ATKINS; PAULA, 2008, p. 250).
Como exemplo de aplicações, esse modelo pode ser base para estudo da estrutura
eletrônica dos metais e também em termodinâmica estatística (movimento translacional
das partículas).
7.1.3.1 Soluções aceitáveis
Na região entre as paredes da caixa (no intervalo [0, L]) a energia potencial V (x) é
nula. Então, a equação de Schrödinger é escrita da seguinte forma
2 d2 Ψ
− = EΨ (ĤΨ = EΨ ) (7.9)
2m dx2
2 d2 Ψ k 2 2
em que Ĥ = − é o operador Hamiltoniano e Ek = é a soma das energias
2m dx2 2m
cinética e potencial.
Note que a equação (7.9) é uma EDO de segunda ordem homogênea e podemos
escrevê-la da seguinte forma
2 d2 Ψ
+ Ek Ψ = 0
2m dx2
e abaixo segue o desenvolvimento de sua solução:
2 d2 Ψ k2 d2 Ψ
+ Ψ = 0 =⇒ + k2Ψ = 0
2m dx2 2m dx2
cuja equação característica e sua solução é
√
λ2 + k 2 = 0 =⇒ λ = ± −k 2
e, portanto, as raízes são λ1 = ik e λ2 = −ik em que k é um número inteiro positivo. Logo,

pelo Teorema 2.4.2 a solução geral é dada por
Ψk (x) = Acos(kx) + Bsen(kx), A, B ∈ R (7.10)
e pela relação de Euler e±ix = cos(x) ± isen(x) podemos reescrever a equação (7.10)

Ψk (x) = Aeikx + Be−ikx = A cos(kx) + isen(kx) + B cos(kx) − isen(kx)
= (A + B)cos(kx) + (A − B)isen(kx)
e sendo C = A + B e D = (A − B)i, então a solução geral é dada por

Ψk = Csen(kx) + Dcos(kx), C, D ∈ R.
k 2 2
Em se tratando de uma partícula livre, qualquer valor de Ek = é uma solução
2m
aceitável. De acordo com Atkins e Paula (2008), quando a partícula não está livre, isto é,
ela está confinada numa determinada região do espaço, as funções Ψ (x) aceitáveis devem
satisfazer certas condições de contorno.
Impostas as condições de contorno à função Ψ (x), a energia da partícula em uma
caixa unidimensional é quantizada de forma que ela se torne uma função de onda aceitável.
Na próxima seção, apresentaremos a aplicação da partícula na caixa bidimensional que é
modelada matematicamente por uma equação diferencial parcial em duas variáveis, x e y,
uma independente da outra sendo possível a sua solução por separação de variáveis.
7.1.4 Partícula na caixa bidimensional

Agora, vamos analisar a partícula confinada em uma superfície retangular com
comprimento L1 na direção x e L2 na direção y. Novamente, a energia potencial V (x) é
nula em todos os pontos com exceção nas paredes a qual é infinita.
A diferença na equação é que temos que considerar duas variáveis, uma vez que
estamos em um espaço bidimensional; logo a equação de Schrödinger é dada por

2 ∂2Ψ ∂2Ψ
− + = EΨ (7.11)
2m ∂x2 ∂y 2
Temos que resolver uma Equação Diferencial Parcial. Para isso, utilizaremos o
método da separação de variáveis. Esse método consiste em dividir a equação inicial em
duas ou mais equações diferenciais ordinárias, uma para cada variável.
7.1.4.1 Resolução por separação de variáveis
No caso da caixa bidimensional, a solução é simples. Escreveremos a função de

onda como um produto de duas funções: uma dependendo apenas da variável x e a outra
da variável y.
A separação consiste em escrever a função de onda da seguinte forma:
Ψ (x, y) = χ(x)ϕ(y)
que resultará nas seguintes equações, cada uma para sua respectiva coordenada
2 d2 χ
− = Ex χ
2m dx2
2 d2 ϕ
− = Ey ϕ
2m dy 2
E = Ex χ + Ey ϕ
Como χ independe de y e ϕ independe de x torna-se viável a separação das

variáveis por meio de um produto de funções. Então, por propriedades das derivadas
parciais, podemos escrever
∂2Ψ ∂ 2 χϕ d2 ϕ
= =χ 2
∂x2 ∂x2 dx
e também
∂2Ψ ∂ 2 χϕ d2 χ
= = ϕ .
∂y 2 ∂y 2 dy 2
Substituindo ambas as equações em (7.11), temos

2 d2 ϕ d2 χ
− χ 2 +ϕ 2 = Eχϕ
2m dx dy
e, dividindo ambos os membros da equação por χϕ obtemos a seguinte expressão

1 d2 χ 1 d2 ϕ 2mE
+ =− 2 .
χ dx 2 ϕ dy 2
Notemos que o segundo membro da equação é uma constante, não variando, assim,
independente dos valores de x e y no primeiro membro da equação. Se variarmos x apenas
o termo que depende dessa variável altera-se e como o termo do segundo membro é uma
2mEx
constante podemos escrevê-la como − 2 . E como raciocínio análogo, podemos realizar

2mEy
o mesmo para y escrevendo-o como − 2 . Portanto, E = Ex + Ey . Logo, concluímos

que
1 d2 χ 2mEx
=− 2
χ dx 2
1 d2 ϕ 2mEy
=− 2
ϕ dy 2
que, como podemos observar, são duas equações diferenciais ordinárias cada uma para sua
respectiva coordenada.
Esse processo resultou em equações diferenciais ordinárias semelhantes à equação
de Schrödinger para a partícula numa caixa unidimensional. É possível fornecer o mesmo
tratamento algébrico para uma partícula em uma caixa tridimensional separando suas
respectivas variáveis e encontrando três equações diferenciais ordinárias.
A seguir, apresentaremos uma aplicação da função de onda que é o processo de
tunelamento no qual é um fenômeno que a partícula atravessa a barreira de potencial.
7.1.5 Processo de tunelamento

Tunelamento quântico ou também chamado efeito túnel ocorre quando uma onda,
em seu percurso, encontra uma barreira de potencial. A intensidade da onda não decresce
a zero no obstáculo mas sim exponencialmente após o contato.
No estudo da mecânica clássica, um elétron com uma determinada energia E

encontra uma barreira com energia P , ele reflete-se totalmente (para P > E). Porém, na
mecânica quântica existe uma certa probabilidade de o elétron atravessar a barreira. É
possível observar esse fenômeno por intermédio de experimentos e suas aplicações são
inúmeras, como por exemplo, os diodos túnel que é um semicondutor extremamente rápido
em que sua operação é na casa dos GHz, microscópios de varredura por sonda, transístores
de efeitos de campo, entre outros.
Através da Equação de Schrödinger é possivel calcular a probabilidade de tune-
lamento de uma partícula com massa m ao encontrar um obstáculo. Da equação (7.2)
k 2 2
cuja solução é Ψk = Aeikx + Be−ikx e E = ⇒ k = (2mE)1/2 descreve a onda antes
2m
da incidência na barreira finita. Para a região da barreira, no intervalo [0, L], a energia
potencial V é constante, então a equação de Schrödinger é
2 d2 Ψ
− + V Ψ = EΨ.
2m dx2
Considerando que as partículas tenham as somas das energias potencial e cinética

E menores do que a energia potencial V constante (E < V ), então podemos escrever a
equação da seguinte forma:
2 d2 Ψ
− + (V − E)Ψ = 0
2m dx2
De fato, nessas condições, V − E resulta em energia positiva, então de acordo com

o Teorema 2.4.1 (i) a solução é dada por
Ψ = Cekx + De−kx C, D ∈ R
k = [2m(V − E)]1/2
Observe que se E > V a solução torna-se complexa de acordo com o Teorema 2.4.2.
Com efeito,
2 d2 Ψ d2 Ψ 2m
− + (V − E)Ψ = 0 =⇒ + 2 (E − V )Ψ = 0
2m dx 2 dx 2
sendo a equação característica dada por

2m
λ2 + (E − V ) = 0
2

2m 2m 2m
cujas raízes são λ1 = i (E − V ) e λ 1 = −i (E − V ). Fazendo ϕ = (E − V ),
2 2 2
temos, como solução geral de acordo com o Teorema 2.4.2, a seguinte função:
Ψ (x) = A cos(ϕx) + B sen(ϕx), A, B ∈ R

7.2. Oscilador harmônico simples (quântico) 129
Há uma série de aplicações do efeito túnel e um exemplo são os diodos túneis

descoberto em 1973 pelo físico Leo Esaki. Este tipo de semicondutor funciona somente na
área de resistência negativa que, ao diminuir a tensão tem como consequência o aumento
da corrente, porém apenas quando tais tensões estão próximas de zero. Ao trabalhar com
tensões fora dessa região considerada, o diodo túnel funciona como um diodo comum.
Na próxima seção, apresentaremos o oscilador harmônico simples para a mecânica
quântica. No Capítulo 4, analisamos os osciladores em sistemas mecânicos e elétricos,
porém o estudo a seguir estudaremos esse tipo de oscilação para sistemas microscópicos.
7.2 Oscilador harmônico simples (quântico)

No estudo da Mecânica Clássica, em algum momento, estudamos o oscilador
harmônico simples. Neste trabalho, o estudamos no Capítulo 4 fazendo uma análise por
meio da Lei de Hooke. A descrição desse tipo de sistema é dado por uma energia potencial
tal que
1 2
V (x) = kx (7.12)
2
onde k é a constante elástica da mola e x é o deslocamento em relação ao equilíbrio do

sistema.
Vimos, também, que as soluções da equação diferencial mx + kx = 0, são funções
k
que oscilam à medida que o tempo varia de acordo com a frequência natural ωn = m .
Porém, é importante ressaltar que esse tipo de oscilador não está presente apenas em um
sistema massa-mola. Estudamos as oscilações dos circuitos elétricos e há, por exemplo, as
oscilações de um pêndulo e de fluidos.
Considere uma partícula de massa m que está sob a ação de um potencial V
conforme a expressão (7.12). Observe que esta expressão tem a forma de uma função
polinomial do segundo grau, ou seja, graficamente esta representa uma parábola, disso
decorre que a característica da energia potencial deste tipo oscilador também é denominada
como energia potencial parabólica (ATKINS; PAULA, 2008). Sabemos que k = mω 2 , então
1
V (x) = mω 2 x2 (7.13)
2
A relevância para a mecânica quântica está nos níveis de energia para sistemas
microscópicos, como nas vibrações das moléculas em gases e também em sólidos, em que a
solução da Equação de Schrödinger para este tipo de potencial é importante. Vale ressaltar
que esta não envolve a força que age sobre a partícula, mas sim a energia potencial sobre
todo o sistema. Então, a constante k pode ser entendida como um valor que indica o
quanto a energia potencial do sistema aumenta, sendo o valor de referência V = 0 e x = 0.
Sabendo que k = mω 2 , então a Equação de Schrödinger (7.8) independente do

tempo pode ser escrita como
2 d2 Ψ 1
− 2
+ mω 2 x2 Ψ = EΨ (7.14)
2m dx 2
onde substituímos (7.12) em (7.8).

Em se tratando da expressão (7.12), a figura abaixo representa graficamente a
energia potencial parabólica do sistema.
Figura 38 – Energia potencial V em função do deslocamento x.
Quando a partícula é afastada da posição inicial x = 0, a mesma passa a oscilar

no intervalo [−A, +A] e E é o nível de energia quando x = ±A e nestes pontos a energia
cinética é nula sendo a energia total igual a energia potencial.
Para solucionar a Equação de Schrödinger para o oscilador harmônico, temos que
fazer algumas substituições convenientes. Mas, antes disso, sabemos que
2 d2 Ψ 1 2
− + kx Ψ − EΨ = 0
2m dx2 2
e, colocando na forma padrão, temos que

# $
d2 Ψ 2mE kmx2
+ − Ψ =0
dx2 2 2
√
km km 2mE
Daí, sejam ρ = =⇒ ρ2 = 2 e λ = , então
2
d2 Ψ 2 2

+ λ − ρ x Ψ =0
dx2
√
e, assim, seja ξ = ρ x. Daí,
dΨ dΨ dξ dΨ √ dΨ
= =⇒ = ρ
dx dξ dx dx dξ
d2 Ψ d2 Ψ
=⇒ = ρ
dx2 dξ 2
sendo esta última implicação obtida ao observar as derivadas de ordem superior como
sendo aplicação de operadores diferenciais. Com essa mudança de variável, a equação de
Schrödinger é reescrita na variável ξ da seguinte maneira:
d2 Ψ (ξ)
2
+ (ρ − ξ 2 )Ψ (ξ) = 0 (7.15)
dξ
2 /2)
Na equação acima, fazemos Ψ (ξ) = e−(ξ η(ξ), então
dΨ (ξ) 2 2 dη(ξ)
= −ξe−(ξ /2) η(ξ) + e−(ξ /2)
dξ dξ
d2 Ψ (ξ) 2
−(ξ 2 /2) d η(ξ) 2 dη(ξ) 2
=⇒ 2
= e 2
− 2ξe−(ξ /2) + (ξ 2 − 1) e−(ξ /2) η(ξ)
dξ dξ dξ
que é o resultado, também, obtido por Donangelo e Capaz (2009). Agora, substituindo o
resultado acima na equação (7.15), temos a seguinte equação diferencial de segunda ordem
do oscilador harmônico:
d2 η(ξ) dη(ξ)
2
− 2ξ + (ρ − 1)η(ξ) (7.16)
dξ dξ
conhecida como Equação de Hermite. As soluções dessa equação diferencial são chamadas
de Polinômios de Hermite.
Definição 7.2.1 (Polinômios de Hermite). Os Polinômios de Hermite são definidos da

seguinte forma:
2
2 dn e−t
Hn (t) = (−1)n et (7.17)
dtn
para todo n ∈ N∗ . Temos também que,
Hn (−t) = (−1)n Hn (t) (7.18)
o que significa dizer que se n é par, Hn é par e se n é ímpar, Hn é ímpar.

A partir de então, podemos definir os cinco primeiros Polinômios de Hermite:
H0 (t) = 1
H1 (t) = 2t
H2 (t) = 4t2 − 2
H3 (t) = 8t3 − 12t
H4 (t) = 16t4 − 48t2 + 12
A equação (7.16) pode ser reescrita da seguinte maneira quando ρ = 2n + 1:

d2 Hn (ξ) dHn (ξ)
2
− 2ξ + 2nη(ξ) (7.19)
dξ dξ
As soluções da equação acima são dadas por:

k
H2k (ξ) = ci ξ 2i no caso par (7.20)
i=0

k
H2k+1 (ξ) = di ξ 2i+1 no caso ímpar (7.21)
i=0
em que ci e di são os coeficientes. Estes são dados, segundo Donangelo e Capaz (2009),
conforme a seguinte relação de recorrência:
4(i − k)
ci+1 = ci para os polinômios de grau 2k. (7.22)
2(i + 1)(2i + 1)
4(i − k)
di+1 = di para os polinômios de grau 2k + 1. (7.23)
2(i + 1)(2i + 3)
Para um aprofundamento do estudo do oscilador harmônico simples e as soluções

das equações de Hermite, sugerimos a leitura de Donangelo e Capaz (2009).
Neste capítulo, apresentamos a Equação de Schrödinger que tem como solução a
função Ψ denominada função de onda. Primeiramente, mostramos a partícula na caixa
unidimensional e, como soluções aceitáveis, a energia potencial é nula na região entre as
paredes da caixa. Para este caso a abstração do modelo está em considerar como hipótese
inicial a equação (7.5) e, a partir daí, a formulação do modelo tem como consequência a
Equação de Schrödinger considerando a energia potencial V nula e, portanto, a equação
2
2 d Ψ
dada é 2m + Ek Ψ = 0. Tendo a equação modelada, passamos à etapa resolução, neste
dx2
caso, analiticamente, de acordo com o Teorema 2.4.2. Apresentamos como aplicação o
processo de tunelamento que é conhecido, também, como efeito túnel. Esse efeito pode ser
encontrado nos semicondutores chamados diodos túneis.
Para a partícula na caixa bidimensional, mostramos a Equação de Schrödinger para

um espaço de duas dimensões. Apesar deste trabalho estar concentrado nas equações dife-
renciais ordinárias, este modelo é dado por equações diferenciais parciais, porém, de sorte,
a resolução deste modelo é dada por separação de variáveis que consiste em transformar o
problema inicial em duas EDO uma para cada variável, isto é, uma dependendo de x e a
outra de y. Este par de equações diferenciais é dado por:
1 d2 χ 2mEx
=− 2
χ dx 2
1 d2 ϕ 2mEy
=− 2
ϕ dy 2
Como última aplicação, apresentamos o oscilador harmônico simples. Como abs-
tração, utilizamos como hipótese a equação (7.8) e, sabendo que a energia potencial é
V (x) = 12 kx2 , substituímos esse valor na equação mencionada e, portanto, temos a formu-
lação do modelo descrito na equação (7.16). Como resolução, transformamos a equação
(7.16) numa equação denominada equação de Hermite cuja solução é dada pelos polinômios
de Hermite.
135
8 Considerações Finais
Neste trabalho, apresentamos modelos matemáticos descritos por Equações Di-

ferenciais Ordinárias para diferentes áreas do conhecimento. A análise dos modelos foi
dada a partir das etapas propostas por Bassanezi (2002) para a obtenção de um modelo
matemático.
No Capítulo 2, foram apresentados alguns conceitos e teoremas sobre as Equações
Diferenciais Ordinárias que foram utilizados para descrever e solucionar os modelos
mostrados nos capítulos subsequentes. Para um aprofundamento do estudo desses tipos de
equações sugerimos a leitura de Boyce e Diprima (2006) e Zill e Cullen (2001) que também
trazem outros modelos matemáticos diferentes dos apresentados neste presente trabalho.
No Capítulo 3, apresentamos o conceito de modelo matemático segundo Bassanezi
(2002) e as etapas que ele descreve para obter o modelo estudado para assim adequá-las a
cada situação apresentada nos modelos estudados.
No Capítulo 4, as vibrações livres amortecidas e não-amortecidas foram apresentadas
como aplicações das EDO na Engenharia Mecânica bem como os circuitos elétricos RC, RL,
LC e RLC na Engenharia Elétrica. Podemos, também, analisar os modelos matemáticos
de ambas as áreas por meio de comparações. Apesar de serem áreas distintas, os modelos
focalizados têm princípios muito comuns que é a oscilação do sistema mecânico e do
sistema elétrico e então foi possível utilizar uma única análise para fornecer respostas para
dois modelos aparentemente distintos. As Tabelas 1 e 2 apresentaram um conjunto de
comparações neste sentido. É importante ressaltar que procuramos fornecer as resoluções
de maneiras distintas, por exemplo, solucionando os problemas de valores iniciais tanto
pela Transformada de Laplace quanto pelos teoremas apresentados no Capítulo 2 que
também contém um breve estudo das transformadas. Um olhar tanto algébrico quanto
gráfico foi utilizado para nos mostrar visões distintas dos fenômenos analisados e após as
aplicações resumimos brevemente as etapas para a obtenção de um modelo proposto por
Bassanezi (2002).
No Capítulo 5, apresentamos um modelo matemático aplicado à Física a respeito
do movimento de queda livre com influência do atrito. Primeiramente, foi trazida uma
abordagem histórica do Cálculo Diferencial e Integral que foi inventado para explicar e
solucionar fenômenos naturais. Sua formalização deu início com Newton e Leibniz, mas
sabemos que seus primórdios estão concentrados, como já citado, nos trabalhos de Eudoxo
e Arquimedes, por exemplo. Esse modelo é abordado sem a resistência de um fluido, por
exemplo, o ar, porém demos um enfoque com essa influência na qual recai numa EDO de
Bernoulli.
136 Capítulo 8. Considerações Finais
No Capítulo 6, foram apresentados dois modelos matemáticos, sendo o primeiro, o

Modelo de Gompertz para crescimento de tumores comparando-o com uma abordagem via
Modelo Logístico e o segundo, Modelo matemático para absorção de drogas. No primeiro,
procuramos abordar o crescimento de tumores de acordo com o estudo de Domingues (2011)
e uma nova abordagem com os mesmos dados pelo Modelo Logístico. A partir disso, tendo
ambos os dados, procuramos compará-los e observar as diferenças existentes. No segundo
modelo, abordamos a queda da concentração da medicação na corrente sanguínea via
modelo simples em que sua solução dá-se por separação de variáveis. Para uma aplicação,
buscamos informações do medicamento Fenobarbital e estudamos algébrica e graficamente
a queda de concentração desta medicação na corrente sanguínea.
No Capítulo 7, um estudo de modelos matemáticos na Química foi apresentado. Es-
tudamos a função de onda de acordo com a Equação de Schrödinger para qualquer sistema.
Primeiramente, mostramos as soluções aceitáveis da partícula na caixa unidimensional e
após isso, apresentamos a partícula na caixa bidimensional; modelada matematicamente
por uma equação diferencial parcial a qual resolvemos por separação de variáveis trans-
formando, assim, o problema em duas equações diferenciais ordinárias, uma para cada
variável em questão. O processo de tunelamento foi uma outra aplicação apresentada.
O efeito túnel consiste em uma partícula que atravessa uma barreira de potencial não
ficando totalmente confinada dentro da caixa como as aplicações da partícula na caixa
unidimensional e bidimensional. O diodo túnel é um semicondutor rápido que trabalha
na casa dos GHz e é um tipo de aplicação do efeito túnel. Finalmente, apresentamos o
oscilador harmônico simples na mecânica quântica que analisa a oscilação em sistemas
microscópicos. Por exemplo, o movimento de vibração de dois átomos em uma molécula
diatômica pode ser representado por esse tipo de oscilador. A modelagem desse sistema
é dada pela Equação de Schrödinger sendo a energia potencial deste sistema fornecida
por V (x) = 12 kx2 . Dado que a função V é um polinômio de grau dois, então tal energia
é denominada de “energia potencial parabólica”. A solução da equação diferencial que
descreve este sistema é dada pelos polinômios de Hermite.
Como já dito anteriormente, procuramos, nesta pesquisa, apresentar alguns modelos
matemáticos em cada área do conhecimento abordada. Sugerimos uma continuação da
pesquisa por meio de outros modelos em outras áreas ou até o aprofundamento nas áreas
apresentadas aqui. Há a possibilidade de fornecer, também, um tratamento numérico para
solucionar as equações diferenciais ordinárias que descrevem os fenômenos apresentados
e não apenas um tratamento analítico como foi mostrado. Trabalhamos com modelos
contínuos, porém há a possibilidade de trabalhar com modelos discretos, por exemplo, o
modelo logístico de diferença dado por yt+1 = yt (r − dyt ), aplicações de sistemas lineares
e não lineares de diferenças para determinar a interação entre duas espécies em que
há um sistema hospedeiro e um parasita. No caso de sistemas de equações diferenciais
no caso contínuo é possível trabalhar concomitantemente com os conceitos de Álgebra
137
Linear, como auto-valores e auto-vetores tanto no caso real como no complexo. Um modelo
matemático descrito por um sistema de equações diferenciais é o predador-presa que
descreve a interação de duas espécies diferentes de animais em um mesmo ambiente ou
até mesmo os chamados modelos de competição em que duas espécies de animais que
ocupam o mesmo ecossistema competem não entre si, mas pela mesma fonte de recursos,
por exemplo, alimentos e habitat.
Por intermédio de nossas questões iniciais, principalmente aquela que verte na
utilização dos conceitos matemáticos profissionalmente, apresentamos modelos descritos
por equações diferenciais que explicam fenômenos que são estudados por profissionais de
cada área. A conceitualização matemática e a sua possível intervenção prática por meio da
modelagem de sistemas fornece ao estudante possíveis subáreas em que há a possibilidade
de estudar projetos e aplicá-los. Um exemplo é o circuito RLC em que é possível observar,
em um osciloscópio, o comportamento da oscilação por causa da dissipação de energia
no resistor, ou seja, por meio de um aparelho usado na prática é possível confirmar e
comparar os cálculos realizados teoricamente.
Do que apresentamos nesta pesquisa, sugerimos, para cada modelo nas áreas do
conhecimento abordadas, a criação de sequências didáticas para o ensino de EDO por
meio da Modelagem Matemática. Como dito na Introdução, não estamos defendendo o
ensino deste tópico matemático sem seu enfoque formal, mas que tenha a possibilidade de
apresentar aplicações das equações diferenciais, pois elas surgiram para resolver problemas
de outras ciências.
139
Referências
ATKINS, P.; PAULA, J. de. Físico - Química. 8. ed. Rio de Janeiro: LTC, 2008.
BARBOSA, J. C. Modelagem na educação matemática: contribuições para o debate

teórico. Emanped, Caxambú, 2001.
BASSANEZI, R. C. Modelagem matemática: uma disciplina emergente nos programas de

formação de professores. Biomatemática, Campinas, n. 9, p. 9–22, set. 1999.
BASSANEZI, R. C. Equações Diferenciais Ordinárias - Um curso introdutório. 1. ed. São

Paulo: UFABC, 2002.
BASSANEZI, R. C. Ensino - aprendizagem com modelagem matemática: uma nova

estratégia. 3. ed. São Paulo: Contexto, 2011.
BATTAGLIN, P. D.; BARRETO, G. Revisitando a história da engenharia elétrica.

Revista de Ensino de Engenharia, Campinas, v. 30, n. 2, p. 49–58, 2011.
BIEMBENGUT, M. S. Modelagem matematica no ensino. 5. ed. São Paulo: Contexto,

2011.
BOYCE, W. E.; DIPRIMA, R. C. Equações diferenciais elementares e problemas de

valores de contorno. Rio de Janeiro: LTC, 2006.
BOYER, C. A. História da Matemática. 2. ed. São Paulo: Edgar Blücher, 2010.
BURAK, D. Modelagem matemática: uma metodologia alternativa para o ensino de

Matemática na 5a série. 188 f. Dissertação (Mestrado em Educação Matemática) —
Universidade Estadual Paulista “Júlio de Mesquisa Neto” - UNESP - Câmpus Rio Claro,
Rio Claro, 1987.
CARVALHO, H. M. de. Modelagem matemática: elementos históricos sobre seu

desenvolvimento em cursos de pós-graduação. 66 f. Dissertação (Mestrado em Matemática
Universitária) — Instituto de Geociências e Ciências Exatas, Câmpus Rio Claro,
Universidade Estadual Paulista “Júlio de Mesquisa Neto” - UNESP, Rio Claro, 2010.
CHAPRA, S. C.; CANALE, R. P. Métodos numéricos para engenharia. São Paulo:

McGraw-Hill, 2008.
DOMINGUES, J. S. Análise do modelo de gompertz no crescimento de tumores sólidos e

inserção de um fator de tratamento. Biomatemática IMECC - Unicamp, Campinas, n. 21,
p. 103–112, 2011.
DONANGELO, R. J.; CAPAZ, R. B. Introdução à Mecânica Quântica. vol. 3. Rio de

Janeiro: Fundação CECIERJ, 2009.
EVES, H. Introdução à História da Matemática. 1. ed. Campinas: Unicamp, 2004.
FICHERA, G.; LACAGNINA, M.; PETRONE, F. Modelling of torsion beam rear

suspension by using multibody method. Multibody System Dynamics, Netherlands, v. 12,
n. 4, p. 303–316, 2004.
140 Referências
FIGUEIREDO, D. G.; NEVES, A. F. Equações diferenciais aplicadas. 3. ed. Rio de

Janeiro: IMPA, 2012.
FRANCO, N. B. Cálculo numérico. São Paulo: Pearson Prentice Hall, 2006.
FRIBERG, S.; MATTSON, S. On the growth rates of human malignant tumors:

implications for medical decision making. Journal of Surgical Oncology, Stockholm, n. 65,
p. 284–297, abr. 1997.
GILLESPIE, T. D. Fundamentals of vehicle dynamics. Warrendale: Society of Automotive

Engineers Inc, 1992.
HALLIDAY, D.; WALKER, J.; RESNICK, R. Fundamentos de física. vol. 2. 8. ed. Rio de
Janeiro: LTC, 2010.
HALLIDAY, D.; WALKER, J.; RESNICK, R. Fundamentos de física. vol. 3. 8. ed. Rio de
Janeiro: LTC, 2010.
INCA. Ações de enfermagem para o controle do câncer: uma proposta de integração

ensino-serviço. 3. ed. Rio de Janeiro: INCA, 2008.
KENNETT, B. A short history of suspension. 2000. Disponível em: <http:

//www.citroenet.org.uk/miscellaneous/suspension/suspension2.html>. Acesso em: 22 mar.
2014.
KRAIGE, L.; MERIAN, J. Mecânica Dinâmica. 4. ed. Rio de Janeiro: LTC, 1997.
LIMA, C. Mecânica Quântica de Schrödinger. Edição de janeiro de 2009. Minas Gerais:

Notas de Aula, 2009.
MODELO. In: DICIONáRIO da Língua Portuguesa. São Paulo: iDicionário Aulete,

2014. Disponível em: <http://aulete.uol.com.br/site.php?mdl=aulete_digital&op=
loadVerbete&pesquisa=1&palavra=modelo>. Acesso em: 25 mar. 2014.
NETO, J. B. Mecânica Newtoniana, Lagrangeana e Hamiltoniana. 1. ed. São Paulo:

Livraria da Fisica, 2004.
OGATA, K. Engenharia de controle moderno. 4. ed. Rio de Janeiro: Pearson - Prentice

Hall, 2003.
ZILL, D. G.; CULLEN, M. R. Equações diferenciais. 3. ed. São Paulo: Pearson Makron
Books, 2001.

Aplicações Nas Ciencias

Enviado por

Direitos autorais:

Formatos disponíveis

Aplicações Nas Ciencias

Enviado por

Dados do documento

Título original

Direitos autorais

Formatos disponíveis

Compartilhar este documento

Compartilhar ou incorporar documento

Opções de compartilhamento

Você considera este documento útil?

Este conteúdo é inapropriado?

Direitos autorais:

Formatos disponíveis

Aplicações Nas Ciencias

Enviado por

Direitos autorais:

Formatos disponíveis

29

3.1 Modelo Matemático

Modelo matemático é um conjunto de símbolos e relações matemáticas

Destacamos a deﬁnição do Professor Bassanezi, pois ele é um dos pioneiros no

Se formularmos a hipótese de que a variação populacional é proporcional à população

a Matemática tem penetrado fortemente na Economia, Química, Biologia,

A partir das necessidades em um determinado momento, os proﬁssionais dessas áreas

velocidades consideradas baixas, mas Einstein modelou o sistema considerando velocidades

A modelagem matemática constitui-se em um conjunto de procedimentos

4 Modelagem Matemática nas En-

4.1 Aplicações das EDO na Engenharia Mecânica

porque focam-se muito no aprendizado por memorização, o que é muito diferente do

4.1.1 Vibração livre não-amortecida

A expressão (4.2) é uma Equação Diferencial de Segunda Ordem e é conhecida

vimos no Capítulo 2, a solução dessa EDO Linear Homogênea de 2a Ordem é caracterizada

Podemos, na equação acima, fazer uma substituição conveniente. Ao analisarmos as

x = Acos (ωn t) + Bsen (ωn t) , A, B ∈ R (4.5)

Figura 2 – Mola em equilíbrio.

Figura 3 – Mola após sofrer uma tensão (esticamento).

Para a análise do modelo de acordo com as etapas de modelagem descritas no

4.1.2 Vibração livre amortecida

mx + cx + kx = 0. (4.6)

Novamente, temos um fenômeno modelado matematicamente por uma Equação Diferencial

x + 2ζωn x + ωn2 x = 0. (4.7)

A equação característica é dada por

1. Superamortecido. Quando ζ > 1 a equação característica fornece duas raízes reais λ1

2. Criticamente amortecido. Quando ζ = 1 a equação característica fornece duas raízes

x(t) = e−ωn t (A + Bt) . (4.9)

3. Subamortecido. ζ < 1 a equação característica fornece duas raízes complexas, uma

x(t) = e−ζωn t [Ccos(ωf t) + Dsen(ωf t)] (4.11)

em que C = (A + B) e D = i (A − B), em que A, B ∈ R

E, por ﬁm, a avaliação é dada quando Δ é estritamente negativo. Neste caso,

x (t) + ωn2 x(t) = 0

L {x } + ωn2 L {x} = 0

Por ﬁm, fazendo os ajustes necessários e aplicando a transformada inversa, obtemos

mx (t) + kx(t) = 0 =⇒ 3x (t) + 588x(t) = 0

=⇒ x (t) + 196x(t) = 0

Portanto, a função que descreve o movimento livre do sistema é x(t) = 0, 10cos(14t)+

Figura 4 – Gráﬁco: deslocamento versus tempo.

em que Em é a energia mecânica, Ec é a energia cinética e Ep é a energia potencial.

mx (t) + ke−αt x(t) = 0.

Exemplo 4.1.2. Para o movimento subamortecido, quando ζ < 1, utilizaremos a Trans-

x (t) + 2ζωn x (t) + ωn2 x = 0

L {x } + 2ζωn L {x } + ωn2 L {x} = 0

s2 L {x} − sx0 − V0 + 2ζωn sL {x} − 2ζωn x0 + ωn2 L {x} = 0

Reordenando os termos, temos

No denominador podemos completar quadrado da seguinte forma:

s2 + 2ζωn s + ωn2 + (ωn2 ζ 2 − ωn2 ζ 2 ) = s2 + 2ζωn s + ωn2 ζ 2 + ωn2 − ωn2 ζ 2

Fazendo os ajustes necessários para poder aplicar a Transformada inversa de

Exemplo 4.1.3. Para o movimento com amortecimento crítico, basta fazer ζ = 1 em

sx0 + 2ζωn x0 + V0 sx0 + 2ωn x0 + V0

Aplicando a Transformada inversa de Laplace, temos

x(t) = e−ωn t [x0 + (ωn x0 + V0 )t] (4.14)

sx0 + 2ζωn x0 + V0 sx0 + 2ζωn x0 + V0

Utilizando as frações parciais, podemos reescrever L {x} da seguinte forma

Encontrando A e B a partir das condições iniciais já estabelecidas e aplicando a

Nesta seção, apresentamos tópicos de onde as equações diferenciais ordinárias podem

4.1.4 Suspensões mecânicas automotivas

Por volta do século XVIII, os processos de fabricação eram mais desenvolvidos,