1

2. EQUAÇÕES DE DIFERENÇAS E EQUAÇÕES DIFERENCIAIS ORDINÁRIAS

Equações de diferenças e equações diferenciais ordinárias representam uma versátil ferramenta para a análise de sistemas dinâmicos,

Possuem uma grande riqueza teórica e uma excelente representação de muitos

sistemas dinâmicos,

Por enquanto, será abordado apenas o caso uni-variável, entretanto, a teoria associada fornece subsídios básicos para uma teoria mais geral envolvendo sistemas multi-variáveis.

2.1 Equações de diferenças Suponha uma variável y(k) de um sistema dinâmico representado uma sequencia de pontos representados de maneira discreta e equiespaçadas no tempo indexado por k. Esta variável y(k) representa um número real associado com cada um destes pontos. Uma equação de diferenças permite a relação do valor y(k), no instante k, com outros pontos, usualmente na vizinhança de y(k). Exemplos de equações de diferenças:

)()1( kayky ,2,1,0k (2.1)

)1()(21)1()2( kykykaykky ,2,1,0k (2.2)

A ordem de uma equação de diferenças é definida pela diferença entre o maior e o menor índice que aparecem na equação. Uma equação de diferenças linear possui a seguinte forma,

)()()()1()()1()()()( 011 kgkykakykankykankyka nn (2.3) onde os termos skan

')( são os coeficientes (ou parâmetros) da equação a diferenças. Se eles não dependem de k, são chamados invariantes no tempo e )(kg denota o termo força ou simplesmente lado direito. Soluções Uma solução de uma equação de diferenças é uma função )(ky que reduz a equação (2.3) a uma identidade. Exemplo 1- Seja a seguinte equação de diferenças )()1( kayky para

2,1,0k , a solução desta equação é,

kCaky )(

2

onde C é uma constante. Para mostrar que kCaky )( é realmente a solução da equação acima, basta substituir como,

11 kkk CaaCaCa (identidade) Exemplo 2 – Seja a seguinte equação de diferenças 1)()1()1( kkykyk para

2,1,0k , a solução para tal equação é a seguinte,

kAky 1)(

sendo A uma constante. Basta substituir a solução na equação a diferenças acima para comprovar a identidade.

Exemplo 3 – A equação de diferenças não linear )(1

)()1(ky

kyky

possui a seguinte

solução,

AkAky

1

)(

sendo A uma constante. Substitua e comprove que a solução acima realmente procede. Exemplo 4 – A equação a diferenças não linear 1)()1( 22 kyky por sua vez não possui solução. 2.2 Teoremas da Existência e Unicidade das Soluções

Como com qualquer conjunto de equações, uma equação de diferenças não necessariamente pode possuir uma solução, e se ela possuir, tal solução pode não ser única. Será abordado agora a um exame geral das questões da existência e unicidade das soluções das equações de diferenças. Condições Iniciais: representa uma característica essencial de uma equação de diferenças. Exemplo: Seja a seguinte equação de diferenças )(2)1( kyky . Então, para 1,0k

)1(2)2()0(2)1(

yyyy

tem-se duas equações para três incógnitas, )0(y )1(y e )2(y .

Agora para 2,1,0k tem-se,

3

)2(2)3()1(2)2()0(2)1(

yyyyyy

tem-se agora três equações para quatro incógnitas.

De maneira geral, existe n incógnitas a mais que equações, sendo n a ordem da equação de diferenças. Para especificar os primeiros n valores de y(k), são usadas as condições iniciais )0(y , )1(y , ..., )1( ny . Exemplo5: Equação de diferenças linear de primeira ordem.

)()1( kayky A solução da equação acima é kaCky )( , impondo a condição inicial

0)0( yy para 0k na solução, implica que Cy 0 , logo,

kayky 0)( Exemplo 6: Mostre que a equação a diferenças linear de segunda ordem )()2( kyky com condições iniciais )0(y e )1(y possui a seguinte solução,

2

)1()0()1(2

)1()0()( yyyyky k

Teorema da Existência e Unicidade das Soluções Seja uma equação de diferenças da forma,

0),(,),1()( kkynkyfnky (2.4) onde f é uma função real arbitrária, definida sobre uma seqüência finita ou infinita de valores consecutivos de k. A equação possui uma e apenas uma solução correspondente para cada especificação arbitrária de n condições iniciais )( 0ky , )1( 0 ky ,...,

)1( 0 nky . 2.3 Equações de Diferenças de Primeira Ordem Seja a equação de diferenças de primeira ordem,

bkayky )()1( (2.5) onde a é um coeficiente constante e b é chamado termo força. A equação acima pode ser usada em muitas situações como: juros e amortizações, modelos populacionais, modelo de oferta e procura, etc. Seja a seguinte condição inicial Cy )0( . Vai-se usar o método recursivo para obter a solução de (2.5). Substituindo diversos valores de k em (2.5), obtém o seguinte,

4

Cy )0(

baCbayy )0()1( babCabayy 2)1()2(

babbaCabayy 23)2()3( Assim, continuando com o procedimento acima, obtém-se a forma geral da solução de (2.5) como,

baaaCaky kkk 1)( 21 (2.6)

Mas a

aaaak

kk

11121 , logo a equação acima se reduz a,

baaCaky

kk

1

1)(

Para 1a implica que kbCky )( . Para 1a implica que baaCaky

kk

11)( .

2.4 Modelo de Oferta e Procura Trata-se de um modelo de oferta e procura. Seja d uma variável que representa a demanda, s a oferta e p o preço de certa mercadoria.

A dependência da demanda com o preço d(p), do ponto de vista do consumidor é a seguinte,

apdpd 0)( (2.7) onde 0d e a são constantes maiores que zero. Isto é, a quantidade de mercadorias comprada pelos consumidores decresce quando o preço aumenta. A dependência da oferta com o preço s(p), do ponto de vista do produtor, pode ser modelada pela como,

bpsps 0)( (2.8) onde 0s é uma constante positiva ou negativa e b uma constante positiva. Isto é, a oferta aumenta com o preço pelo produtor. A Figura 1 mostra um gráfico típico de oferta e procura usando as equações (2.7) e (2.8). Neste caso foram usados os seguintes parâmetros 100000 d , 120a ,

1000 s , 70b , sendo que o preço máximo do bem é assumido igual a 100 unidades financeiras.

A condição de equilíbrio sugere que a oferta seja igual a procura. O preço nesta condição de equilíbrio é chamado preço de equilíbrio.

Seja que no instante k exista um preço p(k). O produtor baseia sua produção futura no instante k+1 como,

5

)()1( 0 kbpsks

A procura no instante k+1 é por sua vez o seguinte,

)1()1( 0 kapdkd

Impondo a condição de equilíbrio para o instante k+1, isto é )1()1( kskd , obtém o seguinte,

)1()( 00 kapdkbps

ou

asdkp

abkp 00)()1(

(2.9)

Figura 1 – Gráfico da oferta e procura de um bem. O preço de equilíbrio é obtido impondo a seguinte condição )1()( kpkp , que resulta em,

basdkp

00)( (2.10)

A solução geral para (2.9) é a seguinte,

)(1)0()( 00 sdbaabp

abkp

kk

(2.11)

6

Se ab , para k , .)( eqpkp . A Figura 2 mostra uma situação onde o preço converge para um valor de equilíbrio usando a Equação (2.11), para uma condição inicial 90)0( p unidades financeiras. Os dados para a realização desta simulação foram os mesmos usados para gerar as curvas de oferta e procura da Figura (1). O código usado (em Matlab) é apresentado também. Compare que o preço de equilíbrio pode ser verificado também pela Figura (1).

Por outro lado, se ab , )(kp diverge.

Figura 2 – Evolução temporal de preço.

clear all clc d0=input('de o valor de d0 procura inicial '); a=input('de o valor do parametro a '); s0=input('de o valor de s0 oferta inicial '); b=input('de o valor do parametro b '); pmax=input('de o valor do preço maximo '); p0=input('de o valor do preço inicial do bem '); N=input('de o nuero de pontos para os graficos de oferta e procura '); dp=pmax/N; for i=1:N p(i)=dp*(i-1); d(i)=d0-a*p(i); s(i)=s0+b*p(i); end Np=input('de o numero de pontos para curva do preço de equilibrio '); for k=1:Np P(k)=(-b/a)^k*p0+((1-(-b/a)^k)/(a+b))*(d0-s0); index(k)=k; end figure(1) plot(p,d,p,s)

7

xlabel('preço') ylabel('oferta e procura') legend('demenda','oferta') grid figure(2) plot(index,P) xlabel('instante k') ylabel('preço') grid 2.5 Equações de Diferenças Lineares Até aqui foi considerado apenas o método recursivo para resolver as equações de diferenças. Este assunto será tratado de maneira mais sistemática a partir deste ponto. Uma equação de diferenças linear com parâmetros variantes no tempo possui a seguinte forma,

)()()()1()()1()()()( 011 kgkykakykankykankyka nn (2.12)

Expressões analíticas para as soluções, Mesmo quando tais soluções não são obtidas, a teoria de equações de diferenças

lineares fornece importantes informações estruturais. Tais informações serão apresentadas em termos de teoremas que seguem.

A Equação Homogênea Se o termo 0)( kg para qualquer instante k, então a equação (2.12) é chamada homogênea. Se 0)( kg para algum k, então (2.12) não é homogênea. A equação de diferenças homogênea associada a (2.12) é, portanto,

0)()()1()()1()()()( 011 kykakykankykankyka nn (2.13) A equação homogênea desempenha um papel fundamental na especificação das soluções para a equação não homogênea. Teorema 1: Seja )(ky uma dada solução para equação (2.12). Então a coleção de todas as soluções desta equação é a coleção de todas as funções da forma, )()()( kzkyky onde z(k) é uma solução da correspondente homogênea associada com (2.13). Demonstração do teorema, vide Luenberger. Exemplo 7: Considere a equação de diferenças linear com parâmetros variantes com o tempo,

1)()1()1( kkykyk 1k

Por inspeção uma solução para esta equação de diferenças é 1)( ky .

8

A homogênea associada é obtida anulando o lado direito (ou termo força) como,

0)()1()1( kkzkzk cuja solução possui a seguinte forma,

kAkz )( Assim, com base no Teorema 1, a solução geral (ou solução total) é obtida como,

kAky 1)( Obs.: o valor da variável A depende da escolha da condição inicial. Teorema 2: Se )(1 kz , )(,),(2 kzkz m são todas soluções da homogênea (2.13), então qualquer combinação linear,

)()()()( 2211 kzckzckzckz mm onde mccc ,,, 21 são constantes arbitrarias, é também uma solução da homogênea (2.13). Conjunto Fundamental de Soluções Seja a equação homogênea (2.13). Um conjunto fundamental de soluções é qualquer conjunto de n soluções de (2.13) linearmente independentes. Teorema 3: Se z(k) é qualquer solução da homogênea (2.13), então z(k) pode ser expressa em termos de um conjunto fundamental destas n soluções na forma,

)()()()( 2211 kzckzckzckz nn para as constantes mccc ,,, 21 e )()(),( 21 kzkzkz n sendo as n soluções linearmente independentes. Exemplo 8: Considere a seguinte equação de diferenças homogênea de segunda ordem,

0)()1(2)2( kzkzkz

Por inspeção, 1)(1 kz e kkz )(2 são soluções para esta equação de diferenças, porém elas não representam um conjunto fundamental de soluções, pois não são linearmente independentes.

Entretanto, kkz 1)(1 e kkz )(2 representa um conjunto fundamental de soluções, pois elas são linearmente independentes. Logo, uma solução arbitrária do problema acima pode ser representada como,

9

sendo que 1c e 2c são obtidos pelas condições iniciais impostas pelo analista. Teorema 4: Suponha )()(),( 21 kzkzkz n um conjunto de soluções linearmente independentes para a equação homogênea (2.13). Então, qualquer solução z(k) de (2.13) pode ser expressa como combinação linear,

)()()()( 2211 kzckzckzckz nn sendo que nccc ,,, 21 são constantes que podem ser obtidos pelas condições iniciais. Solução da Não Homogênea Para encontrar a solução da equação não homogênea (2.12) procede-se da seguinte forma cumprindo os dois seguintes passos:

Encontre um conjunto de n soluções linearmente independentes para a equação homogênea associada e

Encontre uma solução particular para a equação não homogênea (2.12) que não necessariamente satisfaça as condições iniciais dadas. Esta solução particular ‘se parece muito com o lado direito (ou termo força) de (2.12)’.

A solução geral para a equação (2.12) é então modificada pela adição da

combinação linear de soluções homogêneas. Se )(ky é solução particular de (2.12) e )()(),( 21 kzkzkz n é um conjunto de

soluções da homogênea associada linearmente independentes, então a solução geral para (2.12) é,

)()()()()( 2211 kzckzckzckyky nn

Observaçoes:

Métodos recursivos podem ser usados para encontrar a solução particular, Métodos analíticos para determinação das soluções particulares são disponíveis

apenas para casos especiais. 2.6 Equações de Diferenças Lineares com Coeficientes Constantes

É possível encontrar todas a n soluções LI para a equação homogênea. A Equação Característica Considere a equação homogênea de diferenças linear com coeficientes constantes,

0)()1()1()( 011 kyakyankyanky n (2.14)

kccckckckz )()1()( 21121

10

A forma da solução adotada para (2.14) é a seguinte,

kkz )( (2.15) sendo uma constante a ser determinada. Substituindo (2.15) em (2.14) dá o seguinte,

001

11

1

kknkn

nk aaa (2.16) Multiplicando (2.16) por k obtém-se o seguinte polinômio,

0011

1 aaa n

nn (2.17)

O polinômio acima é chamado de polinômio característico ou equação característica. Note que um conjunto de n raízes distintas do polinômio característico leva a um conjunto de n soluções linearmente independentes do tipo k

ii kz )( , como sugeridas no Teorema 4. Teorema 5: Uma condição necessária e suficiente para que a seqüência geométrica

kkz )( seja solução de (2.14) é que a constante satisfaça a equação característica (2.17). Observações: 1) a equação característica é um polinômio de grau n, 2) as raízes do polinômio característico podem ser todas distintas ou repetidas, nestas notas é considerado apenas o caso de raízes distintas e 3) tais raízes distintas implica em um conjunto LI de soluções e 4) as raízes da equação característica podem ser real ou complexas. Exemplo 8: Equação de primeira ordem.

)()1( kyaky

aa 0

kCaky )( (solução) sendo que a constante C é determinada pela condição inicial 0)0( yy Exemplo 9: Equação de primeira ordem não homogênea.

bkayky )()1( Como uma solução particular vai-se tentar Cky )( . Substituindo na equação acima dá o seguinte,

abCbCaC

1

, portanto, a solução geral é a seguinte,

11

abCaky k

1)(

sendo que o termo C é determinado usando a condição inicial 0)0( yy . Exemplo 10: Resolva a equação de segunda ordem com termo força

kkykyky 3)(2)1(3)2( para as condições iniciais 0)0( yy e 1)1( yy . Solução homogênea

21023 212 e

kkk cccckz 221)( 2121

Solução particular Adotando-se a solução particular do tipo kCky 3)( e substituindo na equação a diferenças de segunda ordem acima dá o seguinte,

21332333 12 CC kkkk , portanto, a solução geral é a seguinte,

kkcckykzky 3

212)()()( 21

Impondo as condições iniciais

Para a condição inicial 0)0( yy , implica que 021 21 ycc

Para a condição inicial 1)1( yy , implica que 121 232 ycc

Que leva ao seguinte sistema de equações algébricas lineares para a determinação simultânea das constantes 1c e 2c ,

2321

2111

1

0

2

1

y

y

cc

.

Respondas o seguinte com base no exercício anterior:1) o sistema é estável? 2)

se as duas condições iniciais forem nulas, isto é, 0)0( y e 0)1( y , a resposta da equação de diferenças acima terá a parte homogênea como parte integrante da resposta total? e 3) quais as condições iniciais que levam a não existência da parte homogênea da resposta total?

Exemplo 11: Seqüência de Fibonacci.

12

Considere a série de números 1, 1, 2, 3, 5, 8, 13, 21, ... .Tal seqüência é conhecida como seqüência de Fibinacci. Seus termos são gerados através equação de diferenças homogênea )()1()2( kykyky com as condições iniciais

.1)2()1( yy Solução

251

25101 21

2

e

Condições iniciais

11

222

211

2211

cccc

que implica 5

15

121 cec

51

251

251)(

kk

ky

Exemplo 12: Resolva a seguinte equação de diferenças 0)()2( kyky com

0)1(1)0( yey .

1012 i , portanto,

kk icicky )()()( 21 Condições iniciais

0)()(1

21

21

iciccc

implica que 2121 cc

O que leva a seguinte solução,

kk iiky )(21)(

21)(

Faça um gráfico para a solução acima. clear all clc Np=input('de o numero de pontos '); for k=1:Np index(k)=k; i=sqrt(-1); y(k)=0.5*(i)^k+0.5*(-i)^k; end plot(index,y,'*',index,y) xlabel('k') ylabel('y(k)')

13

Figura 3 – Gráfico da resposta do sistema.

Exemplo 13. A Ruína do Jogo.

Considere dois jogadores A e B, sendo que A é o convidado e B é a banca (o dono do jogo).

Seja p a probabilidade que o jogador A ganhe uma ficha de B. Assim, pq 1 é a probabilidade de B ganhar uma ficha de A.

Sejam a e b as quantidades iniciais de fichas de A e de B respectivamente. Qual a probabilidade do jogador A ganhar todas as fichas de B? Esta pergunta será respondida usando um modelo dinâmico que será apresentado

a seguir. Seja k o número de fichas de A, então, bak 0 . O número de fichas de B é, portanto, kba . Seja )(ku a probabilidade do jogador A, eventualmente, ganhar todas as fichas de B, tendo consigo uma quantidade de fichas k. Seja o modelo abaixo,

)1()1()( kuqkupku com as seguintes condições iniciais 0)0( u (probabilidade de A ganhar todas as fichas de B tendo nenhuma ficha) e 1)( bau (probabilidade de A ganhar todas as fichas de B tendo todas as fichas). Pode-se reescrever o modelo acima como,

0)1()()1( kqukukpu

14

Trata-se de uma equação de diferenças linear, homogênea com coeficientes invariantes no tempo. A equação característica associada ao modelo acima é,

02 qp cujas raízes são 11 e pq2 . A solução da equação dinâmica acima é portanto,

k

pqccku

21)(

Impondo as condições iniciais, pode-se determinar as constantes 1c e 2c através do seguinte conjunto de equações algébricas linear,

1

0

21

21

cpqc

ccba

Que leva a solução do modelo como,

ba

k

pq

pq

ku

1

1)(

Para o jogo de roleta, por exemplo, onde existem 18 divisões vermelhas, 18 pretas e 1

verde. Neste caso 3718

p , 3719

q . Para 100a fichas e 1000b fichas, qual a

probabilidade do jogador A ganhar todas as fichas de B?

241100

100

1029.3

18191

18191

)100(

u .

2.7 Equações Diferenciais Ordinárias

são descritas em termos de tempo contínuo, y(t) representa uma função definida dentro de um intervalo de tempo 10 ttt , equação diferencial ordinária é uma equação que conecta uma função y(t) e

algumas de suas derivadas. Exemplos de EDO’s:

)()( tyadt

tdy (EDO linear de primeira ordem sem termo força)

15

)cos()())(()(2

2

tdt

tdytysendt

tyd (EDO não linear de segunda ordem com

termo força). Qual é o método numérico usado para resolver EDO’s não linear? A ordem de uma EDO é definida pela derivada de maior ordem.

Linearidade Uma EDO linear possui a seguinte forma,

)()()()()()()()()( 011

1

1 tgtytadt

tdytadt

tydtadt

tydta n

n

nn

n

n

(2.18)

A EDO pode possuir, em geral, coeficientes variantes no tempo. Se eles não dependem de t, são chamados invariantes no tempo e )(kg denota o termo força ou simplesmente lado direito. Condições Iniciais É usualmente necessário especificar um conjunto de condições iniciais para especificar uma única solução para uma EDO.

Exemplo: considere a EDO )()( tyadt

tdy . A resposta é a seguinte,

taCety )(

sendo que a constante arbitrária C é determinada pela condição inicial 0)0( yy . Substituindo, para o instante 0t o termo 0)0( yy na forma da solução acima, encontra-se que,

0yC Portanto,

taeyty 0)(

Exemplo: Considere a EDO de segunda ordem 0)(2

dt

tyd , a forma de sua solução é a

seguinte,

BtAty )(

Impondo as condições iniciais encontra-se que )0(yA edt

dyB )0( , portanto,

16

tdt

dyyty )0()0()(

Teorema da Existência e Unicidade É consideravelmente muito mais complicado estabelecer uma prova para a existência e unicidade de uma EDO que para equações de diferenças. Suponha os coeficientes )(tai com ni ,,1,0 e a função g(t) da equação (2.18) sejam contínuas sobre um intervalo Tt 0 . Então, para qualquer conjunto de valores ib , com 1,,1,0 ni , existe uma única solução para a EDO (2.18), satisfazendo as condições iniciais,

11

1

1

0

)0(

)0()0(

nn

n

bdt

yd

bdt

dyby

2.8 Equações Diferenciais Ordinárias Lineares Novamente aqui a equação homogênea desempenha um papel fundamental na solução geral. Equações Homogêneas A equação homogênea é obtida fazendo 0)( tg na Equação (2.18), que dá o seguinte,

0)()()()()()()()( 011

1

1

tytadt

tdytadt

tydtadt

tydta n

n

nn

n

n

Teorema 1: Seja )(ty uma dada solução da EDO (2.18). Então, a coleção de todas as soluções desta equação é a coleção de todas as funções da seguinte forma,

)()()( tztyty onde )(tz é uma solução da equação homogênea correspondente. Teorema 2: Se )(,),(),( 21 tztztz m são m soluções da equação diferencial homogênea, então qualquer combinação linear destas m soluções,

)()()()( 2211 tzctzctzctz mm é também é uma solução, onde mccc ,,, 21 são constantes arbitrárias.

17

Exemplo 14: Equação de primeira ordem. Seja baydtdy

uma EDO,

Por inspeção, abty )( é uma solução, a homogênea correspondente é a

seguinte EDO,

aydtdy

cuja solução possui a seguinte forma,

atCetz )( sendo C uma constante arbitrária. Portanto, a solução geral para esta EDO proposta acima é a soma das soluções particular e a homogênea como,

abCety at )(

Soluções Fundamentais O conjunto de n soluções homogêneas de (2.18) )(,),(),( 21 tztztz n obtidas a partir das condições iniciais,

.,0

1,1)0(contcasoikse

dtzd k

Define-se um conjunto de soluções fundamentais o seguinte,

0)0(,0)0(,0)0(,1)0(12

1

dtzd

dtzd

dtdzz

n

0)0(,0)0(,1)0(,0)0(12

1

dtzd

dtzd

dtdzz

n

1)0(,0)0(,0)0(,0)0(12

1

dtzd

dtzd

dtdzz

n

Teorema 3: Se )(tz é qualquer solução da equação homogênea de (2.18), então )(tz pode ser expressa em termos de n soluções fundamentais na forma,

)()()()( 2211 tzctzctzctz nn onde mccc ,,, 21 são constantes arbitrárias.

18

Teorema 4: Suponha que )(,),(),( 21 tztztz m é um conjunto de soluções da equação homogênea associada com (2.18) linearmente independentes. Então, qualquer solução z(t) pode ser expressa como uma combinação linear,

)()()()( 2211 tzctzctzctz nn para algumas constantes arbitrárias mccc ,,, 21 . Equação Característica Para uma EDO com coeficientes constantes, a equação homogênea correspondente pode ser resolvida considerando uma equação característica. Este método é baseado em assumir que as soluções da forma tetz )( existem para alguma constante a ser determinada. Especificamente, considere a EDO homogênea,

0)()()()(011

1

1

tyadt

tdyadt

tydadt

tyda n

n

nn

n

n

Supondo como solução tetz )( e substituindo na equação acima, obtém-se o seguinte,

0011

1

tttnn

tnn eaeaeaea

Como te nunca é nulo, portanto,

0011

1 aaaa n

nn

n (2.19) A equação (2.19) acima é chamada polinômio característico. As raízes do polinômio característico são chamadas valores característicos. Se as raízes do polinômio característico são distintas, n soluções diferentes são obtidas para a solução homogênea. Nestas notas será considerado apenas o caso onde a equação (2.19) apresenta raízes distintas. Exemplo 15: EDO de primeira ordem. Considere a seguinte EDO

aydtdy

0 a (polinômio característico)

a (raiz do polinômio característico) taCety )( (solução da EDO)

Exemplo 16: EDO de segunda ordem com forçamento. Sistema mecânico tipo massa, mola amortecedor com forçamento harmônico.

19

Vai-se analisar o comportamento dinâmico de um sistema massa mola amortecedor de um único grau de liberdade conforme segundo o modelo físico mostrado na Figura (4).

A coordenada y(t) denota posição do bloco de massa m a partir da posição de equilíbrio estático.

A coordenada y(t), suas primeira e segunda derivadas dt

tdy )( e 2

2 )(dt

tyd são

assumidas positivas quando apontas para baixo. Além disso, os termos y(t) dt

tdy )( e

2

2 )(dt

tyd serão denotados, por simplicidade, respectivamente, por yeyy , .

O termo força será adotado como )cos()( 0 tFtF , onde 0F e denotam,

respectivamente, amplitude e freqüência (em rad/s) da força externa ou excitação do sistema mecânico.

Figura 4 – Modelo físico e diagramas de corpo livre e força resultante.

Assumindo que a força elástica da mola varie linearmente com o deslocamento y e que a força de amortecimento varie também linearmente com a velocidade y , com base na segunda lei de Newton, comprando o diagrama de corpo livre com o diagrama de forças resultantes (não mostrado na Figura 4), leva ao seguinte resultado,

ymtFyckyymFext )( Equação do movimento Da segunda lei de Newton obtém-se a equação do movimento abaixo,

)cos(0 tFykycym (2.20)

20

Trata-se de uma EDO de segunda ordem linear com coeficientes constantes com forçamento harmônico cuja solução da homogênea será determinada usando o método da equação característica. Solução da EDO homogênea associada ou vibração livre (sem forçamento) A solução da homogênea (ou movimento livre de forçamento externo) é analisada e possui a seguinte forma,

0 ykycym (2.21) cuja solução leva ao polinômio característico da seguinte forma,

02 kcm (2.22) As raízes do polinômio característico são as seguintes,

mk

mc

mc

mk

mc

mc

2

1

2

1

22

22

(2.23)

que levam a solução homogênea da seguinte forma,

tmk

mc

mct

mk

mc

mc

tt

tmk

mc

mct

mk

mc

mc

tt

ecececectz

e

ecececectz

22

21

22

21

22

22

22

112211

22

2

22

121

)(

)(

(2.24)

onde 1c e 2c são constantes arbitrárias a serem determinadas pelas condições iniciais impostas ao movimento. Definições: Freqüência natural, constante de amortecimento crítico, fator de amortecimento e freqüência natural amortecida A fim de simplificar a representação da solução homogênea e da particular em seguida são definidas as seguintes quantidades:

mk

n (freqüência natural em rad/s)

21

cc : constante de amortecimento crítico, definida como o valor da constante de amortecimento para qual o radical da equação (2.23) torna-se nulo, isto é,

02

2

mk

mcc ou nc m

mkmc 22

ccc

(fator de amortecimento)

Fazendo o termo mc

2 da equação (2.23) como n

c

c mc

cc

mc

22

, as

expressões de 1 e 2 tornam-se o seguinte,

n

n

mk

mc

mc

mk

mc

mc

122

122

22

1

22

1

(2.25)

Substituindo a equação (2.25) na equação (2.24), obtém-se uma forma mais

conveniente e mais simplificada para a resposta da homogênea como,

tttt nn ecececectz 12

1121

2221)(

Para 10 , o termo dn 21 é chamado de freqüência natural amortecida. Com base nos valores assumidos pelo fator de amortecimento são analisados quatro casos para a solução homogênea (válidos também, é claro, para o caso de vibração livre), Caso 1: Quando 0 Neste caso tem-se o caso de vibração livre não amortecidas. As expressões das raízes do polinômio característico, neste caso, são,

nimk

)(2,1 (2.26)

sendo 1i .

Agora impondo as condições iniciais 0)0( yy e 0)0( yy , determina-se a solução para as constantes 1c e 2c . Para tanto é necessário a determinação das expressões de )(tz e )(tz que são as seguintes,

22

tt

tt

ecectz

ecectz

21

21

2211

21

)(

)(

(2.27)

Impondo as condições iniciais 0)0( yy e 0)0( yy na equação (2.27) acima, o

calculo de 1c e 2c é obtido através da solução do seguinte sistema de equações lineares,

0

0

2121

1

2121

2

0

01

21̈2

1

0

0

2

1

21¨

1

111

11

yy

yy

cc

ouyy

cc

donde obtém-se imediatamente as expressões para as constantes 1c e 2c como,

21

00

21

1

21

00

21

2

2

1

yy

yy

cc

(2.28)

Substituindo as expressões de 1 e 2 da equação (2.26) na equação (2.28) acima, encontram-se as expressões de 1c e 2c como,

00

00

00

00

2

1

21

21

21

21

yi

y

yi

y

iiyy

iii

iiyy

iii

cc

n

n

nnnn

n

nnnn

n

(2.29)

Agora substituindo as expressões de 1 e 2 da equação (2.26) e de 1c e 2c da

equação (2.29) acima na Equação (2.27) obtém-se,

23

titititin

n

titi

n

titi

ti

nn

ti

nn

ti

n

ti

n

tt

nnnn

nnnn

nn

nn

eeyeeiyitz

eeyi

eeytz

ou

eyi

yieyi

yitz

eyi

yeyi

yecectz

00

00

0000

000021

21

2)(

21

21)(

21

21

21

21)(

21

21

21

21)( 21

(2.30)

Lembrando que,

)sin()cos( tite nnti n

e substituindo na equação (2.30) obtém-se, finalmente, as expressões de z(t) e

)(tz como,

)()cos()( 00 tsenytytz n

nn

(2.31)

e

)cos()()( 00 tytsenytz nnn (2.32) Das equações (2.30) e (2.31) acima nota-se que para o caso de vibração livre sem amortecimento o bloco do sistema vibra na freqüência natural n . A figura (5) mostra a resposta de um sistema massa mola livre com massa igual a 1 kg e rigidez igual a 1000 N/m para condições iniciais 0)0( y e 1)0( y m/s. O intervalo de discretização usado foi de 0005.0t seg.

24

Figura (5) – Resposta de um sistema massa mola livre.

clear all clc m=1; % massa do bloco em kilogramas k=1000; % rigidez da mola em N/m wn=sqrt(k/m); % calculo da frequencia natural Np=2000; % numero de pontos dt=0.0005; % intervalo de tempo em segundos y0=0; % condiçao inicial de posiçao yp0=1; % condiçao inicial de velocidade for i=1:Np t(i)=dt*(i-1); z(i)=y0*cos(wn*t(i))+(yp0/wn)*sin(wn*t(i)); zp(i)=-wn*y0*sin(wn*t(i))+yp0*cos(wn*t(i)); end figure(1) plot(t,z,t,zp) xlabel('seg') ylabel('z(t) e zp(t)') legend('z(t)','zp(t)') Caso 2: Sistema subamortecido, quando ccc ou 10 Neste caso as raízes do polinômio característico são as seguintes,

25

n

n

i

ei

22

21

1

1

Impondo as condições iniciais 0)0( yy e 0)0( yy e usando procedimento análogo ao do caso anterior, mostre que a resposta livre do sistema subamortecido, é dada pela seguinte equação,

tsenyytyetz

ou

tsenyytyetz

dd

nd

t

nn

nn

t

n

n

000

2

2002

0

cos)(

11

1cos)(

(2.33)

Dá para mostrar também que a Equação (2.33) acima pode ser escrita em uma representação mais compacta em termos de amplitude e fase como,

tsenYetz dtn)( (2.34)

onde a amplitude Y e a fase são calculadas como,

d

n

d

n

yyytg

e

yyyY

0

001

2

0020

Com base nas equações (2.33) ou (2.34), nota-se que a resposta livre de um

sistema massa mola amortecedor subamortecida ocorre, portanto, na freqüência natural amortecida e apresenta um decaimento exponencial.

A figura (6) mostra a resposta de um sistema massa mola subamortecido livre com massa igual a 1 kg, constante de amortecimento igual a 20N/ms e rigidez igual a 1000 N/m para condições iniciais 0)0( y e 1)0( y m/s. O intervalo de discretização usado foi de 0005.0t seg.

26

Figura (6) – Resposta de um sistema massa mola subamortecido livre.

clear all clc Np=2000; % numero de pontos m=1; c=20; k=1000; wn=sqrt(k/m) cc=2*m*wn; qsi=c/cc wd=(1-qsi^2)*wn y0=0; yp0=1; dt=0.0005; % intervalo de discretizaçao for k=1:Np t(k)=(k-1)*dt; z(k)=exp(-qsi*wn*t(k))*(y0*cos(wd*t(k))+((yp0+qsi*wn*y0)/wd)*sin(wd*t(k))); end plot(t,z) xlabel('seg') ylabel('z(k)') Caso 3: Sistema criticamente amortecido, quando ccc ou 1 Neste caso as duas raízes do polinômio característico são iguais como,

27

nc

mc 221 (2.35)

Mostre que para condições iniciais do tipo 0)0( yy e 0)0( yy , a expressão para a resposta livre do sistema criticamente amortecido é a seguinte,

tn

netyyytz 000)( (2.36) A equação (3.36) mostra que tal sistema apresenta um movimento aperiódico, isto é, não vibra como também mostra a figura (8).

Figura (6) – Resposta de um sistema massa mola criticamente amortecido livre.

clear all clc Np=2000; % numero de pontos m=1; k=1000; wn=sqrt(k/m) qsi=1 wd=(1-qsi^2)*wn y0=0; yp0=1; dt=0.0005; % intervalo de discretizaçao for k=1:Np t(k)=(k-1)*dt; z(k)=exp(-wn*t(k))*(y0+(yp0+wn*y0)*t(k)); end

28

plot(t,y) xlabel('seg') ylabel('z(k)') Caso 4: Sistema superamortecido, quando ccc ou 1 Neste caso, as raízes do polinômio característico são reais de distintas como,

n

n

e

1

1

22

21

(2.37)

Usando as equações (2.24) e impondo as condições iniciais 0)0( yy e

0)0( yy , e adotando procedimento análogo ao usado para os casos do sistema sem amortecimento e subamortecido, mostre que a resposta homogênea, neste caso, é a seguinte,

t

n

nt

n

n eyy

eyy

tz 21

12

1

12

1)(

2

02

0

2

02

0

(2.38)

Figura (7) - Resposta de um sistema massa mola superamortecido livre.

clear all clc Np=2000; % numero de pontos

29

m=input('de o valor de m '); k=input('de o valor de k '); wn=sqrt(k/m); c=input('de o valor de c '); cc=2*m*wn; qsi=c/cc wd=(1-qsi^2)*wn; lambda1=(-qsi+sqrt(qsi^2-1))*wn lambda2=(-qsi-sqrt(qsi^2-1))*wn y0=input('de o valor de y0 '); yp0=input('de o valor de yp0 '); dt=0.0005; % intervalo de discretizaçao A=[ 1 1; lambda1 lambda2]; b=[y0 yp0]'; C=inv(A)*b; for k=1:Np t(k)=(k-1)*dt; D=[ exp(lambda1*t(k)) exp(lambda2*t(k)); lambda1*exp(lambda1*t(k)) lambda2*exp(lambda2*t(k))]; Z=D*C; z(k)=real(Z(1)); zp(k)=real(Z(2)); end plot(t,z,t,zp) xlabel('seg') ylabel('z(t),zp(t)') legend('z(t)','zp(t)')

O código computacional acima é bem geral e vale para o calculo da resposta homogênea para qualquer tipo de sistema massa, mola amortecedor para qualquer fator de amortecimento, exceto para o caso criticamente amortecido. Os parâmetros adotados para obtenção do movimento mostrado na Figura (7) são mostrados abaixo quando da execução do código acima. de o valor de m 1 de o valor de k 1000 de o valor de c 100 qsi = 1.5811 lambda1 = -11.2702 lambda2 = -88.7298 de o valor de y0 0.25 de o valor de yp0 1

Solução particular

30

A forma da solução particular é adotada como,

)cos()( pp tYty (2.39) onde que pY e p são constantes a serem determinas, sendo, respectivamente, a amplitude e a fase de )(ty . Substituindo a equação (3.39) com suas derivadas de primeira e segunda ordens na equação (2.20), obtém-se o seguinte,

)cos(cos 02 tFtsenctmkY ppp (2.40)

Usando as relações trigonométricas a seguir,

)()cos()cos()(

)()()cos()cos(cos

ppp

ppp

senttsentsen

sentsentt

e substituindo em (3.40), e igualando os coeficientes de )cos( t e )( tsen em ambos os lados dá o seguinte,

0)cos()(

)()cos(

2

02

ppp

ppp

csenmkY

FsencmkY

(2.41)

As equações acima representam um sistema de duas equações algébricas (não lineares) para duas incógnitas pY e p . Somando o quadrado das partes superior e inferior da Equação (3.41) encontra-se,

21

2222

0

cmk

FYp

(2.42)

e usando a equação de baixo, determina-se a expressão da fase como,

21

mk

ctgp (2.43)

A fim de reescrever as Equações (2.42) e (2.43) de forma mais compacta, define-se deflexão estático devido 0F como,

kF

st0

e razão de freqüências como,

n

r

31

Usando os termos acima, as Equações (2.42) e (2.43) podem ser reescritas como,

21

222

21

1

nn

st

pY

(2.44)

2

1

12

rrtgp (2.45)

Segue abaixo, na Figura (8) os gráficos do fator de ampliação e da fase para o sistema mecânico em estudo com forçamento harmônico. Os dados usados estão dados abaixo da listagem do programa em Matlab.

Figura 8 – Gráficos de fator de ampliação e fase

clear all clc m=1; c=input('de o valor de c '); k=1000; disp('frequencia natural') wn=sqrt(k/m) cc=2*m*wn; disp('fator de amortecimento') qsi=c/cc disp('frequencia natural amortecida') wd=sqrt(1-qsi^2)*wn r=input('de o valor maximo para r ');

32

Np=input('de o numero de pontos para plotar a fase em funcao de r '); dtr=r/Np; for i=1:Np r(i)=(i-1)*dtr; fa(i)=1/(sqrt((1-r(i)^2)^2+(2*qsi*r(i))^2)); y=2*qsi*r(i); x=1-r(i); fi(i)=atan2(y,x); end subplot(2,1,1), plot(r,fa) xlabel('r') ylabel('fator de ampliaçao') grid subplot(2,1,2), plot(r,fi*180/pi) xlabel('r') ylabel('fase') grid O código acima foi rodado usando os parâmetros conforme lista de comandos abaixo, de o valor de c 10 frequencia natural wn = 31.6228 fator de amortecimento qsi = 0.1581 frequencia natural amortecida wd = 31.2250 de o valor maximo para r 4 de o numero de pontos para plotar a fase em funcao de r 2000

Observe das curvas do fator de ampliação st

pY

e da fase p o seguinte:

o valor máximo de st

pY

ocorre quando 221 r ou,

equivalentemente, quando 221 n .

O valor máximo de st

pY

(quando 221 r ) é dado por

2121

st

pY

Para qualquer valor de implica que 090)1( rp .

33

Mais detalhes interessantes sobre as curvas de fator de ampliação (ou ganho) e fase podem ser encontrados na literatura de vibrações mecânicas.

Resposta total De maneira geral, a solução total para a posição do bloco a partir da posição de

equilíbrio estático y(t) é dada por,

)()()()( 2121 tyecectytzty tt (2.46)

Para calcular o valor das constantes 1c e 2c arbitrárias é necessário neste agora impor as duas condições iniciais de posição e velocidade 0)0( yy e 0)0( yy . Para tanto, a expressão da velocidade do bloco é dado por,

)()()()( 212211 tyecectytzty tt (2.47)

A determinação das 1c e 2c é obtida resolvendo simultaneamente as equações (2.46) e (2.47) para o instante 0t que resulta no seguinte sistema de duas equações algébricas lineares como,

)0()0(11

)0()0(11

)0()0(

0

01

212

1

0

0

2

1

2102211

021

yyyy

cc

yyyy

cc

yyccyycc

(2.48)

de onde são calculados os valores de 1c e 2c para serem substituídos em (2.46) e (2.47). Para este caso especial de movimento harmônico, onde )cos()( pp tYty , tem-se o seguinte,

)()(

)cos()(

pp

pp

tsenYty

tYty

)()cos(11

)()cos(11

0

01

210

01

212

1

pp

pp

pp

pp

senYyYy

senYyYy

cc

A Figura (9) mostra os gráficos do movimento harmonicamente forçado para o caso de um sistema subamortecido.

Note que a resposta da homogênea ocorre na frequência natural amortecida, a resposta da particular na frequência da força harmônica externa que excita o sistema e a resposta total ocorre na frequência natural amortecida e na frequência de excitação. Entretanto, a resposta homogênea é transitória, devido seu decaimento exponencial.

34

Após o término dessa fase transitória da homogênea, a resposta total ocorre apenas na frequência de excitação.

A listagem do programa com parâmetros usados são apresentados no que segue.

Figura (9) – Vibração harmonicamente forçada.

clear all clc Np=2000; % numero de pontos m=1; % massa do bloco k=1000; % igidez da mola wn=sqrt(k/m) % calculo da frequencia natural c=5; % constante de amortecimento cc=2*m*wn; % calculo do amortecimento critico qsi=c/cc % calculo do fator de amortecimento wd=(1-qsi^2)*wn % calculo da frequencia natural amortecida lambda1=(-qsi+sqrt(qsi^2-1))*wn; % calculo da raiz do polinomio caracteristico lambda2=(-qsi-sqrt(qsi^2-1))*wn; % " y0=0; % condiçao inicial de posiçao yp0=1; % condiçao inicial de posiçao F0=10; % amplitude da força externa w=20 % frequencia da força externa dt=0.001; % intervalo de discretizaçao r=w/wn; % razao de frequencias deltast=F0/k; % deflexao estatica devido F0 fa=1/(sqrt((1-r^2)^2+(2*qsi*r)^2)); % calculo do fator de ampliaçao

35

Yp=fa*deltast; % calculo da amplitude da particular fy=2*qsi*r; fx=1-r; fs=atan2(fy,fx); % calculo da fase A=[ 1 1; lambda1 lambda2]; b=[y0-Yp*cos(fs) yp0-Yp*sin(fs)]'; C=inv(A)*b; for k=1:Np t(k)=(k-1)*dt; D=[ exp(lambda1*t(k)) exp(lambda2*t(k)); lambda1*exp(lambda1*t(k)) lambda2*exp(lambda2*t(k))]; Z=D*C; z(k)=real(Z(1)); % soluçao da homogenea zp(k)=real(Z(2)); % solucao da derivada temporal da homogenea yp(k)=Yp*cos(w*t(k)-fs); % soluçao da particular ypp(k)=-Yp*w*sin(w*t(k)-fs); % soluçao da derivada temporal da particular yt(k)=z(k)+yp(k); % soluçao total ytp(k)=zp(k)+ypp(k); % soluçao da derivada temporal da soluçao total end subplot(3,1,1), plot(t,z) xlabel('seg') ylabel('z(t)') subplot(3,1,2), plot(t,yp) xlabel('seg') ylabel('yp(t)') subplot(3,1,3), plot(t,yt) xlabel('seg') ylabel('yt(t)') Exemplo 17 – Resolva de maneira análoga ao exercício anterior a EDO

)(0 tsenFykycym .

Exemplo 18 – Resolva a EDO teconsFykycym tan0 para condições iniciais de posição e velocidade. Faça gráficos dos resultados conforme apresentado no Exemplo 16. 2.9 Série de Fourier

Apenas funções reais serão tratadas aqui. Sugere-se fortemente recorrer à literatura de matemática para mais detalhes

relacionados com a parte teórica.

36

Uma função x(t) periódica com período T pode ser expressa através de um conjunto de funções exponenciais )(tk como,

k

kk txtx )()( (2.49)

onde,

)()cos()( 000 tksenitktkietk sendo,

T 2

0 (frequência fundamental em rad/s)

2/

2/

,2,1,0,1,2,)( 0

T

T

tkik kparadtetxx (coeficientes de Fourier)

A série de Fourier relaciona, portanto, uma função contínua e periódica no

domínio do tempo com funções discretas no domínio da frequência.

Funções Reais Para funções reais, a série de Fourier pode ser escrita como uma série

trigonométrica, como,

1

0 2sin2cos2

)(k

kk tTkbt

Tkaatx (2.50)

onde, os coeficientes da série de Fourier são calculados como,

T

txT

a0

0 )(2

,2,12cos)(2

0

kdtt

Tktx

Ta

T

k

,2,12)(2

0

kdtt

Tksentx

Tb

T

k

Exemplo 19 - Calcule a expressão da função onda quadrada de amplitude igual a 1 e período de 2 segundos mostrada na Figura (10) abaixo em série de Fourier.

37

Figura (10) - Onda quadrada

,2,1,00)()2(1

22

22cos122cos02

2/2/

2/

0

kksenksenk

a

tTksen

kdtt

Tk

Tdtt

Tk

Ta

k

T

T

T

T

T

k

Da mesma forma, calculam-se os coeficientes kb ,

2,1,0)cos()2cos(12cos2

222

2/2/

kkk

kt

Tk

kdtt

Tksen

Tb

T

T

T

Tk

kbk

2 , para k ímpar

0kb , para k par

38

Figura (11) – Coeficientes ak e bk

Figura (12) – Expansão da função onda quadrada em série de Fourier

% Exemplo numerico - serie de Fourier clear all % 1 - O sinal x(t) original

39

T=2; % periodo do sinal Np=input('de o numero de pontos do sinal '); dt=T/Np; % intervalo de discretizaçao t=(0:1:Np-1)*dt; for i=1:Np/2 x(i)=0; x(i+Np/2)=1; end % 2 - O sinal x(t) obtido via serie de Fourier nt=input('de o numero de termos da expansao da serie de Fourier '); % 2.1 - Coeficientes a(k) e b(k) a0=1; for k=1:nt a(k)=(sin(2*pi*k)-sin(pi*k))/(pi*k); b(k)=-(cos(2*pi*k)-cos(pi*k))/(pi*k); end % 2.2 - Determinaçao de x(t) por expansao em serie de Fourier for i=1:Np sum=0; for j=1:nt sum=sum+a(j)*cos(2*j*pi*t(i)/T)+b(j)*sin(2*j*pi*t(i)/T); end xf(i)=a0/2+sum; end % 3 - Figuras w0=2*pi/T; w=(1:1:nt)*w0; figure(1) subplot(2,1,1), plot(0,a0,'+',w,a,'+') xlabel('rad/s') ylabel('a(w)') grid subplot(2,1,2), plot(w,b,'+') xlabel('rad/s') ylabel('b(w)') title('Coeficientes de fourier') grid figure(2) plot(t,x,t,xf) xlabel('seg.') ylabel('x(t)') grid >> fourier1 de o numero de pontos do sinal 1000 de o numero de termos da expansao da serie de Fourier 5

40

Exemplos: calcule a expansão em série de Fourier para as funções abaixo Figura – dente de serra

Figura – seno truncado

210

10)()(

tparatparatsen

tx

2.10 Solução de EDO’s para Forçamento Periódico

O uso da série de Fourier para expansão da função periódica é a chave para a solução do problema

Para sistemas lineares usa-se o principio da superposição, solução particular é composta pela soma de dois efeitos: resposta a uma força constante e resposta a várias forças harmônicas (do tipo senoidais e cossenoidais).

Truncamento da série de Fourier para uma quantidade finita de componentes espectrais introduz erro.