Eq Diferenciais Ordinarias

EQUAÇÕES DIFERENCIAIS ORDINÁRIAS E APLICAÇÕES

Aluno: Marcelo Luiz Freitas Fogal

Orientador: Ernandes Rocha de Oliveira

Ilha Solteira, 1999

UNIVERSIDADE ESTADUAL PAULISTA “JÚLIO DE MESQUITA FILHO”

Campus de Ilha Solteira Departamento de Matemática

UNESP - FEIS

Estudo de Equações Diferenciais Ordinárias

1



UNESP - FEIS

Introdução

Este trabalho consiste de um estudo de equações diferenciais

ordinárias e de seus métodos de determinação de suas soluções. É bem conhecido

que muitos fenômenos que interessam às Engenharias e outras ciências podem ser

estudadas através de modelos matemáticos nos quais aparecem de modo

importante equações ordinárias. Processos contínuos que envolvam a análise de

taxa de variação, são geralmente descritos por meio da noção da derivada. Entende-

se também que o estudo de modelos matemáticos simples, porém significativos,

permite ao iniciante na matéria compreender melhor o poder e o limite dos métodos

matemáticos utilizados, além disso tais modelos podem servir como um primeiro

passo na busca de formação matemática necessária para que se possa desenvolver

uma confiança na formulação e exploração de novos modelos.


2



UNESP - FEIS

Objetivos

O objetivo principal deste trabalho é o de estudar, tanto do ponto de

vista qualitativo, quanto do ponto de vista quantitativo, as principais classes de

equações diferenciais ordinárias de primeira ordem e as de segunda ordem lineares,

além de alguns problemas significativos modelados por tais equações.

A execução deste trabalho deverá proporcionar uma experiência inicial

com a teoria de equações diferenciais ordinárias, modelagem matemática e o uso de

um software para realizar simulações e visualização de soluções.


3



UNESP - FEIS

Observações Históricas

O desenvolvimento das equações diferenciais está diretamente ligado

com o desenvolvimento da matemática. A fim de se ter uma certa perspectiva

histórica, vamos traçar algumas tendências principais na história do problema e

identificar as personalidades contribuidoras mais eminentes.

O estudo das EDO's inaugurou-se no início do cálculo, com Isaac

Newton (1642-1727) e Leibniz (1646-1716), no século XVII. Embora Newton tenha

trabalhado relativamente pouco no campo das equações diferenciais, o

desenvolvimento que proporcionou ao cálculo e à elucidação dos princípios básicos

da mecânica construíram a base para as aplicações que se fizeram no século XVIII,

por Euler. Newton classificou as equações diferenciais de primeira ordem de acordo

como as formas dy/dx = f(x), dy/dx = f(y) e dy/dx = f(x,y).

Leibniz chegou aos resultados fundamentais do cálculo por via

independente, embora um pouco posterior a Newton, mas foi o primeiro a publicá-

los, em 1684. Leibniz tinha plena consciência do poder de uma boa notação

matemática e a notação que usamos para a derivada (dy/dx), e para a integração,

foram introduzidas por ele. Descobriu o método de separação das variáveis em

1691, a redução de equações homogêneas e equações separáveis, e o

procedimento de resolução de equações lineares de primeira ordem em 1694.

Os irmãos Jakob (1654-1705) e Johann (1667-1748) Bernoulli,

contribuíram muito para o desenvolvimento de métodos de resolução de equações

diferenciais e para ampliar o campo de aplicação destas equações. Com o auxílio do

cálculo formularam como equações diferenciais muitos problemas de mecânica e os

resolveram. Por exemplo, Jakob Bernoulli resolveu a equação diferencial y’ = [a³/(b²y

- a³)]½ em 1690 e, no mesmo artigo usou pela primeira vez o termo "integral" no

sentido moderno. Em 1694, Johann Bernoulli resolveu a equação dy/dx = y/ax,

embora não se soubesse na época, que d(ln x) = dx/x.

Euler teve especial interesse na formulação de problemas de

mecânica em linguagem matemática e o desenvolvimento de métodos de resolução

destes problemas matemáticos, identificou a condição de exatidão das equações

diferenciais de primeira ordem em 1734-1735, desenvolveu a teoria dos fatores de


4



UNESP - FEIS

integração e apresentou a solução geral das equações lineares com os coeficientes

constantes em 1743.

No que se refere às equações diferenciais elementares, Joseph-Louis

Lagrange (1736-1813) mostrou em 1762-1765 que a solução de uma equação

diferencial homogênea de ordem n é uma combinação linear de n soluções

independentes. Depois, em 1774-1775, publicou o desenvolvimento completo do

método da variação de parâmetros. Lagrange também é conhecido pelo seu

tratamento fundamental nas equações diferenciais parciais e no cálculo das

variações.

A equação de Laplace é fundamental em muitos ramos da física

matemática, e Pierre-Simon de Laplace estudou-a profundamente em suas

investigações da atração gravitacional. A transformada de Laplace também recebeu

o nome em sua homenagem, embora sua utilidade para a solução de equações

diferenciais só tenha sido reconhecida muito mais tarde.

As diversas equações diferenciais que resistiram à resolução por meios

analíticos levaram à investigação de métodos numéricos de aproximação. Na altura

de 1900, métodos de integração numérica, muito eficientes, já tinham sido

elaborados, mas a implementação destes métodos estava severamente restringida

pela necessidade de execução de cálculos manuais ou com equipamento de

computação muito primitivo. Nos últimos cinqüenta anos, o desenvolvimento de

computadores cada vez mais poderosos e versáteis ampliou a gama de problemas

que podem ser investigados com eficiência por meio de métodos numéricos. Durante

o mesmo período, desenvolveram-se integradores numéricos muito refinados e

robustos, que se encontram em todos os centros de computação científica.

Uma outra característica das equações diferenciais no século XX foi a

criação de métodos geométricos ou topológicos, especialmente para as equações

não-lineares. Assim, embora as equações diferenciais sejam um tema antigo, a

respeito do qual seja grande o conhecimento, tornou-se no final do século XX uma

fonte de problemas fascinantes e importantes, ainda não resolvidos.


5



UNESP - FEIS

Classificação das EDO's

Muitos problemas significativos da engenharia, das ciências físicas e

das ciências sociais, formulados em termos matemáticos, exigem a determinação de

uma função que obedece a uma equação que contém uma ou mais derivadas da

função desconhecida. Estas equações são equações diferenciais.

Temos as equações diferenciais ordinárias e equações diferenciais

parciais. Uma das classificações se baseia em a função desconhecida depender de

uma só variável ou de diversas variáveis. No primeiro caso, na equação diferencial

só aparecem derivadas ordinárias e a equação é uma equação diferencial ordinária.

No segundo caso, as derivadas são derivadas parciais, e a equação é uma equação

diferencial parcial.

Como um exemplo para equações diferenciais ordinárias, temos:

)()(

tKRdt

tdR−= , onde K é uma constante conhecida. Um exemplos

típico de equação diferencial parcial é a equação do potencial:

0),(),(

2

2

2

2

=+dy

yxud

dx

yxud

A equação do potencial aparece em muitos problemas de eletricidade

e magnetismo.

Sistemas de equações diferenciais é uma outra classificação que

depende do número de funções desconhecidas que estão envolvidas. Quando forem

duas ou mais as funções desconhecidas, é necessário ter um sistema de equações.

A ordem de uma equação diferencial é a ordem de derivada de maior

ordem que aparece na equação. De uma forma mais geral, a equação F[t, u(t),

u’(t),…, u(n)(t)] = 0 é uma equação diferencial de ordem n. Esta equação pode ser

escrita como F(t, y, y’,…, y(n)) = 0.

Por exemplo, y’’’ + 2ety’’ + yy’ = t4 é uma equação diferencial de terceira

ordem em y = u(t). Admitindo-se que é sempre possível resolver uma dada equação

diferencial ordinária na derivada de ordem mais elevada e ter y(n) = f(t, y, y’, y’’,…,

y(n-1)).

Uma solução desta equação no intervalo α<t<β, é uma função φtal que

φ’, φ’’,…, φ(n) existem e satisfazem a :


6



UNESP - FEIS

φ(n)(t) = f[t, φ(t), φ’(t),…,φ(n-1)(t)], para todo t em α<t<β.

Uma classificação importante das equações diferenciais é a que as

divide em lineares e não-lineares. A equação diferencial ordinária F(t, y, y’,…, y(n)) =

0 é linear se F for uma função linear das variáveis y, y’,…,y(n). Definição semelhante

aplica-se às equações diferenciais parciais. Assim, a equação diferencial ordinária

linear, de ordem n, é:

a0(t)y(n) + a1(t)y

(n-1) +…+ an(t)y = g(t).

Uma equação que não tenha essa forma é uma equação não-linear.

Um problema simples que leva a uma equação não-linear é o do pêndulo. O ângulo

θ que um pêndulo de comprimento L faz com a direção vertical obedece à equação

não-linear:

0sen2

2

=+ θθ

L

g

dt

d.

A teoria matemática e as técnicas correspondentes para a resolução

das equações lineares estão muito desenvolvidas, mas em contraposição, para as

equações não-lineares a situação não é tão satisfatória. Em boa parte faltam

técnicas gerais para a resolução de equações não-lineares, e a teoria associada a

estas equações também é mais complicada do que a teoria das equações lineares.

Por isso, é bom que muitos problemas importantes levem a equações diferenciais

ordinárias lineares ou, pelo menos em primeira aproximação, a equações lineares.

Campos de Direções

Uma interpretação geométrica das equações diferenciais e das

respectivas soluções são os campos de direções. O ponto de vista geométrico é

especialmente útil no caso de equações de primeira ordem, isto é, de equações com

a forma:

),( ytfdt

dy= (a1)

A solução da equação (a1) é uma função y = φ(t), a solução geométrica

de uma solução é o gráfico de uma função. Geometricamente a equação (a1) afirma

que, em qualquer ponto (t,y), o coeficiente angular dy/dt da solução neste ponto é

dado por f(t,y). Podemos representar graficamente esta situação traçando um

pequeno segmento de reta, no ponto (t,y), com o coeficiente angular f(t,y). O


7



UNESP - FEIS

conjunto de segmentos de reta é o campo de direções da equação diferencial (a1).

O campo de direções pode ser visualizado pelo desenho de pequenos segmentos

de reta num conjunto representativo de pontos do plano ty. Embora este traçado seja

tedioso de ser feito manualmente, é tarefa simples para um computador, pois exige

somente o cálculo repetido de f(t,y) para diferentes valores de t e de y. Com o

desenho de campo de direções, pode-se perceber muitas vezes, o comportamento

qualitativo das soluções ou observar regiões no plano que tenha interesse especial.


8



UNESP - FEIS

Equações Diferenciais de Primeira Ordem

Equações Lineares

Se a função da equação dy/dt = f(t,y) depende linearmente da variável

dependente y, então a equação pode ser escrita na forma:

)()( tgytpdt

dy=+ , e é chamada de equação diferencial linear de

primeira ordem. Vamos admitir que p e g são funções conhecidas e contínuas no

intervalo α < x < β. Por exemplo, a equação diferencial:

2

3

2

1=+ y

dt

dy, é uma equação linear particularmente simples, com as

funções p(x) = 21 e g(x) = 2

3 , ambas constantes.

Temos como exemplo a resolução da equação 2

3

2

1=+ y

dt

dy e

determinar como as soluções se comportam para grandes valores de x. Para

resolver esta equação, observamos que se y ≠ 3 podemos rescrever a equação na

forma 2

3−−=

y

dt

dy, daí então

2

1

3−=

−y

dtdy

.

Uma vez que o primeiro membro da equação2

1

3−=

−y

dtdy

, é a derivada

de 3ln −y , temos 2

13ln −=−y

dt

d. Segue-se então que se:

ct

y +−=−2

3ln , onde c é uma constante de integração arbitrária.

Portanto, com a forma exponencial de ambos os membros, obtemos:

)2exp()exp(3 tcy −±=− , onde )exp(cc ±= também é uma constante

arbitrária não nula.

Fator integrante

Inicialmente, escrevemos )exp(rtr

Ky +−= (I) na forma


9



UNESP - FEIS

crtr

Kye

rt +−−=− )exp( (II), e depois derivando os dois membros em

relação a t, vem:

)exp()exp()'( rtKrtryy −=−− (III), que é equivalente a equação

Krydt

dy+= (IV).

Observe que agora podemos resolver a equação (IV) invertendo os

passos precedentes. Transpondo o termo ry para o lado esquerdo da equação e

multiplicando por e-rt, obtemos a equação (III). Note que o lado esquerdo de (III) é a

derivada de ye-rt, de modo que a equação se torna:

)exp())'exp(( rtKrty −=− (V)

Finalmente integrando os dois membros da equação (V), obtemos a

equação (II) e portanto a solução (I).

Lidando agora com a questão da equação )()( tgytpdt

dy=+ (VI). Nosso

objetivo é multiplicar a equação diferencial (VI) por um fator integrante apropriado e

assim colocá-la em uma forma integrável. Para determinar esse fator integrante

primeiro multiplicamos a Eq.(VI) por uma função µ(t), ainda indeterminada. Temos

agora:

µ(t)y’ + µ(t)p(t)y = µ(t)g(t) (VII)

Devemos agora reconhecer o lado esquerdo da Eq.(VII) como a

derivada de alguma função. O fato é que existem dois termos e um dos termos é

µ(t)y’ sugere que o lado esquerdo da equação (VII) pode ser a derivada do produto

µ(t)y. Para que isto seja verdade, o segundo termo do lado esquerdo da Eq.(VII),

µ(t)p(t)y, deve ser igual a µ’(t)y. Isto significa que µ(t) deve satisfazer a equação

diferencial:

µ’(t) = p(t)µ(t) (VIII)

Se admitirmos que µ(t) é positiva, podemos escrever a equação (VIII)

como:

)()(

)('tp

t

t=

µ

µ, ou )()(ln tpt

dt

d=µ , então integrando ambos os termos,

∫ += Kdttpt )()(ln µ .


10



UNESP - FEIS

Pela escolha da constante K arbitrária como zero, temos a função u

mais simples possível, ou seja ∫= dttpt )(exp)(µ (IX).

De forma que µ(t) é positiva para todos os t conforme admitimos.

Depois de determinar-mos o fator integrante u(t), voltamos a equação

)()( tgytpdt

dy=+ e multiplicamos por µ(t), obtendo assim a equação (VII). Como µ

satisfaz à Eq.(VIII), a Eq.(VII) se reduz a

[µ(t)y]’ = µ(t)g(t) (X)

Integrando ambos os membros da equação (X) obtemos:

∫ += cdttgtyt )()()( µµ , ou )(

)()(

t

cdttgty

µ

µ∫ += (XI)

Uma vez que y representa qualquer solução da Eq.(VI), concluímos

que toda solução da Eq.(VI) está incluída no segundo membro da Eq.(XI). Portanto,

esta expressão é uma solução geral da Eq.(VI). Observe que para encontrar a

solução dada pela Eq.(XI) são necessárias duas integrações, uma para ter µ(t) pela

Eq.(IX) e outra para determinar y pela Eq.(XI).

Para se determinar o fator integrante µ(t), é necessário ter certeza que

a equação diferencial tem exatamente a forma (VI).

Interpretando geometricamente a Eq.(XI) nota-se que ela é de uma

família infinita de curvas, uma para cada valor de c. Estas curvas são as curvas

integrais da equação diferencial. Muitas vezes é importante selecionar um membro

particular da família de curvas integrais, o que se faz pela identificação de u ponto

particular (t0, y0) por onde deve passar uma das curvas da solução. Isto se escreve

como y(t0) = y0, e é conhecida como uma condição inicial.

Exemplo

Determinar o valor de y0 para o qual a solução do problema de valor

inicial permaneça finita quando t → ∞. Dado:

y’ –y = 1 + 3sent, y(0) = y0

Inicialmente determinamos o fator integrante:

∫ −= )exp()( dttµ , )exp()( tt −=µ


11



UNESP - FEIS

Agora multiplicando todos os termos da equação pelo fator integrante

temos:

tttyttyt sen3)exp()exp()exp()(')exp( −+−=−−−

Identificando o lado esquerdo da equação como sendo a derivada do

produto e integrando ambos os lados temos:

∫∫ +−+−−=−

ctttdt

tytdsen)exp(3)exp(

))()(exp(, sendo c uma constante.

Dividindo ambos os lados por exp(-t), temos:

)exp()sen(cos2

31)( tcttty ++−−=

Como nossa condição inicial é y(0) = y0, substituindo t por zero na

equação acima temos:

cy +−−=2

31)0( , portanto cy +−=

2

50 e logo

2

50 += yc

Com isso, )exp()2

5()sen(cos

2

31)( 0 tyttty +++−−= , e como queremos

que a solução do problema permaneça finita, devemos igualar o termo multiplicativo

da exponencial a zero, tendo assim:

02

50 =+y , o que nos leva a

2

50 −=y .


12



UNESP - FEIS

Aplicações das equações lineares de primeira ordem

Exemplo 1

A pessoa A abre uma conta remunerada com 25 anos, deposita

R$2000,00 por ano durante 10 anos e depois disso não faz mais nenhum depósito. A

pessoa B espera até os 35 anos para abrir a sua conta remunerada, mas deposita

R$2000,00 por ano durante 30 anos. Nos dois casos não há nenhum investimento

inicial.

(a) Supondo uma taxa de juros de 8% ao ano, qual será o saldo das

duas contas quando os titulares tiverem 65 anos?

(b) Para uma taxa de juros constante mas não especificada r,

determine o saldo das duas contas quando os titulares tiverem 65 anos em função

de r.

Resolução: (a)

Pessoa A: 100 ≤≤ t

=

+=

0)0(

2000)(100

8)('

a

aa

S

tStS , queremos saber Sa(10) = ?

Resolvendo a equação através do fator integrante, obtemos:

)1)(exp(2000

)( −= rtr

tS a , agora substituindo em t o tempo de 10 anos,

temos: )1)10(exp(2000

)10( −= rr

S a 52,30638$R= , portanto, passados 10 anos a

pessoa A possuirá este valor.

Agora temos: 3010 ≤< t

rtSatS a )()(' = , achando a solução desta equação obtemos:

)exp()( rtKtS a = , sendo K uma constante, que em nosso caso, é

exatamente Sa (10) = R$ 30638,52. Substituindo temos:

81,337733$)30exp(52,30638)30( RtS a ==

Pessoa B:


13



UNESP - FEIS

=

+=

?)30(

2000)(100

8)('

b

b

S

tSbtS

Como já havíamos encontrado a solução desta equação para a pessoa

A, temos:

4,250579$)1)100

8(exp(

1008

2000)1)(exp(

2000)( Rtrt

rtSb =−=−= .

(b)

Expressando em função de r, temos:

Pessoa A: )1)10)(exp(30exp(2000

−rrr

Pessoa B: )1)30(exp(2000

−rr

Uma forma de analisarmos melhor esse problema seria utilizando o

Matlab, e plotarmos as curvas para ambas pessoas:

r=0.001:0.001:0.08;

Sa=(2000./r).*(exp(10*r)-1).*exp(30*r);

Sb=(2000./r).*(exp(30*r)-1);

plot(r,Sb,r,Sa,'r') % 'r' dá cor vermelha a segunda curva

grid on

Exemplo 2

Um estudante universitário faz um empréstimo de R$8000,00 para

comprar um carro. A taxa de juros cobrada pelo banco é de 10% ao ano. Supondo

que os juros sejam capitalizados continuamente e que o estudante amortize a dívida

continuamente a uma taxa anual constante K para que o empréstimo seja pago em

três anos. Determine também o total de juros pagos durante esses três anos.

Resolução:

=

−−=

0)3(

)1)(exp()exp()( 0

S

rtr

KrtStS , onde K é a constante que queremos

determinar.


14



UNESP - FEIS

0)1)3(exp()3exp(8000)3( =−−= rr

KrS , resolvendo esta equação

obtemos:

39,3085$5,387,10798 RKK =⇒=

Determinando agora os juros, temos:

J(t) = Kt – S0, onde S0 é o empréstimo.

17,125680003.84,3085)3( =−=J

Portanto os juros foram de R$ 1256,17.

Exemplo 3

O comprador de uma casa não pode pagar mais do que R$800,00 por

mês de prestação. Suponha que a taxa de juros seja 9% ao ano e que o prazo de

pagamento seja de 20 anos. Suponha também que os juros sejam capitalizados

continuamente e que os pagamentos também sejam feitos continuamente.

(a) Determine o preço máximo que este comprador pode pagar pela

casa.

(b) Determine o total de juros pagos pelo comprador se ele comprar

a casa nas condições do item (a)

Resolução: (a)

Partindo da mesma equação do exemplo anterior:

=

−−=

0)20(

)1)(exp()exp()( 0

S

rtr

KrtStS , onde K é a quantia que o comprador

pagará por ano.

Então: 800.12 = 9600 = K

0)1)20(exp(9600

)20exp()( 0 =−−= rr

rStS , resolvendo:

78,89034)8,1exp(

06,5386290 ==S


15



UNESP - FEIS

Portanto o máximo que o comprador poderá pagar pela casa será R$

89034,78.

(b)

J(t) = Kt – S0 = 9600t – S0

J(20) = 192000 – 89034,78

J(20) = 102965,22

Portanto, o total de juros pagos pelo comprador é de R$ 102965,22.

Exemplo 4

Quais seriam as respostas do problema anterior se o prazo do

empréstimo fosse aumentado para 30 anos?

Resolução: (a)

=

−−=

0)30(

)1)(exp()exp()( 0

S

rtr

KrtStS onde K é a quantia que o comprador

pagará por ano.

Então: 800.12 = 9600 = K

0)1)30(exp(9600

)30exp()( 0 =−−= rr

rStS , resolvendo:

08,99498)7,2exp(

72,14805040 ==S

Portanto o máximo que o comprador poderá pagar pela casa será R$

99.498,08.

(b)

J(t) = Kt – S0 = 9600t – S0

J(30) = 288000 – 99498,08

J(30) = 188501,92

Portanto, o total de juros pagos pelo comprador é de R$ 188501,92.

Exemplo 5


16



UNESP - FEIS

Uma bola com massa de 0,25 Kg, é lançada para cima, com a

velocidade inicial de 20 m/s, do terraço de um edifício com 30 m de altura. Desprezar

a resistência do ar. Calcular a altura máxima que a bola atinge acima do solo.

Resolução:

Partindo de: gamgma −=⇔−=

Como a aceleração é a derivada da velocidade, temos:

gdt

dv−= , agora integrando ambos os lados

dtgdt

dvtt

∫∫ −=00

tvtvgtvtv 8,9)0()()0()( −=⇒−=−⇒ , temos que a

velocidade no instante zero é de 20 m/s.

Temos que quando a bola sobe sua velocidade final é zero, então

substituindo:

sttt 041,28,9

208,9200 =⇒=⇒−=

Sabemos que a velocidade é a derivada do espaço, então:

tdt

dx8,920 −= , integrando ambos os lados

dttdtdt

dxttt

∫∫∫ −=000

8,920 29,420)( tttx −=⇒

Para calcularmos a altura máxima percorrida pelo bola devemos

substituir o tempo gasto para a bola subir. 2)041,2.(9,4)041,2.(20)041,2( −=x

mx 41,20)041,2( = , mas como queremos a altura que a bola atinge

acima do solo e o edifício possuí 30 m, devemos:

Altura máxima = 20,41 + 30 = 50,41 m

Portanto a altura máxima atingida pela bola é de 50,41 m.

Exemplo 6


17



UNESP - FEIS

Vamos admitir as mesmas condições do problema anterior, porém com

uma força resistiva devido ao ar, expressa por 30v , onde v está em m/s. Achar a

altura máxima, acima do solo, atingida pela bola?

Resolução:

=

−−=

20)0(30

v

vmg

dt

dvm

=

−−=⇒

20)0(5,7

8,9

v

v

dt

dv

Achando agora o fator integrante

)152exp()75

10exp()( ttt ==µ , multiplicando todos os membros da

equação pelo fator integrante

)152exp(8,9)15

2exp(15

2)15

2exp( tvtdt

dvt −=+ , identificando o lado

esquerdo da equação como a derivada de um produto e integrando ambos os lados,

temos:

∫ ∫−=t t

dttvtdt

d

0 0

)152exp(8,9))15

2(exp( , resolvendo

)1)152(exp(

2

15)8,9()0()()15

2exp( −−=− tvtvt

5,93)152exp(5,73)()15

2exp( +−= ttvt

)152exp(5,935,73)( ttv −+−=

Para calcularmos o tempo que a bola levou para subir, igualamos a

velocidade final a zero

sttt 805,1

15

2)78,0ln()15

2exp(5,935,730 =⇒−=⇒−+−=

Como sabemos que a derivada do espaço é a velocidade, encontramos

então a equação do espaço para calcularmos a altura máxima que a bola alcançou

)1)152(exp(

2

15)5,93(5,73)( −−−−= tttx


18



UNESP - FEIS

25,701)152exp(25,7015,73)( +−−−= tttx , como queremos a altura

máxima e temos o tempo de subida da bola, substituímos o tempo de subida na

equação do espaço:

25,701)15

)805,1(2exp(25,701)805,1(5,73)805,1( +−−−=x

Portanto x(1,805) = 17,3 m, mas como a bola foi lançada de um edifício

de 30 m, teremos que a altura máxima alcançada pela bola foi de 47,3 m.

Exemplo 7

Um pára-quedista que pesa 180 lb (incluindo-se o equipamento) cai

verticalmente para baixo, de uma altura de 5.000 ft, e abre o pára-quedas depois de

10 segundos de queda livre. Vamos admitir que a força de resistência do ar seja

v75,0 quando o pára-quedas estiver fechado, e v12 quando o pára-quedas estiver

aberto, onde a velocidade v está em pés por segundo.

(a) Achar a velocidade do pára-quedista quando o pára-quedas se

abre.

(b) Achar a distância percorrida em queda livre, até o pára-quedas

se abrir.

(c) Qual a velocidade limite v1 depois de o pára-quedas se abrir?

(d) Estimar o tempo de permanência do pára-quedista no ar depois

de o pára-quedas se ter aberto.

Resolução: (a)

Temos que a aceleração da gravidade é de g = 32ft/s² e que a massa

do pára-quedista é de m = 5,625 slugs. Sendo K1 a constante de resistência antes

do pára-quedas ser aberto.

=

−=

0)0(

1

v

vKmgdt

dvm

=

−=⇒0)0(

1

vm

vKg

dt

dv

, 100 <≤ t

gm

vK

dt

dv=+ 1 , achando agora o fator integrante


19



UNESP - FEIS

)exp()( 1m

tKt =µ , multiplicando todos os termos pelo fator integrante

gm

tK

m

vK

mtK

dt

dv

mtK )exp()exp()exp( 1111 =+ , identificando o lado

esquerdo da equação como a derivada de um produto e integrando ambos os lados,

obtemos:

∫ ∫=t t

gdtm

tKv

mtK

dt

d

0 0

11 )exp())(exp( , resolvendo

)1)(exp()0()()exp( 1

1

1 −=−m

tK

K

mgvtv

mtK , como v(0) = 0

)exp()( 1

11m

tK

K

mg

K

mgtv −−= , agora substituindo os valores para o

tempo de 10 s, temos:

7,176)10( =v , então no instante t = 10, o pára-quedista estará a uma

velocidade de 176,7 ft/s.

(b)

Queremos a distância percorrida pelo pára-quedista nestes 10

segundos. Como conhecemos a equação da velocidade e sabemos que a derivada

do espaço é a velocidade, então:

=

−−=

0)0(

)exp( 1

11

x

mtK

K

mg

K

mg

dt

dx

, queremos x(10) = ?

Integrando ambos os lados, temos:

∫∫ ∫ −−=tt t

mtK

K

mg

K

mg

dt

dx

0

1

10 0 1

)exp(

−

−

−−= 1)exp()( 1

111m

tK

K

m

K

mg

K

mtgtx , agora substituindo os

valores, obtemos:

47,074.1)10( =x , então no instante t = 10, o pára-quedista terá

percorrido uma distância de 1.074,47 ft.

(c)


20



UNESP - FEIS

−=∞→ )exp(

1)(

222lim

mtKK

mg

K

mgtv

t

)exp(7,176 2m

tK−+ , como t está

tendendo ao infinito, então as duas últimas parcelas tenderão a zero, restando

apenas:

sftK

mgtv

t

/1512

625,5)32()(

2lim ===

∞→

Portanto, a velocidade limite será de 15 ft/s.

(d)

Temos que )exp(7,176)exp()( 22

22m

tKm

tK

K

mg

K

mgtv

−+

−−= , e como

sabemos que o espaço é a derivada da velocidade, temos:

)exp(7,176)exp( 22

22m

tKm

tK

K

mg

K

mg

dt

dx −+

−−= , agora integrando

ambos os lados

dtm

tKm

tK

K

mg

K

mg

dt

dxtt

∫∫

−+

−−=

0

22

220

)exp(7,176)exp( , e resolvendo, temos:

−

−

−+

−

−

−−= 1)exp(7,1761)exp()( 2

2

2

222m

tK

K

m

mtK

K

m

K

mg

K

mtgtx

Para encontrarmos o tempo que o pára-quedista ainda permanecerá no

ar, devemos saber qual a distância que ele ainda deve percorrer e para isso

devemos subtrair o que ele já percorreu da distância total:

5000 – 1074,47 = 3925,53, agora substituindo na equação do espaço

essa distância encontrada, acharemos o tempo que ele ainda permanecerá no ar.

−

−

−+

−

−

−−= 1)exp(7,1761)exp(52,3925 2

2

2

222m

tK

K

m

mtK

K

m

K

mg

K

mtg

Substituindo os valores, obtemos t = 256,6 segundos.

Portanto o pára-quedista permanecerá no ar por 256,6 segundos.

Exemplo 8


21



UNESP - FEIS

O nuclídeo radioativo tório 234 desintegra-se a uma taxa proporcional à

quantidade presente. Se 100 mg deste material reduzirem-se a 82,04 mg em uma

semana, achar uma expressão que dê a quantidade presente em qualquer instante.

Achar também o intervalo de tempo necessário para que a massa do material decaia

à metade do seu valor original.

Resolução:

Seja Q(t) a quantidade de tório 234 presente em qualquer instante, t,

onde Q está em miligramas e t em dias. A observação que o tório 234 se desintegra

a uma taxa proporcional à quantidade presente significa que a taxa temporal de

variação dQ/dt é proporcional a Q. Então Q satisfaz a equação diferencial.

rQdtdQ −= , onde a constante r > 0 é a constante de desintegração. A

condição inicial é que: Q(0) = 100.

Temos também uma segunda condição Q(7) = 82,04.

A solução geral da equação diferencial é rtctQ

−= exp)( ,onde c é uma

constante arbitrária. A condição inicial exige que c = 100 e, portanto rttQ

−= exp100)( .

Satisfazendo a segunda condição, temos:

r7exp10004,82 −= , e daí 102828,07

8204,0ln −=−= diasr .

Essa é a constante de desintegração, agora substituindo r na equação,

temos:

mgtQt02828,0exp100)( −= , que dá o valor de Q(t) em qualquer instante.

Chamamos de meia vida o intervalo de tempo que o material se reduz

a metade. Seja α o tempo para que Q(t) seja igual a 50 mg. Pela equação

rttQ

−= exp100)( temos: αr−= exp10050 . Com o valor de r já encontrado temos:

.5,2402828,0

2lndias≅=α

Durante o processo de decaimento, a massa de um material radioativo,

como o tório234, reduz-se e acaba se aproximando de zero.

Exemplo 9


22



UNESP - FEIS

No instante t = 0 um tanque contém Q0 Kg de um certo sal dissolvido

em 100 litros de água. Uma solução de sal em água, com 0,25 Kg de sal por litro de

água, entra no tanque à razão de r litros/minuto e uma solução homogênea sai do

tanque a uma mesma vazão. Formule o problema de valor inicial que descreve esse

processo. Determine a quantidade de sal Q(t) presente no tanque em um dado

instante t e também a quantidade limite Ql que está presente depois de um longo

tempo. Se r = 3 e Q0 = 2Ql, determine o intervalo de tempo T após o qual a diferença

entre a quantidade de sal e Ql é menor que 2%. Determine também a vazão em

litros/minuto para que o valor de T não seja maior do que 45 minutos.

Resolução:

O sal não é criado nem destruído no tanque, portanto qualquer

variação da quantidade de sal se deve unicamente aos fluxos de entrada e saída.

Como as vazões de entrada e saída são iguais, o volume de água no tanque

permanece constante em 100 litros; como a solução é homogênea, a concentração é

a mesma em todo o tanque e é dada por 100/)(rQ . Assim, a taxa com que o sal

deixa o tanque é 100/)(trQ Kg/minuto. Assim, a equação diferencial que descreve o

processo é 1004

rQr

dt

dQ−= . A condição inicial é 0)0( QQ = .

Podemos prever que a longo prazo a solução presente no tanque será

substituída pela solução que está entrando. Em conseqüência, a quantidade de sal

no tanque deve tender para 25 Kg.

A solução geral da equação diferencial é: 100exp25)( rt

ctQ−+= , onde c é uma constante arbitrária. Para satisfazer

a condição inicial devemos tomar 250 −= Qc , portanto a solução do problema de

valor inicial é 100

0 exp)25(25)( rtQtQ

−−+= , também

( ) 1000

100 expexp125)( rtrtQtQ

−− +−=

Supondo agora que 5020 == lQQ , então a equação

1000 exp)25(25)( rt

QtQ−−+= fica:

rtQ

03,0exp2525)( −+= .


23



UNESP - FEIS

Como 2% de 25 é 0,5, para encontrarmos o tempo T para que o valor

de Q(t) = 25,5, substituímos t = T e Q = 25,5 na equação

( ) 1000

100 expexp125)( rtrtQtQ

−− +−= , resolvendo para T temos:

.min4,13003,0

50ln≅=T

Agora para determinarmos r de forma que T = 45, voltamos a equação 100

0 exp)25(25)( rtQtQ

−−+= , igualamos t = 45, Q0 = 50, Q(t) = 25,5 e resolvemos em

r:

.min/69,850ln45

100galr ≅

=

Exemplo 10

Na investigação de um homicídio, ou de uma morte acidental, é muitas

vezes importante estimar o instante da morte. Vamos escrever uma forma

matemática para abordar este problema.

A partir de observações experimentais, sabe-se que a temperatura

superficial de um corpo se altera com uma taxa proporcional à diferença de

temperatura entre a do corpo e a temperatura ambiente. É o que se conhece como a

Lei de Newton do resfriamento. Assim, se θ(t) for a temperatura do corpo num

instante t, e se T for a temperatura constante do ambiente, então θ deve obedecer a

equação diferencial linear:

)( TKdt

d−−= θ

θ, onde K > 0 é a constante de proporcionalidade. O

sinal negativo desta equação provém do fato de se o corpo for mais quente que o

meio (θ > T) então ele se torna mais frio com o tempo.

Admitindo-se que no instante t = 0 descobre-se um cadáver e que sua

temperatura é medida e igual a θ0. Admitindo-se também que no momento da morte

tm a temperatura do corpo fosse θm (37º C).

A resolução da equação )( TKdt

d−−= θ

θ, com a condição inicial θ(0) =

θ0 é:

KtTTt

−−+= exp)()( 0θθ .


24



UNESP - FEIS

Fazendo uma Segunda medida de temperatura do corpo, num instante

t1, podemos determinar o valor de K. Suponhamos que θ = θ1. Substituindo-se estes

valores na equação KtTTt

−−+= exp)()( 0θθ , encontramos:

KtTT

−−=− exp)( 01 θθ , e daí

T

T

tK

−

−−=

0

1

1

ln1

θ

θ, onde θ0, θ1, T e t1 são grandezas conhecidas.

Para determinar tm, fazemos t = tm e θ = θm na equação

KtTTt

−−+= exp)()( 0θθ e resolvemos em tm, então temos:

T

T

Kt m

m−

−−=

0

ln1

θ

θ, onde K já é conhecida.

Agora como exemplo admitimos que a temperatura do corpo seja 30ºC

no instante da descoberta e 23ºC duas horas depois. A temperatura ambiente é

20ºC. Pela equação T

T

tK

−

−−=

0

1

1

ln1

θ

θ, temos:

2030

2023ln

2

1

−

−−=K

16020,0 −≅ h , pela equação T

T

Kt m

m−

−−=

0

ln1

θ

θ

.881,02030

2037ln

6020,0

1htd −≅

−

−−=

Portanto o corpo foi descoberto aproximadamente 53 minutos após a

morte.

Comandos para utilizar o Matlab (Utilizamos os micros contendo a

instalação do MATLAB no Departamento de Matemática da UNESP de Ilha Solteira,

para realização deste trabalho).

Para melhores detalhes deve-se consultar a Edição do Estudante do

MATLAB versão 5.

criação de vetores:

=>x=0:0.1:20; %cria um vetor x de dimensão 1x200

=>y=x.^2+2; %eleva cada elemento do vetor x ao quadrado (por isso o

ponto!)


25



UNESP - FEIS

e soma 2 a cada elemento

o comando plot:

=>plot(x,y)

=>title('gráfico da função y=x^{2}+2')% põe título no gráfico

=>xlabel('0 \leq x \leq 20 ')%rotula o eixo x

=>ylabel('y')%rotula o eixo y

Criação de variáveis simbólicas e o comando ezplot:

=>z='sin(x)'; %cria a variável simbólica z=sen(x) (por isso os plics)

=>w=diff(z); %calcula a derivada simbólica de z

=>f=diff(z,2) %calcula a deriva segunda de z

=>ezplot(z,[0 ,2*pi]) %faz a gráfico da função simbólica z

Comando grid on:

=>grid on %faz linhas de grade

Comando dsolve:

=>q=dsolve('Dq+3*q=t+exp(-2*t)') %encontra a solução geral da

equação diferencial y'+3y=t+exp(-2t)

=>r=dsolve('Dr+3*r=t+exp(-2*t)','r(0)=0')%encontra a solução da edo

que em t=0 vale 0

=>figure (2), ezplot(r,[0 , 10])

=>figure(2),plot(x,y,'r')%plota na cor vermelha (red)

Exemplo 1

Determine a solução do problema de valor inicial


26



UNESP - FEIS

y'-y/2=exp(-t)

y(0)=-1

y1=dsolve('Dy-y/2=exp(-t)','y(0)=-1');

y2=dsolve('Dy-y/2=exp(-t)','y(0)=1');

y3=dsolve('Dy-y/2=exp(-t)','y(0)=1/2');

y4=dsolve('Dy-y/2=exp(-t)','y(0)=0');

y5=dsolve('Dy-y/2=exp(-t)','y(0)=-1/2');

y6=dsolve('Dy-y/2=exp(-t)','y(0)=-2');

ezplot(y1)

hold on

ezplot(y2)

ezplot(y3)

ezplot(y4)

ezplot(y5)

ezplot(y6)

axis([0 5 -4 4])%Define a variação dos eixos

grid on

hold off

Exemplo 2

Achar a solução do problema de valor inicial

y'+2ty=t, y(0)=0

Plotar também o campo de direções para a equação.

y=dsolve('Dy+2*t*y=t','y(0)=0')

ezplot(y,[-1,4])

grid on

dfield5 %campo de direções

Dfield5, foi um programa encontrado na Internet. É um excelente

programa para visualizar campos de direções.


27



UNESP - FEIS

Equações de Variáveis Separáveis

Podemos usar x em vez de t para designar a variável independente de

uma equação diferencial, isto para ser mais conveniente. Neste caso, a equação

geral de primeira ordem assume a forma

),( yxfdx

dy=

Se esta equação é não-linear, ou seja, se f não é uma função linear de

variável dependente y, não existe um método geral para resolver a equação.

Primeiramente, reescrevemos a equação ),( yxfdx

dy= na forma

0),(),( =+dx

dyyxNyxM , podemos conseguir isso fazendo M(x,y) = -

f(x,y) e N(x,y) = 1. Caso M seja uma função apenas de x e N seja uma função

apenas de y, a equação 0),(),( =+dx

dyyxNyxM se torna

0)()( =+dx

dyyNxM , essa equação é dita separável pois pode ser

escrita na forma diferencial 0)()( =+ dyyNdxxM . Esta forma também é mais

simétrica e tende diminuir as diferenças entre as variáveis dependentes e

independentes.

Exemplo: Resolver o problema de valor inicial

)1(2

243 2

−

++=

y

xx

dx

dy, 1)0( −=y

Inicialmente escrevemos a equação na forma

dxxxdyy )243()1(2 2 ++=− , agora integrando o primeiro membro em

relação a y e o segundo membro em relação a x, temos:

cxxxyy +++=− 222 232 , de forma que c é uma constante arbitrária.

Fazendo x = 0 e y = -1, isto para satisfazer a condição inicial, temos que a constante

c = 3. Deste modo, uma solução para o problema de valor inicial é dado por

3222 232 +++=− xxxyy


28



UNESP - FEIS

O Teorema da existência da unicidade

Este teorema afirma que, sob certas condições de f(x,y), o problema de

valor inicial 00 )(),,(´ ytyytfy == , (1)

tem uma única solução num certo intervalo que contém o ponto t0.

Em alguns casos, a existência de um problema de valor inicial (1) pode

ser comprovada diretamente, pela resolução do problema e obtenção de uma

fórmula que dá a solução.

Em primeiro lugar, observamos que basta considerar o problema no

qual o ponto inicial (t0,y0) é a origem; isto é, o problema 0)0(),,(´ == yytfy .

Se o ponto inicial for outro, poderemos sempre fazer uma mudança

preliminar de variáveis correspondente a uma translação dos eixos do sistema de

coordenadas, e levar o ponto dado (t0,y0) para a origem. O teorema de existência e

unicidade pode, portanto, ser enunciado da forma:

Se f e yf ∂∂ / forem contínuas no domínio retangular

byatR ≤≤ ,: ,então há um intervalo aht ≤≤ no qual existe uma solução única

)(ty φ= do problema de valor inicial (2).

Suponhamos, que há uma função )(ty φ= que satisfaz o problema de

valor inicial; então [ ])(, ttf φ é uma função contínua exclusiva de t. Então podemos

integrar y´=f(t,y) do ponto inicial t=0, até um valor arbitrário de t, para ter:

[ ]dsssft

t

∫=0

)(,)( φφ (3), onde usamos a condição inicial 0)0( =φ .

Uma vez que (3) envolve a integral de uma função desconhecida f, denominada uma

equação integral. Esta equação não é uma fórmula para o problema de valor inicial,

mas proporciona outra relação satisfeita por qualquer solução da equação (2).

Inversamente, suponhamos que há uma função contínua )(ty φ= que satisfaz a

equação integral (3), o que mostra que a condição inicial é satisfeita. Além disso,

uma vez que o integrando da Eq.(3) é contínuo, pelo teorema fundamental do

cálculo, que [ ])(,)´( tftft =φ . Portanto, o problema de valor inicial e a equação

integral são equivalentes, no sentido de que qualquer solução de um é também

solução do outro. É mais conveniente mostrar que há uma única solução da


29



UNESP - FEIS

equação integral, num certo intervalo ht ≤ . A mesma conclusão valerá para o

problema de valor inicial.

Um método de mostrar que a equação integral (3) tem uma única

solução é conhecido como o método das aproximações sucessivas, ou método de

iteração de Picard. Ao se usar esse método, principia-se pela escolha de uma função

inicial f0, ou arbitrariamente ou da forma a se aproximar da solução do problema de

valor inicial. A escolha mais simples é 0)(0 =tφ ; e então 0φ satisfaz à condição inicial

da Eq.(2), embora, possivelmente, não satisfaça à equação diferencial. A

aproximação seguinte 1φ se obtem pela substituição de )(sφ por )(0 sφ no segundo

membro da Eq.(3), e simbolizando o resultado deste operação por )(1 tφ . Assim,

[ ]dsssft

t

∫=0

01 )(,)( φφ . Analogamente, 2φ se obtém de 1φ :

[ ]dsssft

t

∫=0

12 )(,)( φφ , e em geral,

[ ]dsssft

y

o

nn ∫=+ )(,)(1 φφ .

Desta forma, geramos uma seqüência de funções

{ } ,...,...,, 10 nn φφφφ = Cada membro desta seqüência satisfaz à condição inicial mas,

em geral, nenhuma satisfaz a equação diferencial.


30



UNESP - FEIS

Equações Lineares de Segunda Ordem

Oscilações Mecânicas e Oscilações Elétricas

Duas áreas significativas de aplicação são a das oscilações mecânicas

e a das oscilações elétricas. Por exemplo, o movimento de um corpo ligado a uma

mola, as oscilações de um eixo acoplado em um volante, a corrente elétrica num

circuito simples em série e muitos outros problemas físicos são descritos pela

solução de um problema de valor inicial da forma

´)0´(,)0(),(´´´ 00 yyyytgcybyay ===++ (1a)

Para se entender melhor, consideremos uma massa m suspensa na

extremidade de uma mola helicoidal vertical de comprimento em repouso l. A massa

provoca uma elongação L da mola, na direção para baixo (positiva). Duas forças

atuam no ponto onde a massa está ligada à mola. A força gravitacional, atua para

baixo e tem módulo mg. Há a força Fm devida a mola, que atua para cima. Admitindo

que a elongação da mola seja pequena, a força da mola é proporcional a L; relação

conhecida como Lei de Hooke. Escrevemos então Fm = -KL, onde K é a constante

da mola, e o sinal menos provém de a força da mola agir para cima. A massa

estando em equilíbrio, as duas forças se anulam, desta forma 0=− KLmg (2a).

Seja u(t) o deslocamento da massa em relação a posição do equilíbrio,

no instante t, medido positivamente para baixo. Então u(t) está relacionada às forças

que atuam sobre a massa pela lei do movimento, de Newton, )()´´( tftmu = (3a),

onde u´´ é a aceleração da massa e f é a força resultante que atua sobre ela. As

duas funções, u e f, são funções do tempo. Ao se determinar f, devem ser

consideradas quatro forças separadas:

• A força peso mg

• A força da mola Fm

• A força de amortecimento Fa

• Uma força externa F(t)

Levando em conta estas forças podemos escrever a lei de Newton (3a)

como

[ ] ).()´()()()()()´´( tFtutuLKmgtFtFtFmgtmu am +−+−=+++= γ


31



UNESP - FEIS

Uma vez que mg-KL=0, pela Eq.(2), deduz-se que a equação do

movimento da massa é

)()()´()´´( tFtKututmu =++ γ (4a)

A equação (4a) tem a mesma forma da equação (1a).

A fim de completar a formulação do problema de oscilação do corpo

precisamos especificar duas condições iniciais, a posição inicial u0 e a velocidade

inicial uo´ da massa:

u(0) = uo, u´(0) = uo´ (5a)

Oscilações Livres não-amortecidas

Na ausência de força externa, F(t) = 0 na Eq.(4a). Suponhamos que

não haja amortecimento, de modo que 0=γ . Esta é uma configuração idealizada do

sistema, que raramente é atingível na prática. Então a equação do movimento (4a)

se deduz a 0´´ =+Kumu (6a).

A solução da Eq.(6a) é tBtAu 00 sencos ωω += (7a),

onde m

K=2

0ω .

As constantes A e B podem ser determinadas se as condições iniciais

tiverem a forma de (5a). Ao se discutir a solução da Eq.(6a), podemos reescrevê-la

com o formato

)cos( 0 δω −= tRu , ou tRtRu 00 sensencoscos ωδωδ += (8a).

Comparando a Eq.(8a) e a Eq.(7a) vemos que A, B, R e δ estão

relacionadas pelas equações

δcosRA = , δsenRB = (9a)

Assim

22BAR += , .tan AB=δ (10a)

Ao se calcular δ é preciso cuidado para escolher o quadrante correto,

o que pode ser feito mediante a verificação dos sinais de cosδ e senδ nas Eq.(9a).

O gráfico da equação (8a) ou da equação equivalente (7a), descreve o

movimento periódico harmônico simples do corpo. O período do movimento é


32



UNESP - FEIS

2

1

0

22

==

K

mT π

ω

π (11a)

A freqüência circular mK=0ω , medida em radianos por unidade de

tempo, é a freqüência natural do sistema. O deslocamento R máximo do corpo em

relação à posição de equilíbrio é a amplitude do movimento. O parâmetro

adimensional δ é a fase, ou ângulo de fase, e mede o deslocamento da onda no

tempo, em relação à posição normal correspondente a 0=δ .

O movimento descrito pela equação (8a) tem uma amplitude constante,

que não diminui com o tempo. Este efeito reflete a ausência de amortecimento, pois

não há forma de o sistema dissipar a energia que recebeu no deslocamento e na

velocidade iniciais. Além disso, para uma dada massa, e uma dada constante da

mola, o sistema sempre oscila com a mesma freqüência 0ω , independente das

condições iniciais. No entanto, as condições iniciais ajudam a determinar a amplitude

do movimento. Pela Eq.(11a), observa-se que T aumenta quando m aumenta, e

assim corpos de massa maior oscilam mais lentamente. Por outro lado, T diminui

quando K aumenta, o que significa que as molas mais rígidas provocam oscilações

mais rápidas.

Oscilações Livres amortecidas

Incluindo o amortecimento a equação diferencial que governa o

amortecimento será

0´´´ =++ Kuumu γ (12a)

Estamos especialmente interessados em examinar os efeitos das

variações do coeficiente de amortecimento γ , dados os valores de massa m e da

constante da mola K. As raízes da equação característica são

−±−=

−±−=

2

2

21

411

22

4,

γ

γγγ Km

mm

Kmrr (13a)

Conforme o sinal de Km42 −γ , a solução u tem uma das seguintes

formas

)exp()exp(,04 212 trBtrAuKm +=>−γ

( )mtBAuKm 2exp)(,042 γγ −+==−


33



UNESP - FEIS

02

)4(),sencos(2exp(,04

2

12

2 >−

=+

−=<−

m

KmtBtA

mt

uKmγ

µµµγγ (14a)

O caso mais interessante é o terceiro, que ocorre quando o

amortecimento for pequeno. Se fizermos A = Rcosδ e B = Rsenδ na Eq.(14a),

obtemos

)cos(2Re δµγ −

−= t

mt

xpu .

O deslocamento u, num gráfico, fica entre as curvas

( )mtxpu 2Re γ−±= . A curva é parecida com uma co-senóide cuja amplitude diminui

à medida que t aumenta. O movimento é uma oscilação amortecida ou uma vibração

amortecida.

O parâmetro µ determina a freqüência de oscilação da massa; µ é

denominada quase-frequência. Comparando µ com a freqüência 0ω do movimento

não-amortecido, encontramos

KmKmmK

mKm

81

41

2/)4( 22

122

12

0

γγγ

ω

µ−≅

−=

−= (15a)

Não é a grandeza γ , apenas, que determina se o amortecimento é

grande ou pequeno, mas a grandeza de 2γ em comparação com 4Km. Quando

Km42γ for pequena, podemos desprezar o efeito do amortecimento no cálculo da

quase-frequência e do quase-período do movimento. Por outro lado, se quisermos

estudar o movimento detalhado do corpo em qualquer instante, nunca poderemos

desprezar a força de amortecimento, mesma que seja muito pequena.


34



UNESP - FEIS

Resultados

Pára-quedista

Problema 1

A partir de um exemplo de aplicação de EDO´s (exemplo 7), fizemos

modificações que achamos mais interessantes, então o problema ficou da seguinte

forma:

Suponha que o tempo que o pára-quedas leva para se abrir totalmente

seja de 5 segundos e que durante esse tempo a variação da resistência seja

descrita por uma equação linear do tipo K(t)= At + B. Admita que a pessoa abra o

pára-quedas 10 segundos após o salto de uma altura de 5000 pés. Calcule a

velocidade no momento em que o pára-quedas esteja totalmente aberto e a

distância que o pára-quedista percorreu neste intervalo de tempo. Sendo m = 5,625

slugs, K1 = 0,75, K2 = 12 e g = 32 ft/s2.

Resolução:

Tendo a equação da velocidade e do espaço em função da velocidade

gvm

tK

dt

dv=+

)( e )(tv

dt

dx=

A função K(t) =

≥

≤≤−

≤≤

15

15105

)(100

2

12

1

tseK

tseKK

tseK

Para :100 ≤≤ t

gm

vK

dt

dv

vm

vKg

dt

dv

=+⇒

=

−= 11

0)0(

O fator integrante neste caso é da forma

=

mtK

t 1exp)(µ

Multiplicando o fator integrante a todos os termos da equação, temos:


35



UNESP - FEIS

∫ ∫

=

⇒

=

+

t t

dtgm

tKv

mtK

dt

dg

mtK

vm

K

mtK

dt

dv

mtK

0 0

111111 expexpexpexpexp

−

−=⇒

−

=−

mtK

K

mg

K

mgtv

mtK

K

mgvtv

mtK 1

11

1

1

1 exp)(1exp)0()(exp

Substituindo os valores, obtemos: sftv /7,176)10( =

Para obtermos x(10) temos:

−

−

−−=⇒

−

−=∫ ∫ ∫ 1exp)(exp 1

1110 0

1

0 11m

tK

K

m

K

gm

K

gmtxdt

mtK

K

gmdt

K

gm

dt

dxt t t

Substituindo os valores, chegamos: x(10) = 1074,47 ft

Agora para o tempo 1510 ≤≤ t , temos:

gvm

tK

dt

dv=+

)(, calculando o fator integrante

= ∫

t

dsm

sKt

10

)(exp)(µ ,

para facilitar chamaremos a expressão entre parênteses de r(t). Como sabemos que

a multiplicação do fator integrante a todos os termos da equação diferencial nos dá a

derivada de um produto, podemos escrever:

)()´( tgv µµ = , agora integrando ambos os termos

( )∫ +=t

cdssrgttv10

)(exp)()( µ , onde c é uma constante de integração;

∫ ∫∫

−+

=

t tt

dsm

sKcdsdu

m

uK

t

gtv

10 1010

)(exp

)(exp

)()(

µ

Para resolvermos esta expressão, utilizamos o comando dsolve

(resolução simbólica), porém o Matlab nos dá a solução exata, mas não consegue

simplificar a expressão.

Passamos então para o método de resolução numérica, utilizando o

comando ODE 45 (utiliza o método de Runge-Kuta de quarta ordem).

Inicialmente indicamos uma função como arquivo de extensão do

Matlab;

Function dv = Resolve(t,y,flag,K1,ti,m)

K = (12 –0,75)*(t-10)/5 + 0,75;

dv = [-y(2); -(((12 – 0,75)*(t – 10)/5 + 0,75)/5,625)*y(2) + 32

A partir daí podíamos usar o comando ODE 45


36



UNESP - FEIS

[t0,y0] = ODE 45(´Resolve´, [10 15], [x(10); v(10)], [ ])

a1 = length (y1(:,1); % Este comando nos indica o tamanho do vetor,

com isso temos o número de linhas que o vetor possui.

X(15) = y1 (a1,1); % Para saber x(15) pedimos a última linha da primeira

coluna do vetor, tendo assim a posição no instante 15.

V(15) = y2 (a1,2); % Para saber v(15) pedimos a última linha da

segunda coluna do vetor, tendo assim a velocidade no instante 15.

Resta apenas plotar os gráficos para uma melhor visualização:

plot (t0,y0(:,2)), grid on % Com isso obteríamos o gráfico da velocidade,

o comando grid on nos auxilia em uma melhor visualização.

plot (t0,y0(:,1)), grid on % Com isso obteríamos o gráfico do espaço, o

comando grid on nos auxilia em uma melhor visualização.

Com o comando zoom on ativado obtemos um zoom que nos auxilia no

momento de verificar valores.


37



UNESP - FEIS

Em seguida resolvemos modificar novamente o problema:

Problema 2

Agora o problema consiste em determinar o instante em que o pára-

quedas deve ser aberto para que se chegue ao solo com uma determinada

velocidade.

Estipulamos uma velocidade de 20 ft/s, que através de verificações,

achamos que é uma velocidade razoável para que o pára-quedista não sofra

nenhuma lesão.

Resolução:

Para sabermos o instante de abertura simulamos o programa do

problema anterior para diversos tempos e com isso chegamos a um tempo de

abertura de 25,95 s.


38



UNESP - FEIS

Com o comando zoom on ativado, foi possível uma melhor

aproximação para que chegássemos a este resultado.

Aquecimento e Resfriamento em uma Sala

Artigo de C•ODE•E – 1995:

Nossa meta é formular e analisar um modelo matemático que descreva

a variação de temperatura de 24 horas dentro de uma sala. A variação interior será

o resultado da variação de temperatura externa e o calor gerado pelas pessoas e

máquinas dentro do edifício. Nós ignoramos o aquecimento e resfriamento do interior

com aquecedores ou ar-condicionado, assim a situação modelada se aplica melhor

em primavera ou outono.

Nós sabemos da Lei de Newton do Resfriamento que se:

a) T(t) = Temperatura dentro do edifício em função do tempo;

b) M(t) = Temperatura de fora (do ar) em função do tempo;

c) H(t) = Ganho de temperatura interno devido a pessoas, maquinaria,

luzes, etc., então

T’ = K(M(t) – T(t)) + H(t)

O fator K depende das propriedades físicas do edifício, como o número

de portas e janelas, e o tipo de isolação, mas não depende de T, M ou t.

Nós assumiremos que a temperatura do ar varie de uma forma

senoidal em um período de 24 horas, supondo o mínimo em t = 0 (meia-noite) e o

máximo em t =12 (meio-dia); de forma que

M(t) = M0 – Bcoswt

onde B é uma constante positiva e

w = 2π/24 = π/12.

Nós também assumiremos que H(t)=H0, de forma que a taxa interna do

ganho de temperatura é constante. Supondo que o ganho de calor devido a luzes e

maquinaria em função do tempo é uma constante, uma suposição razoável.

1) Manualmente, mostre que a equação pode ser escrita na forma de

equação linear:

T’ + KT = Q(t), onde K é uma constante e

Q(t) = K(B0 – Bcoswt)


39



UNESP - FEIS

Note que B0 é outra constante, derivada como uma combinação de

M0,H0 e K.

2) Assuma isso à meia-noite, 21 de setembro (o começo de queda), o

ar-condicionado é desligado e a temperatura do edifício é 22,22º C, assim as

hipóteses feitas na dedução do modelo podem ser aplicadas. Supondo que durante

os próximos 3 dias, a temperatura varie da meia-noite ao meio-dia de 10º C a 32,22º

C de forma senoidal suposta acima. Supondo a constante K = 1/6. Considerando H0

=1 para um ganho de calor modesto. Devem ser determinados M0 e B.

Determinando estes a partir das informações dadas. Determine a equação para o

edifício a partir destas circunstâncias. Escreva o problema de valor inicial

correspondente com todos os parâmetros convertidos aos valores numéricos

apropriados.

3) Usando o comando dsolve do Matlab resolver a equação de

primeira-ordem linear, tente resolver a equação diferencial que usa esta " fórmula "

automática fácil. Especule em por que Matlab devolve a resposta que faz. Você pode

usar esta resposta?

4) Use Matlab para resolver a equação diferencial.

a) Escreva a solução em forma matemática.

b) Faça o gráfico da temperatura do edifício e a temperatura do ar em

um período de 3 dias (72 horas).

c) Inspecionando seu gráfico, ache para cada um dos 3 dias o máximo

e o mínimo de temperaturas no edifício e o tempo que eles acontecem.

5) Explique a diferença de tempo entre as variações de temperatura

externas e internas.

6) Baseado em sua solução de parte (4a), você poderia ter predito os

resultados de parte (5) sem referência para os gráficos? Examinando sua solução de

(4a), que parâmetros deveria variar para se fazer o edifício mais confortável para os

ocupantes. Explique suas recomendações e o resultado gráfico.

Resolução:

Temos que: ( ) )()()()´( tHtTtMKtT +−= , então:


40



UNESP - FEIS

)()()()´( tHtKMtKTtT ++−= , onde )()()( tftHtKM =+ ; chamaremos

agora ´)´( ytT = , então desta forma teremos:

)(´ tfKyy =+ , onde K = constante.

Temos que: função seno [0, 2π], queremos agora transformar em uma

outra função com outra amplitude e período [0, 24], isto para um dia.

O termo que multiplica [0, 24] que nos dá a função propriamente dita é

12

π, então com isso teremos

t

12sen

πe agora chamando

12

πω = , teremos tωsen .

Como a hora mais quente do dia é às 12 horas e a mais fria são às 24

horas deslocamos a função, assim:

( ) ( ) [ ]ttBtBtB ωωωϖωωω cos6sen6cossen6sen)6(sen −=−=− , agora

substituindo o valor de ω e simplificando, temos

[ ] tBtB ωω coscos −=−

Dependendo da temperatura desejada somamos B0 a equação:

( ) ( )( )K

HMBtBBKtftBMtM 0

0000 cos)(cos)( +=⇒−=⇒−= ωω

Temos que a temperatura às 24 hs é de 10º C , que a temperatura às

12 hs é de 32,22º C e que a temperatura externa da sala é de CT º22,22)0( = e

também sabendo o valor das constantes 61=K e 10 =H , podemos:

BMBMM +=⇒+= 00 22,32)12( (1)

BMBMM −=⇒

−= 00 1024

12cos)24(

π (2)

Resolvendo os sistemas (1) e (2), temos: CM º11,210 = e CB º11,11= ,

agora a partir destes podemos ter 11,27

611

11,21 00 =⇒+= BB .

Agora podemos resolver utilizando Matlab:

Criando um arquivo de função:

function dt=resolve(t,y,flag,B0,B,K)

dt=[-K*y+K*(B0-B*cos(pi*(t-2)/12))];

Utilizando agora o comando ODE 45 (resolução numérica):


41



UNESP - FEIS

x=0:0.1:24;

K=1/6;

H0=1;

M_max=32.22;

M_min=10;

M0=(M_max +M_min)/2;

B=(M_max -M_min)/2;

B0=M0+H0/K;

%T=dsolve('DT+K*T=K*(B0-B*cos((pi*t)/12))','T(0)=22.22')

M=M0-B*cos(pi*(x-2)/12);

[t,T]=ode45('resolve',[0 24],22.22,[],B0,B,K);

Depois de resolvido vamos agora plotar dois gráficos em uma única

janela:

Subplot (2, 1, 1), plot (t, T), grid on, zoom on

Subplot (2, 1, 2), plot (x, M), grid on, zoom on

Com esses comandos obtemos os gráficos:


42



UNESP - FEIS

Resolvemos modificar o problema da seguinte forma: A constante H0

corresponde ao aquecimento interno da sala devido as pessoas e máquinas, então

estipulamos horários em que estas pessoas não estariam na sala e com isso a

constante Ho = 0. Os horários estipulados foram os seguintes:

=≤≤

=≤≤

=≤≤

=≤≤

=≤≤

02418

11814

01412

1128

080

0

0

0

0

0

Ht

Ht

Ht

Ht

Ht

A temperatura máxima às 12 hs e a mínima às 24 hs seriam

respectivamente de 36º e 25º C e a temperatura externa da sala de 26º C.

Inicialmente criamos os arquivos de funções:

function dt=resolve(t,y,flag,B0,B,K)

dt=[-K*y+K*(B0-B*cos(pi*(t-2)/12))];

function dw=resolve1(t,y,flag,B1,B,K)

dw=[-K*y+K*(B1-B*cos(pi*(t-2)/12))];

Agora podemos utilizar o comando ode45:

x=0:0.1:24;

K=1/6; % fluxo de calor entre o meio e o edifício

H0=0;

H1=1; % produção de calor no interior

M_max=36;

temp=26

M_min=25;

M0=(M_max+M_min)/2;

B=(M_max-M_min)/2;

B0=M0+H0/K;

B1=M0+H1/K;

M=M0-B*cos(pi*(x-2)/12);

[t0,T0]=ode45('resolve',[0 8],temp,[],B0,B,K);


43



UNESP - FEIS

a0=length(T0(:,1));

temp0=T0(a0,1);

[t1,T1]=ode45('resolve1',[8 12],temp0,[],B1,B,K);

a1=length(T1(:,1));

temp1=T1(a1,1);

[t2,T2]=ode45('resolve',[12 14],temp1,[],B0,B,K);

a2=length(T2(:,1));

temp2=T2(a2,1);

[t3,T3]=ode45('resolve1',[14 18],temp2,[],B1,B,K);

a3=length(T3(:,1));

temp3=T3(a3,1);

[t4,T4]=ode45('resolve',[18 24],temp3,[],B0,B,K);

a4=length(T4(:,1));

temp_final=T4(a4,1)

plot(t0,T0(:,1),t1,T1(:,1),t2,T2(:,1),t3,T3(:,1),t4,T4(:,1)), grid on

zoom on

A partir desse comando executado obtemos o seguinte gráfico:


44



UNESP - FEIS

Quanto tempo leva para um oscilador harmônico parar de oscilar?

Artigo de C•ODE•E – 1995:

Eu pedi para meus estudantes que determinassem em uma nova

classe de equações diferenciais, quanto tempo que um oscilador linear com uma

certa massa leva para voltar ao repouso depois que entrasse em movimento por um

certo deslocamento de inicial. A resposta que nós achamos nos surpreendeu.

Supondo que o retorno ao descanso é somente assintótico então

levaria o tempo ao infinito, nós concordamos que nossos instrumentos não eram

capazes de medir qualquer deslocamento menor que 0,01. As unidades atuais de

medida permaneceram não especificadas. Isto significa que uma vez que o oscilador

está dentro de um deslocamento de 0,01 de sua posição de equilíbrio, nós podemos

considerar que é em repouso. O tempo para isto acontecer será agora finito. Depois

de alguma discussão adicional aos fatores que influenciam o tempo para devolver o

oscilador ao repouso, nós decidimos enfocar na constante de amortecimento. Nós

formulamos por seguir perguntas específicas. Dado o problema de valor inicial:

0)0(',1)0(,02

2

===++ xxxdt

dxb

dt

xd

1) Determine o valor de amortecimento crítico por este oscilador. Faça

a série de tempo de deslocamento (x por t) para este valor de b, e de uma estimativa

de quanto tempo leva para o oscilador retornar ao descanso, como nós definimos.

2) Que valor da constante de amortecimento b minimiza o tempo para

fazer o oscilador retornar ao repouso?

3) Para um determinado valor de b, deixe T(b) = o tempo para retornar

ao repouso. Esboce um gráfico da função T no intervalo 0 <b <2.

A questão 1 é projetada para adquirir estudantes que pensam ao longo

de uma direção viável. O amortecimento crítico acontece quando b é 2, e resolvendo

a equação diferencial para este b, estime a função de deslocamento ttteetx

−− +=)( .

Uma calculadora gráfica pode ser usada para plotar esta função e resolver a

equação =)(tx 0,01. Nós achamos que o tempo para o oscilador retornar ao repouso

é aproximadamente de 6,638.

A questão 2 pode ser respondida de dois modos. A maioria de meus

estudantes eram rápidos para perceber aquele aumentando de amortecimento sobre


45



UNESP - FEIS

o valor crítico e faria o oscilador se orientar mais lentamente a equilíbrio de forma

que nós deveríamos procurar minimizar o tempo em algum lugar a baixo do valor

crítico 2. Considerando que nós estávamos usando o programa MDEP que nos

permitiu plotar convenientemente, sugeri que nós começássemos a procurar

graficamente através da plotagem do deslocamento x(t) para vários valores de b.

Nós poderíamos examinar cada instante, então determinaríamos quando a curva

entraria no eixo horizontal saltada primeiro por x = ±0,01 e permaneceria lá. A Figura

1 mostra a série de tempo para b = 1,5, 1,6, e 1,7. Deste gráfico nós podemos

calcular T(1,5) = 6,6, T(1,6) = 6,4, e T(1,7) = 4,5. O gráfico também mostra que os

valores corretos estão entre b = 1,6 e b = 1,7.

Figura 1: Série de intervalos para valores de b mostrarem diferentes T

para que o oscilador retorne ao repouso.

Nós podemos estreitar nossa procura agora para este intervalo e

podemos continuar plotando. Esta aproximação feita nos permite calcular o tempo

mínimo de amortecimento constante para duas casas decimais ou até mesmo três

se nós formos persistentes. Porém, há uma alternativa. Se nós olharmos o gráfico,

talvez depois de somar algumas curvas, nós podemos ver cuidadosamente que o

valor de b o qual olhamos ocorre quando o intervalo de tempo é tangente à linha x =

-0,01. Veja Figura 2.

Figura 2: Somatória de curvas para obter b onde T para retornar ao

repouso é mínimo

Isto nos dá um modo para achar b analiticamente. Nós precisamos

achar o valor de b para o qual dx/dt = 0 e x(t) = -0,01. Para estes próximos passos

um sistema de cálculo de computador será extremamente útil. Resolvendo nosso

sistema original analiticamente nós adquirimos

−+

−

−

=−

−

tbeb

tbsinbe

tx

bt

bt

22

2

22

42

1cos

4

42

1

)(

e então

2

22

4

42

1

2b

tbsine

dt

dx

bt

−

−

−=

−


46



UNESP - FEIS

Se nós fixássemos dx/dt = 0 e resolver para t, nós adquirimos t

= ( )24/2 bK −π . Desde que nós queiramos dx/dt = 0 e x(t) = -0,01, nós levamos K =

1. Se nós substituirmos este valor de t na expressão para x(t), nós adquirimos x(t)

= ( )24/exp bb −−− π . Agora iguale a -0,01 e resolve em b para adquirir b = 1,65217.

De fato para achar o tempo por este valor de b, nós temos que resolver a equação

x(t) = 0,01 para t, porque este é o tempo quando a curva se desloca entre o eixo

horizontal saltada por x =± 0,01. O tempo mínimo é 4,19161.

A questão 3 é a mais intrigante e difícil. Novamente pode se tentar

achar graficamente o valor de T(b) para valores de b dentro do intervalo 0 <b <2.

Porém, a menos que nós escolhamos um grid bom de valores de b, este método

negligenciará o fato provavelmente T é uma função descontínua! Realmente a

função T tem uma descontinuidade de salto a cada valor de b onde o deslocamento

x(t) é tangente a linha horizontal x = ±0,01. A Figura 2 provê a perspicácia da que

nós precisamos. Como os aumentos de b de 1,5 para 1,65217 é possível ver que o

tempo para retornar ao repouso é determinado pela interseção da curva de

deslocamento com a linha horizontal x = -0,01. Esta interseção move-se

continuamente à esquerda. Porém, ao ponto b = 1,65217, o tempo salta e é agora

determinado pela interseção com a linha horizontal x = 0,01. Podem ser achados

estes pontos de descontinuidade da mesma maneira que o b mínimo foi achado na

questão 2. A fórmula é

,...,2,1,

)01,0ln(1

22

=

+

= K

K

bK

π

e os primeiros valores são 1,65217, 1,18231, 0,878036.

Plotando a função T avaliamos T a cada bK e então avaliamos uma

quantidade de pontos entre estes valores. O resultado é o gráfico mostrado abaixo.

Foi criado utilizando MATLAB.

Figura 3: Descontinuidades de tempo para o oscilador retornar ao

repouso, plotados em b.

Resolução:

Através da resolução numérica, criamos um arquivo de função:

function dw=pendulo(t,y,flag,b)

dw=[y(2); -b*y(2)-y(1)];


47



UNESP - FEIS

Agora utilizando o comando ode45, resolvemos:

[t0,y0]=ode45('pendulo',[0 20],[1;0],[],1.5);

[t1,y1]=ode45('pendulo',[0 20],[1;0],[],1.6);

[t2,y2]=ode45('pendulo',[0 20],[1;0],[],1.7);

plot(t0,y0(:,1),'b',t1,y1(:,1),'r',t2,y2(:,1),'g'),axis([4 12 -0.03 0.01]),grid on,

xlabel('t'), ylabel('x(t)'), title('posições para b=1,5, 1,6 e 1,7')

A partir deste comando executado, obtemos o seguinte gráfico para as

diferentes posições de b:


48



UNESP - FEIS

Conclusões

Através da realização desta Iniciação Científica foi possível

compreender o poder e o limite dos métodos matemáticos utilizados, além de que

tais métodos serviram como um primeiro passo para uma certa formação

matemática, a qual me deu confiança para formular e explorar novos modelos.

Este trabalho também me propiciou uma experiência inicial com a

teoria de equações diferenciais ordinárias, modelagem matemática e a utilização do

Matlab o qual nos tornou possível realizar simulações e visualizações de soluções.

As discussões com o professor orientador foram de suma importância,

pois os modelos matemáticos apresentados neste trabalho foram formulados e

explorados através destas discussões.


49



UNESP - FEIS

Publicações originadas a partir do presente trabalho

FOGAL, Marcelo Luiz Freitas; OLIVEIRA, Ernandes R. Equações

Diferenciais Ordinárias e Aplicações. In: VII REUNIÃO DE INICIAÇÃO CIENTÍFICA

DA FEIS, 1999, Ilha Solteira. 1999.

FOGAL, Marcelo Luiz Freitas; OLIVEIRA, Ernandes R. Simulação de

um Salto com Pára-quedas. In: IX REUNIÃO DE INICIAÇÃO CIENTÍFICA DA FEIS,

2001, Ilha Solteira. 1999.

FOGAL, Marcelo Luiz Freitas; OLIVEIRA, Ernandes R. Simulação de

um Salto com Pára-quedas. In: XIII CONGRESSO DE INICIAÇÃO CIENTÍFICA DA

UNESP, 2001, Bauru. 1999


50



UNESP - FEIS

Referências bibliográficas

Willian E. Boyce and Richard C. DiPrima, Equações Diferenciais Elementares e

Problemas de Valores de Contorno, LTC – Livros Técnicos e Científicos, 1998.

Rodney Carlos Bassanezi and Wilson Castro Ferreira, Equações Diferenciais com

Aplicações, HARBRA, 1998.

Consortium for Ordinary Differential Equations Experiments – C. ODE.E., 1995.

The Student Edition of MATLAB - Vesion 5 User's Guide - The Mathworks Inc.

Documents

Eq Diferenciais Ordinarias